Abnormally Distributed

統計解析担当のライフサイエンス研究者 -> データサイエンティスト@コンサル

Coursera:AI for Medicine Specializationの受講記録③

随分間が開いてしまったが、AI for Medicine Specializationの3つ目のコースも受講完了したので、レビューしておく。

Course3. AI for Medical Treatment

治療のためのAIというタイトルで、治療効果の分析の他、テキストデータの分析や機械学習モデルの解釈方法がメインテーマとなっている。 キーワードはこんな感じ。

  • RCTデータを用いた効果検証
  • BERTによる質問応答
  • SHAP値
  • GradCam

Week1

ランダム化比較試験(RCT)のデータを用いた治療効果推定がテーマ。因果推論の基本が分かっていれば、特に目新しい内容はなさそう。
ルービンの因果モデルでは平均治療効果は、ある個人について治療した場合と治療しなかった場合の結果の差を全被験者で平均したものだが、実際は個人について両方の結果を観測することはできないため、個人レベルでの治療効果の推定は不可能としている(因果推論の根本問題)。ただし、被験者が治療群とコントロール群にランダムに割り当てられる状況であれば(RCT)、介入の有無以外の変数については両群で平均的には差がないとみなせるため、両群での結果の平均値の差分を取ることで、平均治療効果(ATE, Average Treatment Effect)が推定できる。講義ではこのあたりが詳しく説明されている。

自分が知らなかった内容は、T-learner, S-learnerに関する部分。これらは条件付き治療効果の推定のための手法。例えば年齢が58歳、血圧が130の人について、ATEがどの程度になるか?という条件付き治療効果の推定の問題を考える。被験者数が限られることから、両群から条件に合う患者を抽出して、通常の治療効果の推定のように平均的な効果の差分を取るアプローチは難しい。そこで、共変量から介入あり、なしの場合のアウトカムを予測するモデルを構築し、予測値の差分として条件付き治療効果が推定できる。T-learnerは治療群と対称群で別々のモデルを構築、S-leanerは治療の有無も説明変数に加えて1つのモデルを構築するアプローチ。予測モデルとしては、決定木モデルを使っている。

また、上記の治療効果推定モデルの性能評価についても触れられている。ここでも、同じ患者で治療有無での結果がないことから、評価方法は工夫する必要がある。ちょっとややこしいが、あまり聞いたことがない方法だったので、参考になった。
具体的には、治療群・対称群から、治療効果が同程度と推定された被験者のペアを作成、ペアのアウトカム(治療により悪化、改善、変化なし)と治療効果の推定値を算出し、全ペアについて、アウトカムと推定値の辻褄がどの程度合っているかにより評価する。以前に出てきたC-indexに似た指標である、C-for-benefitを評価指標としている。

Week2

自然言語処理がテーマで、臨床報告からの疾患ラベルの抽出と、質問応答が取り上げられている。
質問応答は質問文と文章が与えられ、回答に対応する部分を文章から抜き出す問題。BERTの大まかな概念や事前学習タスクについて説明がある。Transformerの詳細などには触れていないが、取り敢えずBERTを動かしたく、モデルの概要を知りたいという人にはわかりやすくて良い説明と思われる。
X線画像に対し、各種疾患の有無のラベルデータを作成したいとする。臨床報告から情報を抽出したいが、書式が統一されておらず、一筋縄ではいかない。課題の1つとして用語の包含関係や同義語が挙げられる。これに対処するため、国際医療用語集(NOMED-CT, Systematized Nomenclature of. Medicine-Clinical Term)という用語集を利用する例を紹介している。
また否定表現の抽出も課題である。疾患が言及されていても、否定文かどうかを正しく判定できないと、間違ったラベルを与えてしまう。 講義ではNegBioという臨床テキストデータからの否定表現の抽出用のライブラリを紹介している。

Week3

機械学習モデルの解釈がテーマであり、SHAPなどの特徴量重要度とGradCamが取り上げられている。
まずは特徴量重要度の算出方法として、①特徴量を削除した際の性能の変化に基づく方法、②特徴量の値をランダムにシャッフルした際の性能の変化に基づく方法(Permutation importance)) が紹介されている。後者の方法は、モデルの再学習が必要にならないのが利点である。個人的には以前からちょっと怪しい方法と思っていたが、このコースで紹介されるくらいなので、ある程度確立された方法なのかもしれない。。
次に、個々のレコードレベルでの特徴量重要度を算出する方法として、シャープレイ値(Shapley values)を取り上げている。自分としては、このSHAPに関する解説が一番勉強になった。
Shapley Valueは元々、ゲーム理論で提唱された概念とのこと。具体的な算出方法は以下の通り。まず全ての特徴量の順列を作成する。次に、特徴量を左から順番に、評価対象の特徴量に達するまで取っていく。その上で、対象の特徴量を含めた場合と含めなかった場合のモデルの出力の差分を取り、全ての順列についての平均を取ると、対象の特徴量に関するShapley valueとなる。全特徴量から対象の特徴量を除外した場合だけでの評価としないのは、特徴量どうしの相関の影響をキャンセルするため。
Shapley valueの便利な性質として、全特徴量のShapley valueの合計が、その個人におけるモデルの出力とベースの出力(特徴量なしの場合)との差分と一致することが挙げられる。要は、Shapley valueは出力値に各特徴量がどの程度貢献しているかを表すことになる。これが、よく見る SHAPのforce plotということになる。
なお、Pythonの SHAPパッケージでは実際に全ての特徴量セットでモデルを再学習している訳ではなく、高速に計算できるアルゴリズムを利用している。

次に、CNNのモデルを解釈する手法として、GradCamを説明している。Course1でも出てきたので、詳しい説明は割愛する。要は(flattenされる直前の)最後の畳み込み層の出力(特徴マップ)がタスクに関連する重要な空間的情報を捉えていると考え、k枚の特徴マップについて、モデルの出力値に対する影響に比例する重み(勾配の平均的な大きさ)による加重平均を取った上で、入力画像と同じサイズにリサイズしてヒートマップを作成し、元の画像と重ね合わせることで、モデルが注視している領域を可視化する方法と言える。

Specializationの総括

統計解析、機械学習の経験がそれなりにある身としては、よく知っている内容も多かったが、新しい学びも多かった。各手法の考え方を、数学的な詳細に立ち入りすぎずにわかりやすく説明しており、ノートブックでの演習を通じて身に付けられるようになっており、課金する価値は十分にあったと思う。

Coursera:AI for Medicine Specializationの受講記録②

AI for Medicine Specializationの2つ目のコースも受講完了したのでレビューしてみる。

Course2. AI for Medical Prognosis

Prognosisは日本語だと「予後」であり、病気の今後の見通しといった意味になる。予後を予測することで、病気のリスクが分かり、予防医療が可能になる。
まずコース全体で取り上げられるキーワードを羅列してみる。

  • 線形モデルランダムフォレスト
  • 生存関数・ハザード関数
  • カプラン・マイヤー法
  • Cox比例ハザードモデル
  • Survival forest

Week1

患者の属性、生活習慣、病歴や検診結果などを説明変数として、線形モデルによりリスクスコアを算出しますという話。 例えば、心房細動の患者が1年以内に脳卒中を発症するリスクを示すCHA2DS2-VASc scoreといった指標が紹介されている。 CHA2DS2–VASc score - Wikipedia

スコア具体的な確率などを表すものではなく、相対的に比較する程度の用途のよう。 一般的に使われているくらいなので、それなりにエビデンスがあるのだろうが、説明変数の線形和でリスクを求めるというのは結構単純だなと感じてしまう。

あとは交互作用 (interaction) や、C-indexという評価指標について解説がある。 C-indexはモデルが出力するリスクと実際の結果が一致している(リスクが高い人の方が、病気を発症している)割合を示す指標。

Week2

Week2では、決定木、ランダムフォレストを用いて予後予測モデルを構築する。
ランダムフォレストの仕組みを理解している人にとっては、モデリングの観点では特に真新しい内容はなさそう。

もう1つ、実際の分析の際にも有用なテーマとして、欠損値の話がある。
欠損がランダムではない場合、欠損があるデータを単純に除外して学習しまうと、バイアスのあるモデルになる可能性があるので注意が必要。

例えば、モデル構築に用いたデータでは、40歳以下の患者で血圧の測定値がほとんど欠損しているとする。ここで、単純に欠損値を含むデータを除外して学習したとする。 一方、新しいテストデータでは、血圧の欠損値がないとする。すると、モデル構築時に用いたデータと新しいテストデータでは年齢の分布が異なることになってしまい、新しいテストデータでは予測性能が悪化してしまう可能性がある。

講義では、欠損値の3つのパターンを説明している。 1. Missing completely at random 患者の属性等に関係なく、完全にランダムに欠損する場合。この場合、単純に欠損のあるレコードを除外してもバイアスは発生しない。

  1. Missing at random" 例えば患者が40歳以下の場合のみ、一定の確率で欠損するような場合。欠損のあるレコードを除外すると、バイアスが発生してしまう。

  2. Missing not at random 観測されていない変数(例えば検査時に機器が故障していたかどうか)に基づき、欠損するかどうかが決まる場合。
    この場合、欠損のあるレコードを除外しても、観測された変数の分布が変わるとは限らない。

データを確認しただけでは、どのカテゴリに属する欠損値か判断することは難しいが、欠損のメカニズムが複数考えられることを認識しておくのが重要だとしている。
また、欠損値補完の話題もあり、平均値補完の他、他の変数との回帰モデルによる補完も紹介されている。欠損値の処理は実務でも注意すべき事項なので、有意義と感じた。

Week3

ここから、いわゆる生存時間解析がテーマになる。個人的には興味があった領域。
Week2までは、固定の期間、例えば1年後に病気を発症するかどうかという二値分類がタスクだったが、生存時間分析では事象の発生までの時間を対象として分析を行う。 Week3ではまず生存関数にいて説明し、生存時間分析に特有の打ち切りデータの取り扱いについて詳しく解説している。 その後、この週の主題となる生存関数の推定量であるカプラン・マイヤー推定量を説明している。

Week4

Week4ではまず、ハザード関数、累積ハザード関数を導入している。 そして患者の属性によるハザードの違いを説明できるCox比例ハザードモデルを解説している。 ただし、個人差を考慮できるといってもハザード関数の形自体は一定のまま、共変量に応じて増減させるだけなので、依然かなり単純なモデルだとは思う。

そこで、非線形の手法として、survival treeというモデルを解説している。Rのパッケージを使ったコードも掲載されているが、プログラミング課題とはなっておらず、参考に紹介している程度。

最後に、生存時間分析モデルの評価方法として、Harrell's C-indexを紹介している。

所感

生存時間分析は少し独特な分野だが、医療だけでなく、マーケティング等の領域でも利用されることがある。
初めて生存時間分析を学ぶ人にも分かりやすいように、丁寧に解説されているように感じたので、医療の分析に限らず生存時間分析に興味がある人にはおすすめ。

ただ、最初の2週間は機械学習を少しやったことがある人なら、いささか簡単すぎる嫌いはある。医療の文脈でのモデリングの紹介としては良いと思うが。。 生存時間分析については、別の記事を書くことにする。

Coursera:AI for Medicine Specializationの受講記録

CourseraでAI for Medicine Specializationを受講してみた。Andrew Ng先生で有名なdeeplearning.aiが提供している。
構成としては、AI for Medical Diagnosis. AI for Medical Prognosis, AI for Medical Treatmentの3つのコースからなり、受講目安時間はそれぞれ19時間、30時間、22時間となっている。 受講の流れは、まず講義ビデオを見た上で、Coruseraが提供する環境上でnotebookを動かしながらデータの理解やモデリングに必要なコーディングを進める。 各週の終わりには選択式のクイズがあり、合格すると評価対象となるnotebookでのプログラミング課題に進める。コードは自動で評価され、すべての週のプログラミング課題で合格できれば修了となる。

授業のビデオは1週あたり30分もあるかどうかで、基本的にはnotebookでの演習に時間を費やすことになる。とはいえ、かなり丁寧な導入もあるので、Pythonでの分析経験がある人の場合、難易度は高くないと思われる。

1. AI for Medical Diagnosis

コンテンツは3週分に分かれている。講義の主題としてはディープラーニングのモデルの解説というよりも、医療データの分析における注意点の解説が中心となっている。

1週目は胸部X線の画像を入力として、複数の疾患の有無を判定するmulti-class classificationが主題。pre-trainedのDensenetによる分類モデルを構築する。
主なテーマは不均衡なデータデットへの対応で、重みつきの損失関数を取り上げている。 また、データセットのサイズに応じた対応として、例えばpre-trained modelを使った転移学習やData augmentationを紹介している。

医療画像分析ならではのトピックとしては、正解ラベルの設定方法があった。複数の医師の多数決や追加検査(CT、バイオプシー)といった方法が考えられるとのこと。

2週目はモデルの様々な評価指標の解説のみで、軽めの内容。 sensitivity, specificity, PPV, NPVなどは確かに医療業界特有の用語なのかもしれない。個人的には特に目新しい内容はなかった。

3週目はMRI画像のセグメンテーションが主題。3D画像には触れたことがなかったので興味深かった。 モデルとしては自分ですら知っているほど定番のU-netを利用。kerasで簡単なU-netの実装を追えるnotebookなどもあり理解が深まったのが良かった。 入力データのサイズが大きいため、ランダムに小さい領域を切り取ってモデルに与えるといった工夫をする必要があるそう。 また、セグメンテーションの評価指標として、Dice coefficientを紹介していた。

全体として、分析のコードは便利なutil関数などが予め用意されており、自分で実装する部分はごく一部だった。自動で採点するシステムなので仕方ないのかもしれないが。 また、コーディング課題もテストケースが通れば良いものではなく、向こうが想定している関数を使わないと誤答になるケースが見られたのはネガティブ要素。 courseraのプラットフォームのスペックがあまり高くないようで、実際にモデルの学習を回せないのも気になる点。コードとデータセットをダウンロードしてcolabなどで動かしやすいようにして欲しかった。

以上のようにコースワークとしてできることは限られているが、参考リンクなどは豊富に提示されているので、この内容をベースに自分でいろいろやってみるのが良さそう。

LDAによるトピックモデル ②gensimによる実装

前回の記事ではLDAの概要や関連手法、確率モデルについて書いた。
今回はPythonのgensimというライブラリを用いて、LDAを実践してみる。
その前に、前回触れていなかったトピックモデルの評価方法について説明する。

評価指標

LDAは教師なしモデルであり、精度等の指標で評価することはできない。
一般的には、PerplexityやCoherenceといった指標が用いられる。

Perplexity

予測性能に関する指標。

定義は下記の通りで、単語ごとの尤度の幾何平均の逆数である。Mは(学習データの)文書数、Ndは文書dの単語数、wdは文書dにおける各単語の出現回数を表す。
教師なし学習では、テストデータの尤度が大きくなるほど、つまりpeplexityが小さくなるほど好ましい。

 \displaystyle perplexity(D _ {test}) = \exp \left(- \frac{\sum _ {d=1} ^ {M} \log p(w _ d)} {\sum _ {d=1} ^ {M} N _ d} \right)

直感的な理解についてはこちらを参照。
トピックモデルの評価指標 Perplexity とは何なのか?

要はPeplexityは選択肢の数を表している。
ある文書の1単語を決める際、何も仮定がない状態ではvocabulary中の全単語が候補になるが、モデルの元では候補となる単語がもっと絞り込めるようになる。この絞り込んだあとの候補数を表すのがPerplexityである。

Coherence

トピックの品質の良さ、人間にとっての解釈しやすさを表す指標。様々な計算方法が提唱されている。 例として、下記に2つのトピックとそのトピックにおいて出現頻度が高い単語を示す。

Topic1: guitar album track bass drum vocal
Topic2: know call name several refer hunter

Topic1は音楽に関するトピックだと明らかだが、Topic2は何を表しているのかはっきりわからない。この場合、Topic1の方がcoherenceが高くなるべきである。

coherenceが高いのは各トピックに含まれる単語の類似度が高い場合である。そこで、coherenceの算出方法としては、各トピックの上位N個の単語を抽出し、その中の各単語ペアについて何らかの類似度を計算して全ペアでの平均値を求めることが多い。
類似度の算出方法については、いくつかの方法が提唱されている。

coherenceの評価基準としては、トピックの良し悪しについて人がつけた正解ラベルとどの程度相関があるかになる。

各手法について、詳しくはこちらを参照。
トピックモデルの評価指標 Coherence 研究まとめ #トピ本

"Exploring the Space of Topic Coherence Measures" M. Röder, A Both, A. Hinneburg, 2015
https://dl.acm.org/doi/10.1145/2684822.2685324
gensimのCoherenceModelはこちらの論文に基づいている。

gensimによる実装

Kaggleのデータセットを用いる。
約20万件のニュース記事の見出しと短い説明文を含む。また、各記事のカテゴリーも記載されており、抽出されたトピックとの関連を確認することもできそう。

コードはKaggleのNotebookとして動かせるようになっている。
LDA_gensim | Kaggle

from time import time
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from IPython.core.display import display
pd.set_option('display.max_colwidth', 200)
%matplotlib inline

# load data
df = pd.read_json('../input/news-category-dataset/News_Category_Dataset_v2.json', lines=True)
df.sample(3)
category headline authors link short_description date
134770 HOME & LIVING Exploiting Familial Ideals for a Car Alissa Fleck, Contributor\nAlissa Fleck is a reporter based in New York City. https://www.huffingtonpost.com/entry/exploiting-familial-ideal-for-car-commercials_us_5b9dc09fe4b03a1dcc8c71ce Is anyone else more than a little crushed when images of gooey, familial perfection splash dazzlingly across the screen to a backdrop of dramatic music swelling to an emotional climax, only to hav... 2014/1/2
5634 POLITICS Why You Shouldn't Freak Out Even If The Stock Market Does Zach Carter https://www.huffingtonpost.com/entry/stock-market-dont-freak-out_us_5a846bc5e4b0774f31d15aaf Everything is fine ... except your paycheck. 2018/2/14
944 ENTERTAINMENT R. Kelly Accusers Say They Had To Ask Permission To Use The Bathroom Sara Boboltz https://www.huffingtonpost.com/entry/r-kelly-accusers-bathroom-permission_us_5af1e468e4b0c4f19327a1f3 And when they did, the women said he demanded they call him "Daddy." 2018/5/8

次に、テキストの前処理をする。
headlineとshort descriptionをつなげたデータを用いる。

from gensim.parsing.preprocessing import preprocess_string

combined_text = df['headline'] + df['short_description']
processed_text = combined_text.map(preprocess_string)

gensimのpreprocess_stringを使うと、テキストデータの一般的な前処理を一通り適用できる。 デフォルトでは下記7つの処理を実行する。

  1. strip_tags(): HTMLタグ等の除去
  2. strip_punctuation(): 句読点の除去
  3. strip_multiple_whitespaces(): 複数の空白の除去
  4. strip_numeric(): 数字の除去
  5. remove_stopwords(): ストップワード(頻出単語)の除去
  6. strip_short(): 短い単語の除去(デフォルトでは3文字未満)
  7. stem_text(): 小文字にして、語幹のみ抽出(porter-stemmerの適用)

実行前
Hugh Grant Marries For The First Time At Age 57The actor and his longtime girlfriend Anna Eberstein tied the knot in a civil ceremony.

実行後
[hugh, grant, marri, time, ag, actor, longtim, girlfriend, anna, eberstein, ti, knot, civil, ceremoni]

Porter stemmerは古典的なステマーだが、age→ag、tied→tiなど元の単語が分かりにくくなるものも多く、正直微妙かもしれない。

次に、各単語に番号を割り振ったディクショナリーを作成する。この際、極端に出現頻度が少ない単語や多い単語は除外する(filter_extremes)。 さらに、ディクショナリーを用いて前処理済みのデータをbag-of-words表現に変換する。

from gensim.corpora.dictionary import Dictionary

# Create a corpus from a list of texts
dictionary = Dictionary(processed_text)
dictionary.filter_extremes(no_below=10, no_above=0.7, keep_n=100000)

corpus = [dictionary.doc2bow(text) for text in processed_text]
corpus[10]
[(7, 1),
 (8, 1),
 (101, 1),
 (102, 1),
 (103, 1),
 (104, 1),
 (105, 1),
 (106, 1),
 (107, 1),
 (108, 1),
 (109, 1),
 (110, 1),
 (111, 1),
 (112, 1)]

このように、文章は(単語id, 出現回数)のリストで表現されている。

from gensim.models import LdaModel, LdaMulticore
lda_model = LdaMulticore(corpus, num_topics=20, id2word=dictionary, random_state=42, workers=3)

print(lda_model.log_perplexity(corpus))
print()

# Get the most significant topics 
for idx, topic in lda_model.print_topics(-1): # -1 for all topics
    print('Topic: {} Word: {}'.format(idx, topic))

log-perprexityは-9.15となった。 また、各トピックにおける単語の出現確率を高い順に出力している。

-9.149718461222365

Topic: 0 Word: 0.013*"tip" + 0.013*"food" + 0.013*"cancer" + 0.011*"life" + 0.010*"sleep" + 0.007*"help" + 0.007*"studi" + 0.006*"research" + 0.006*"need" + 0.006*"risk"
Topic: 1 Word: 0.030*"dai" + 0.024*"parent" + 0.018*"children" + 0.018*"need" + 0.012*"time" + 0.011*"know" + 0.009*"help" + 0.009*"kid" + 0.009*"child" + 0.009*"learn"
Topic: 2 Word: 0.021*"marriag" + 0.011*"anim" + 0.009*"gai" + 0.008*"appl" + 0.007*"kid" + 0.006*"famili" + 0.006*"time" + 0.006*"sex" + 0.006*"peopl" + 0.005*"coupl"
Topic: 3 Word: 0.027*"new" + 0.018*"york" + 0.016*"citi" + 0.015*"yoga" + 0.013*"exercis" + 0.010*"summer" + 0.008*"practic" + 0.006*"spiritu" + 0.005*"chees" + 0.005*"want"
Topic: 4 Word: 0.026*"fashion" + 0.025*"photo" + 0.017*"year" + 0.014*"video" + 0.011*"old" + 0.010*"style" + 0.009*"men" + 0.008*"school" + 0.008*"game" + 0.007*"week"
Topic: 5 Word: 0.027*"year" + 0.011*"babi" + 0.011*"best" + 0.010*"dai" + 0.010*"weight" + 0.010*"mother" + 0.010*"old" + 0.010*"mom" + 0.009*"hotel" + 0.009*"spring"
Topic: 6 Word: 0.013*"peopl" + 0.008*"kill" + 0.008*"polic" + 0.007*"cocktail" + 0.006*"pari" + 0.006*"china" + 0.006*"year" + 0.006*"video" + 0.005*"said" + 0.005*"human"
Topic: 7 Word: 0.023*"obama" + 0.013*"recip" + 0.010*"presid" + 0.009*"flavor" + 0.009*"dish" + 0.008*"trump" + 0.007*"grill" + 0.007*"late" + 0.007*"jimmi" + 0.006*"question"
Topic: 8 Word: 0.018*"red" + 0.015*"morn" + 0.011*"decor" + 0.009*"carpet" + 0.009*"ic" + 0.009*"welcom" + 0.008*"new" + 0.007*"fruit" + 0.007*"cream" + 0.006*"video"
Topic: 9 Word: 0.025*"love" + 0.015*"divorc" + 0.015*"life" + 0.012*"live" + 0.011*"photo" + 0.010*"time" + 0.010*"dai" + 0.010*"thing" + 0.008*"wai" + 0.008*"like"
Topic: 10 Word: 0.039*"photo" + 0.024*"wed" + 0.016*"check" + 0.013*"twitter" + 0.012*"huffpost" + 0.011*"look" + 0.011*"week" + 0.011*"facebook" + 0.011*"style" + 0.010*"video"
Topic: 11 Word: 0.014*"new" + 0.014*"super" + 0.012*"bowl" + 0.010*"workout" + 0.009*"resolut" + 0.008*"state" + 0.007*"america" + 0.006*"media" + 0.006*"polit" + 0.005*"trump"
Topic: 12 Word: 0.022*"studi" + 0.014*"new" + 0.010*"tumblr" + 0.010*"diseas" + 0.008*"peopl" + 0.007*"kate" + 0.007*"heart" + 0.007*"come" + 0.006*"time" + 0.006*"american"
Topic: 13 Word: 0.017*"new" + 0.011*"dai" + 0.010*"like" + 0.010*"video" + 0.008*"dad" + 0.008*"look" + 0.007*"photo" + 0.007*"eat" + 0.007*"gift" + 0.007*"year"
Topic: 14 Word: 0.030*"health" + 0.022*"care" + 0.011*"tax" + 0.008*"patient" + 0.008*"parti" + 0.007*"work" + 0.006*"want" + 0.006*"plan" + 0.006*"help" + 0.006*"republican"
Topic: 15 Word: 0.013*"click" + 0.013*"video" + 0.009*"clinton" + 0.009*"com" + 0.008*"star" + 0.008*"new" + 0.008*"photo" + 0.008*"hillari" + 0.007*"luxuri" + 0.007*"dog"
Topic: 16 Word: 0.014*"world" + 0.012*"healthi" + 0.012*"recip" + 0.011*"food" + 0.010*"wai" + 0.008*"chang" + 0.007*"time" + 0.006*"eat" + 0.006*"wine" + 0.006*"dii"
Topic: 17 Word: 0.031*"photo" + 0.011*"dress" + 0.010*"bodi" + 0.010*"look" + 0.010*"art" + 0.008*"like" + 0.008*"design" + 0.007*"diet" + 0.007*"green" + 0.006*"wai"
Topic: 18 Word: 0.014*"medit" + 0.011*"vintag" + 0.008*"bank" + 0.008*"san" + 0.007*"new" + 0.006*"ebai" + 0.006*"home" + 0.005*"googl" + 0.005*"sandi" + 0.005*"gun"
Topic: 19 Word: 0.016*"pinterest" + 0.015*"new" + 0.009*"water" + 0.007*"loss" + 0.007*"fit" + 0.007*"american" + 0.006*"park" + 0.006*"wall" + 0.006*"drink" + 0.005*"evolut"

WordCloudでトピックを可視化する。

from wordcloud import WordCloud
from PIL import Image
 
fig, axs = plt.subplots(ncols=4, nrows=int(lda_model.num_topics/4), figsize=(15,10))
axs = axs.flatten()
 
def color_func(word, font_size, position, orientation, random_state, font_path):
    return 'darkturquoise'

 
for i, t in enumerate(range(lda_model.num_topics)):
 
    x = dict(lda_model.show_topic(t, 30))
    im = WordCloud(
        background_color='white',
        color_func=color_func,
        random_state=0
    ).generate_from_frequencies(x)
    axs[i].imshow(im)
    axs[i].axis('off')
    axs[i].set_title('Topic '+ str(t))
        
plt.tight_layout()
plt.savefig('/kaggle/working/wordcloud.png')

f:id:kibya:20200517145648p:plain

ある程度は解釈できそうなトピックが見られる。 例えばTopic0はcancer, food, sleep, research, riskなどから、生活習慣と疾病リスクに関する研究結果の記事などが想像される。
またTopic1はparent, children, kid, learnなどから子供や教育に関するトピックだと思われる。

最後に、coherenceを算出する。top_topicsメソッドを用いると、各トピックのcoherenceと出現確率の高い単語を、coherenceの高い順に表示できる。
デフォルトでは、coherenceはu-massと呼ばれる手法で算出する。

topics = lda_model.top_topics(corpus, topn=10)

print('coherence: top words')
for topic in topics:
    print('{:.3f}: {}'.format(topic[1], ' '.join([i[1] for i in topic[0]])))
coherence: top words
-2.436: photo wed check twitter huffpost look week facebook style video
-2.708: love divorc life live photo time dai thing wai like
-2.939: dai parent children need time know help kid child learn
-3.486: new dai like video dad look photo eat gift year
-3.514: tip food cancer life sleep help studi research need risk
-3.548: health care tax patient parti work want plan help republican
-3.575: fashion photo year video old style men school game week
-3.992: photo dress bodi look art like design diet green wai
-4.208: world healthi recip food wai chang time eat wine dii
-4.223: year babi best dai weight mother old mom hotel spring
-4.578: marriag anim gai appl kid famili time sex peopl coupl
-4.654: studi new tumblr diseas peopl kate heart come time american
-4.664: new super bowl workout resolut state america media polit trump
-5.238: peopl kill polic cocktail pari china year video said human
-5.576: click video clinton com star new photo hillari luxuri dog
-6.067: red morn decor carpet ic welcom new fruit cream video
-6.260: new york citi yoga exercis summer practic spiritu chees want
-6.344: obama recip presid flavor dish trump grill late jimmi question
-6.711: pinterest new water loss fit american park wall drink evolut
-8.871: medit vintag bank san new ebai home googl sandi gun

coherenceが上位のトピックは確かに関連する単語が多いが、下位のトピックは一見して無関係の単語が多くなっており、たしかに人間の感覚と一致しているようだ。

LDAによるトピックモデル ①概要と確率モデル

トピックモデルはテキストデータの潜在的意味の解析に用いられる統計モデルの総称である。 ここでは一番有名なモデルである、Latent Dirichlet Allocation (LDA)について説明する。

LDAとは

文書データをLDAで処理することで、文書のトピックの分類や、トピックごとの単語の出現頻度を知ることができる。

LDAでは、単語の羅列である文章の背後に潜在的なトピック(スポーツ、政治、教育など)があると考え、トピックに基づいて文書中の単語が生成されると考える。 例えばスポーツに関する文章なら「試合」「得点」「選手」などの単語が多く、政治に関する文章なら「議会」「法案」「選挙」などの単語が多いと考えられる。 1つ注意点として、LDAは教師なし学習なので、トピックを教師データとして与えて学習する必要はない。

また、LDAでは文章における単語の順序を考慮せず、単語の出現回数のみを考える(Bag of words)。

トピックモデル とLDAの歴史

トピックモデルの起源は1990年に発表されたLatent Semantic Indexing (LSI, LSA, 潜在意味解析)にさかのぼる。
LSIでは共起行列(文書-単語行列)の特異値分解により、トピックを表す潜在変数を取得する手法。例えばK個のトピックを得たい場合は、特異値が大きいものからK個を選択すれば良い(行列の低ランク近似)。
ただし、LSIは確率モデルではなく、文書データの生成プロセスを直接的にモデリングする手法ではない。

1999年にHoffmanらが提唱したProbabilistic Latent Semantic Analysis (PLSA) は確率モデルである。
PLSAは混合モデルであり、文書ごとのトピックの確率 p(z \mid d)を混合比率、トピックごとの単語の生成確率 p(w \mid zを混合する確率分布としている。 ただし、新規の文書dについて、 p(z \sim d)を生成することができないため、完全な生成モデルとは言えない。つまり学習に用いた文書以外にはうまく適用できる保証がない。

f:id:kibya:20200516153355p:plain
PLSA

一方、LDAは完全な生成モデルであり、2003年、David Bleiらの論文で発表された。 その後、LDAについては教師ありモデル、単語の順番を考慮したモデルなどへの拡張や、最適なトピック数の推定、変分ベイズによる効率的な推論方法など様々な研究の対象となっている。

確率モデル

まず、全体のイメージを下図で確認して欲しい。

f:id:kibya:20200506141706p:plain
LDAのイメージ

次に、具体的な確率モデルを説明する。 まず、トピックの総数をKとする(この数値は自分で設定する必要がある)。 文書はK種類のトピックがトピック比率 (topic proportion) \theta_d)と呼ばれる重みをつけて混合されたものと考える。
例えば、国会で新しい教育制度の法案が可決されたことに関するニュース記事だとすると、トピック「政治」の割合0.7、トピック「教育」の割合0.2、その他のトピックの割合0.1のようになるだろう。

 \theta _ d)はK次元のベクトルであり、各成分の和は1になる。
このような値を生成する確率分布としては、ディリクレ分布(ベータ分布の多次元版)が用いられる。

 p(\theta _ d) = {\rm Dirichlet} (\theta _ d \mid \alpha)

文書中の各単語はK種類のトピックのいずれかから生成されると考える。 まず、文書dにおけるi番目の単語 w _ {d,i}がどのトピックから生成されたかを示すトピック割り当て (topic assignment) z _ {d, i}を導入する。  z _ {d, i}はK次元のベクトルであり、いずれか1つの要素が1、他の要素は0である。

このような値を生成する確率分布としては、カテゴリ分布(ベルヌーイ分布の多次元版)が用いられる。

 p(z _ {d, i} \mid \theta _ d) = {\rm Categorical} (z _ {d, i} \mid \theta _ d)

次に、単語 w _ {d,i} z_ {d, i}で指定されたトピックから生成される。 ただし、単語 w _ {d,i}はいわゆるOne-of-K表現したベクトルである。すなわち、語彙の総数をVとすると  w _ {d,i}は長さVのベクトルであり、該当する単語のインデックスのみ1、他は0となっている。

LDAでは、単語の出現確率がトピックにより異なると仮定しており、トピックkにおける単語の出現比率はV次元ベクトル \beta _ kで表される。  \beta _ kの各成分は0から1の範囲の値であり、成分の和は1になる。
そこで、 w _ {d,i}もカテゴリ分布により生成されると考えることができる。この際、 z_ {d, i}のk番目の成分が1のとき、単語  w _ {d,i}はトピックkから生成されていることに注意する。

 p(w _ {d,i} \mid z _ {d, i}, \boldsymbol{\beta}) = \prod _ {k=1} ^{K} {\rm Categorical} (z _ {d, i} \mid \beta _ k) ^ {z _ {d, i, k}}

ただし、 \boldsymbol{\beta} = \{ \beta_1, \beta_2, \dots, \beta_k \} としている。

また、 \beta _ kの事前分布としては、ディリクレ分布を用いる。
 p(\beta _ k) = {\rm Dirichlet} (\beta _ k \mid \eta)

グラフィカルモデルで示すと、下記の通り。

f:id:kibya:20200506092410p:plain
LDAのグラフィカルモデル

実践

scikit-learnのLatentDirichletAllocationやgensimのLdaModelを使うと簡単に実装することができる。
長くなってしまうので、この記事はここまで。 Tensorflow Probabilityを使ったベイズ推論について、別の記事で説明する。

参考文献

  1. “Latent Dirichlet Allocation” D. Blei, A. Ng, M. Jordan, 2003 (http://www.jmlr.org/papers/v3/blei03a)
    LDAを提唱した元論文。

  2. Chapter5でLDAの変分ベイズによる推論方法を解説している。

  3. "Indexing by latent semantic analysis" S. Deerwester, S. Dumais, T. Landauer, G. Furnas, and R. Harshman, 1990 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.108.8490 LSI (LSA)の元論文。

  4. "Probabilistic latent semantic indexing" T. Hoffman, 1999(https://dl.acm.org/doi/10.1145/312624.312649) PLSIの元論文。

  5. PLSAの仕組みとその有用性についての記事 www2.deloitte.com

Optunaの概要

仕事でOptunaを使う機会があったので、論文を軽く読んでまとめてみる。

arxiv.org

Optunaとは

OptunaはPFNが開発したハイパーパラメータのチューニングツールであり、近年、勾配ブースティングやニューラルネットワークのチューニングによく使われている。
既存のツールに対するOptunaの特徴として、下記3点が挙げられている。

  1. Define-by-runで動的にパラメータ探索空間を設定可能
  2. 効率的なsampling, pruningアルゴリズム
  3. 設定が簡単で、軽量な実験から分散処理による大規模な計算まで対応可能

1. Define-by-run

define-by-runは深層学習フレームワークの分類に使われる用語だが、筆者らはパラメータチューニングツールにも同様の概念を提唱している。
従来のツールはdefine-and-run、すなわちパラメータの探索範囲を最初に設定した上で、最適なパラメータの探索を行うものだった。

一方、define-by-runでは動的に探索範囲を設定することができる。これは例えばニューラルネットワークのレイヤー数とユニット数を同時にチューニングする際などに役に立つ。

2. 効率的なsampling, pruningアルゴリズム

パラメータ最適化の効率は探索戦略(探索すべきパラメータを決める)と性能評価戦略(学習曲線から現在探索中のパラメータの良さを判断し、捨てるかどうか判断する)により決まる。

論文には探索(サンプリング)アルゴリズムの詳細は記載されていないが、TPE (tree-based Parzen estimator)*1, CMA-ES(Covariance Matrix Adaptation - Evolution Strategy)*2, GP-BO(Bayesian optimization) *3などを用いているとのこと。 サンプリング手法にはパラメータの相関を利用するrelational samplingと独立にサンプルするindependent samplingがある。TPEはindependent sampling、CMA-ES, GP-BOはrelational samplingである。

また、性能評価については、見込みの悪い試行を途中で終わらせること(pruning、枝刈り)が重要である。 OptunaではAsynchronous Successive Halving (ASHA)*4と呼ばれるアルゴリズムを採用している。ASHAでは試行中のパラメータセットの暫定的な順位に基づき、分散処理の個々のworkerが非同期的にpruningを行う手法とのこと。

使い方

Optunaではパラメータの最適化をobjective functionの最小化・最大化の問題と捉える。
最適化のプロセスをstudy、objective functionの1回の評価をtrialと呼ぶ。 パラメータの探索空間はtrialオブジェクトのメソッドにより動的に構築される。

以下のコードでは(モデルのチューニングではないが)、簡単な例として2次関数の最小化を行っている。

import optuna

def objective(trial):
    x = trial.suggest_uniform('x', -10, 10)
    return (x - 2) ** 2

study = optuna.create_study()
study.optimize(objective, n_trials=100)

study.best_params  # E.g. {'x': 2.002108042}}



以下はより実践的な、TensorFlowでニューラルネットワークのパラメータチューニングを行う場合のコード抜粋。 完全なコードは下記リンク先参照。 https://colab.research.google.com/drive/1MQp5rn6Dxhdv4r9PuJOyGf1sdsS2PHF1

ここでは、隠れ層の数、ユニット数、weight decayといったモデルのパラメータに加え、optimizerの種類とそのパラメータも最適化している。 従来のツールでこうしたチューニングを行うためには、最初にかなり複雑なパラメータ探索範囲を記述する必要があるが、Optunaのdefine-by-runの思想により、簡潔に実装することができる。

def create_model(trial):
    # We optimize the numbers of layers, their units and weight decay parameter.
    n_layers = trial.suggest_int("n_layers", 1, 3)
    weight_decay = trial.suggest_loguniform("weight_decay", 1e-10, 1e-3)
    model = tf.keras.Sequential()
    model.add(tf.keras.layers.Flatten())
    for i in range(n_layers):
        num_hidden = int(trial.suggest_loguniform("n_units_l{}".format(i), 4, 128))
        model.add(
            tf.keras.layers.Dense(
                num_hidden,
                activation="relu",
                kernel_regularizer=tf.keras.regularizers.l2(weight_decay),
            )
        )
    model.add(
        tf.keras.layers.Dense(CLASSES, kernel_regularizer=tf.keras.regularizers.l2(weight_decay))
    )
    return model


def create_optimizer(trial):
    # We optimize the choice of optimizers as well as their parameters.
    kwargs = {}
    optimizer_options = ["RMSprop", "Adam", "SGD"]
    optimizer_selected = trial.suggest_categorical("optimizer", optimizer_options)
    if optimizer_selected == "RMSprop":
        kwargs["learning_rate"] = trial.suggest_loguniform("rmsprop_learning_rate", 1e-5, 1e-1)
        kwargs["decay"] = trial.suggest_uniform("rmsprop_decay", 0.85, 0.99)
        kwargs["momentum"] = trial.suggest_loguniform("rmsprop_momentum", 1e-5, 1e-1)
    elif optimizer_selected == "Adam":
        kwargs["learning_rate"] = trial.suggest_loguniform("adam_learning_rate", 1e-5, 1e-1)
    elif optimizer_selected == "SGD":
        kwargs["learning_rate"] = trial.suggest_loguniform("sgd_opt_learning_rate", 1e-5, 1e-1)
        kwargs["momentum"] = trial.suggest_loguniform("sgd_opt_momentum", 1e-5, 1e-1)

    optimizer = getattr(tf.optimizers, optimizer_selected)(**kwargs)
    return optimizer


def objective(trial):
    # Get MNIST data.
    train_ds, valid_ds = get_mnist()

    # Build model and optimizer.
    model = create_model(trial)
    optimizer = create_optimizer(trial)

    # Training and validating cycle.
    with tf.device("/cpu:0"):
        for _ in range(EPOCHS):
            learn(model, optimizer, train_ds, "train")

        accuracy = learn(model, optimizer, valid_ds, "eval")

    # Return last validation accuracy.
    return accuracy.result()

とりあえず、特徴量を手当たり次第に作ってみてLightGBMの投入し、Optunaでハイパーパラメータ調整する、なんて作業は 誰か他の人にやってもらいたい。。

説明変数のある時系列予測(状態空間モデル)

概要

時系列データは実務で取り扱う機会が多い。その際、予測対象の時系列データ(目的変数)と似たような変動をする他の時系列データ(説明変数)を手がかりに予測したい場合がある。

説明変数を用いて目的変数を予測する最も単純な方法として、線形回帰がある。ただ、線形回帰では時系列データの特徴を捉えきれないという問題があり、高い精度は期待できない。

この記事では時系列分析の1手法である状態空間モデルによる予測を行い、線形回帰に比べて何が優れているかを説明する。

データ

ここでは、下図のように変動するX, Yという2つの変数を考える。 Yが予測対象の変数、Xは説明変数である。 Xを利用して、Yをリアルタイムに予測したい。つまり、t時点のYを予測する際には、t-1時点までのYの値とt時点までのXの値が利用可能とする。

f:id:kibya:20200104224855p:plain

XとYの散布図は下図のようになる。

f:id:kibya:20200104225925p:plain

XとYにははっきりと相関があり、XはYの予測に役立ちそうである。

線形回帰

ここで単純な予測モデル、すなわち線形回帰することを考える。
これはt-1時点までのX, Yを用いて回帰直線を作成し(切片αと傾きβを求め)、t時点のXの値を代入することでYの予測値を算出する。

 {
\begin{eqnarray} 
y_t &=& \alpha+ \beta x_t + \varepsilon_t \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\end{eqnarray} 
}

下図は、t=10までのデータ(グレー)を用いて、回帰直線を作成し、t=11のYを予測している。

f:id:kibya:20200104230756p:plain

全時点の予測結果をプロットすると下図の通り。

f:id:kibya:20200104231634p:plain

そこまで悪い予測ではないが、途中に連続して予測を大きく外してしまうケースが見られる。一期先予測誤差(MAE)は2.38程度となった。 これは線形回帰では、過去のデータ点をすべて同等に扱っており、直近のデータは類似しているというデータからはっきり読み取れる傾向を無視しているからである。
この直近のデータは類似している、という性質こそ時系列データの特徴であり、予測の際に考慮すべき情報なのだ。

状態空間モデル

時系列予測の1手法として、状態空間モデルがある。
特に今回扱うのは、一般的に「回帰成分を含むローカルレベルモデル」、と呼ばれるモデルで、式で表すと下記の通り。
線形回帰と違うのは、回帰直線の傾きと切片が、1つ前の時点の傾き、切片にノイズが加わって得られる、としている点である。

 {
\begin{eqnarray} 
y_t &=& \alpha_{t} + \beta_{t} x_t + \varepsilon_t \\
\alpha_{t} &=& \alpha_{t-1} + \eta_{\alpha,t} \\
\beta_{t} &=& \beta_{t-1} + \eta_{\beta,t} \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\eta_{\alpha,t} &\sim& N(0, \sigma^2_{\eta_{\alpha}}) \\
\eta_{\beta,t} &\sim& N(0, \sigma^2_{\eta_{\beta}}) \\
\end{eqnarray} 
}

回帰直線の切片と傾きが「状態」と呼ばれる潜在変数である。各時点の状態は1時点前の状態によって決まる。
つまり予測の際は直近のデータの影響が大きくなる。

状態の推定にはカルマンフィルタと呼ばれる手法が用いられる。詳しくはこちらの記事を参照。 logics-of-blue.com

カルマンフィルタの要点としては、観測値が得られるたびに状態の推定値が補正され、補正された状態を用いて次の時点における状態と観測値の予測が行われる、ということになる。

状態空間モデルによる一期先予測の結果は下図の通り。

f:id:kibya:20200107075456p:plain

一期先予測誤差(MAE)は1.81程度となり、線形回帰よりはっきり改善している。

予測結果の比較

2つのモデルの違いを見るため、t=47, 48における予測結果を比べてみる。

まず、線形回帰の場合。この場合は過去の全ての観測値を同等に扱って予測を行う。結果、t=47, 48のいずれでも実際の観測値より大幅に小さい値を予測してしまっている。

f:id:kibya:20200107080906p:plainf:id:kibya:20200107080925p:plain

次に、状態空間モデルの場合。t=47では予測誤差が大きいのは同じだが、t=48の予測誤差は比較的小さくなっている。 これは、t=47の観測値を用いた状態の補正が行われているためである。つまり、直近の観測値は予測値よりだいぶ大きかった、という情報をモデルに与えることで、次の時点の予測値を大きくするように補正が働く

f:id:kibya:20200107081003p:plainf:id:kibya:20200107081013p:plain

アニメーションも載せておく。

f:id:kibya:20200111084459g:plainf:id:kibya:20200111084530g:plain
線形回帰(左)、状態空間モデル(右)

まとめ

時系列データの分析には状態空間モデルなど時系列モデリングの適用をしましょう。