機械学習の小部屋

機械学習入門～病床機能報告公表データの分析～

Medical Information Analysis
－医療情報分析入門－
「はじめの一歩」

Aコース
機械学習入門－病床機能報告公表データの分析－
テキスト

この動画は、2019.12.22に開催されたセミナーの模様を撮影したものです。テキストを見ながら視聴してください。ただし時間の関係でテキストに記載されている内容の多くを省略しています。また、説明の都合や講師のミスで説明の順序がテキストと違ったり、間違った説明をしているところもあります。そのあたりについてはテキストしたがって正しい操作を行ってください。

[9]Mining fda drug labels using an unsupervised learning technique

論文概要

医薬品添付文書（FDA-approved drug labels）にLDAを適用して医薬品の安全性や効能の観点からグルーピングを行ったという論文。
医薬品添付文書が解析対象のドキュメント。ただし、そのまま使うのではなくMedDRAを使って標準的なADR用語に変換している。そして安全性や効能のパターンをLDAを用いてトピックとして抽出するというシナリオ。

モデル選択

トピック数は天下り的に100としている。根拠は示されていない。

評価方法

トピックを構成する医薬品に安全性と薬効の観点から共通性があるかどうかで評価を行っている。安全性については添付文書に囲み警告文書（BW: Boxed Warning）があるものがどの程度あるか、効能については解剖治療化学分類法（ATC: Anatomical Therapeutic Chemical Classification System）で分類される医薬品がどの程度あるか調べて、それらが7割以上あれば有効であるとしている。

読後感

評価方法が理解しづらかった。
医薬品にトピックを対応させる際、その医薬品のトピック含有率が最大のトピックを割り当てている。つまりこの段階で他のトピックの役割を切り捨てている。
次に、こうして各医薬品に割り当てられたトピックのうち、少なくとも10以上の医薬品を持つトピックのみ残して他は解析対象から除外した。
これは意味のある統計解析とバイアスを取り除くためと述べている。

Bisgin H, Liu Z, Fang H, Xu X, Tong W. Mining FDA drug labels using an unsupervised learning technique--topic modeling. BMC Bioinformatics. 2011;12 Suppl 10(Suppl 10):S11. Published 2011 Oct 18. doi:10.1186/1471-2105-12-S10-S11

医療分野でトピックモデルを利用した研究（その４）

[8]Evaluating disease similarity using latent Dirichlet allocation

テキストで記述された疾患からLDAを用いて疾患間の類似性を求める研究。

使用するテキストはOMIM（Online Mendelian Inheritance in Man：ヒトの遺伝子と遺伝子によって規定される表現型のカタログ）。

LDAによって（疾患ラベルを割り当てた）1,311のドキュメントのトピック分布を求め、トピック分布間のKL-divergenceで疾患間の類似性を測っている。

評価方法：この手法で求めた疾患間の類似性をDO（Disease Ontology）による基準（DOの階層で3階層以内に同一の祖先がいれば同じ疾患とみなすという基準）をGold Standardとして感度特異度分析を行った。

疾患1	疾患2	KL-divergence
D₁	D₂	d_1,2
D₁	D₃	d_1,3
・・・	・・・	・・・
D₁	D₁₃₁₀	d_1,1310

図1．疾患間の類似性

ある疾患Diとそれ以外の疾患Dj (i != j)の組み合わせについてトピック分布のKL-divergenceを計算し（図1）、KL-divergenceの昇順（類似している準）に並べ替え、上からn個がPositive、残りがNegativeとしてGold Standardを使ってnを変えながら感度特異度分析を行う。これによってAUCを1つ得ることができる。

図2．ROC曲線（原論文のFigure 2）

これをすべてのiについて行えば、全部で1,311個のAUCが得られる。

モデル選択（トピック数の決定）：トピック数を変えながら感度特異度分析を行って最大AUCを与えるトピック数を最適なトピック数とした。

図3．平均AUCのトピック数依存性（原論文のFigure 3）

なお、図3において抽出単語数を5k, 7k, 10kと変えて実験を行っている。ここで、抽出単語とはあらかじめ設定した基準を上回る出現頻度の単語の集合のことである。

読後感

テキストを対象とした典型的なLDAの応用例。疾患に関するテキストをコーパスとし、LDAによってテキストのトピック分布を求めている。これによって疾患のベクトル表現を得ている。これは広大な単語空間を小数次元のトピック空間に次元削減するのと等価で、これもLDAの一側面を表している。
興味深いのは、テキストを単語空間のベクトルとみなして主成分分析（PCA）を行い、次元削減したものとLDAを比較していることである。
トピック数と同数の主成分を残し、その主成分空間上のユークリッド距離で疾患間の類似性を定義し、LDAの場合と同じやり方でDOをGold Standardに用いて感度特異度分析を行ったところ、なんとほぼ同様な特性を示した。

図4．LDAとPCAにおけるAUCのトピック数（主成分数）依存性

図4では、トピック数がほぼ20あたりでLDAと主成分数が同数のPCAが同程度のAUCを実現している。この結果からトピック数は主成分分析における主成分の数とほぼ同じ役割をしていることがわかる。すなわちデータを要約する上で最適な次元数を表しているものと考えられる。
この知見は当然と言えば当然かもしれない。なぜならLDAであれPCAであれ、次元削減モデルは同じ最適な次元を導き出すだろうから。もし、モデルによって異なる次元数が得られたなら、むしろその方が奇妙である。

James M Frick, Rajarshi Guha, Tyler Peryea, Noel T Southall: Evaluating disease similarity using latent Dirichlet allocation. preprint and has not been peer-reviewed, https://www.biorxiv.org/content/early/2015/11/03/030593, 2016, (access: 2019/1/15).

医療分野でトピックモデルを利用した研究（その３）

[7]Predicting inpatient clinical order patterns with probabilistic topic models vs conventional order sets

【論文概要】

トピックモデルを使って入院患者に発行するオーダパターンを予測するという論文。入院患者がドキュメント、オーダが単語、そしてオーダパターンがトピックになるのかな・・・。
評価はモデルが予測するオーダパタンと人手で作成した既存のオーダセットと比較して行っている。
評価指標として精度（Precision）と適合率（Recall）を次のように定義する。

精度（Precision）: モデルが予測したオーダのうち実際に発行されたオーダの割合
適合率（Recall）: 発行されたオーダのうちモデルが予測したオーダの割合

ここでモデルは20,000人の患者からランダムに抽出した10,655人（～50％）の患者のオーダを用いて構築した。
評価は20,000人の患者からランダムに抽出した4,820人（～25％）の患者の各々に対して入院後最初のオーダセットが発行されてから時間 t 以内に発行された実際のオーダ（これを「正しい」オーダセットと呼ぶ）とモデルが予測したオーダについて表1に示す混同行列を作成して精度と適合率を計算した。

表1．混同行列

		そのオーダは発行されたか
		Yes	No
そのオーダは予測されていたか	Yes	p	q
そのオーダは予測されていたか	No	r	s

精度 = p ÷ (p + q)
適合率 = p ÷ (p + r)

なお、あらかじめ人手によって作成された既存のオーダセットに含まれるオーダはすべて発行されるとは限らないので、実際に発行された「正しい」オーダセットと既存のオーダセットは一般に異なる。

ここで、p + r は時間 t 以内に発行された「正しい」オーダセットの数で、そのうちモデルによって予測されたオーダの数が p である。また、 p + q はモデルが予測したオーダの数で、そのうち時間 t 以内に発行されたオーダの数が p である。s はトピックモデルに基づくオーダセットにも「正しい」オーダセットにも入っていないオーダ数であるが、これがなくても精度と適合率の計算には支障はない。

ある患者（これを k とする）に対してトピックモデルに基づくオーダセットを以下のようにして作成する。
まず患者 k のトピック含有率P(Topic_j|Patient_k)を推定し、トピック j にオーダ i が含まれる確率P(Item_i|Topic_j)を用いて、患者 k にオーダ i が発行される確率P(Item_i|Patient_k)を次式で計算する。

P(Item_i|Patient_k) = Σ_j P(Item_i|Topic_j)×P(Topic_j|Patient_k)

確率P(Item_i|Patient_k)の降順にN個のオーダを並べたものがトピックモデルに基づくオーダセットとなる。ここで、N は既存のオーダセットに含まれるオーダ数である。

トピック含有率P(Topic_j|Patient_k)は患者 k に発行された実際のオーダデータを使ってトピックモデルが推定する。その際、患者 k は一般にモデルを作成するのに使用した訓練データではなく検証用データである。従ってモデルから検証用データに対するトピック含有率を求める必要があるが、Rのtopicmodelsパッケージを調べてもそのようなメソッドは見当たらない。しかし、PythonのGensimのmodels.ldamodelライブラリにはget_document_topicsという与えられたドキュメントからトピック分布を得るメソッドがある。

次に、表1の混同行列においてトピックモデルによる予測の当否の代わりに既存のオーダセットに含まれているかどうかで区分するとそのオーダセットの精度と適合率が計算できる。これとトピックモデルによる精度と適合率を比較することによりトピックモデルの評価を行う、というのがこの論文での評価方法である。

図1（論文ではFigure 4）に比較した結果を示す。

図1．オーダセットの精度の比較

図1の横軸の t は最初にオーダセットが発行されてからの時間（検証時間：followup verification time）である。t によって精度や適合率がどのように変化するのかを見ているのだろう。

最上段のグラフは「正しい」オーダの数の時間変化である。時間とともに「正しい」オーダ（Possible Items）の数は増加する。なぜなら実際に発行されたオーダを「正しい」オーダとしているからである。また、実際に発行されるオーダの中にはモデルが予測したオーダセットに含まれるオーダ（図中のTopic Models）や既存のオーダセットに含まれるオーダ（図中のOrder Sets）があり、その数も時間とともに増えている。

二段目のグラフはトピックモデルと既存のオーダセットの適合率の時間変化である。実際に発行されるオーダの中に、最初のころはトピックモデルが予測するオーダや既存のオーダセットに含まれるオーダの割合は高いだろうから適合率は高いが、時刻とともにその割合も減ってくるだろうから適合率も低くなっていくものと考えられる。

三段目のグラフはトピックモデルと既存のオーダセットの精度の時間変化である。
t が小さいうちは実際に発行されるオーダの数は N に比べて少ないので精度（p / N）は低いが、時間とともに発行されるオーダの数が増えると、その中に含まれるトピックモデルが予測するオーダや既存のオーダセットに含まれるオーダの割合も増えていくので精度も上がっていく。

最下段のグラフはAUC（Area Under Curve）をプロットしたものである。AUCを算出するにはROC曲線をプロットする必要があるが、そのためには多数の（精度, 適合率）の組が必要である。論文のFigure 3の説明には次の記述がある。

Average discrimination accuracy (ROC AUC) when predicting additional clinical orders occurring within t followup verification time of the invocation of a pre-authored order set during the first 24 hours of hospitalization for 4820 validation patients.

また、figure 4 (D)の説明には次の記述がある。

evaluating full ranking of possible orders scored by topic model or included/excluded by order sets

さらに、本文中には次のように書いてある。

Prediction of these subsequent orders is evaluated by the area under the receiver operating characteristic curve (c-statistic) when considering the full score-ranked list of all possible clinical orders. Existing order sets will have N suggested orders to choose from, so we evaluated those N items vs the top N score-ranked suggestions from the topic models toward predicting subsequent orders by precision (positive predictive value) at N and recall (sensitivity) at N.

N を大きくしていくと実際に発行されたオーダはほとんどモデルが予測するオーダの中にあるだろう。すると、表1において p は r に比べて大きくなるので適合率は1に近づく。一方、N を大きくしていくとモデルが予測するオーダにあって実際に発行されないオーダの数も増えるだろうから q は p に比べて大きくなり、精度は0へ近づく。

一方、N を小さくしていくと実際に発行されたオーダのうちモデルが予測したオーダの数は減るので r に比べて p は小さくなり、適合率は小さくなる。その反面、N 自体が小さくなるのでモデルが予測するオーダのトップ N は実際に発行される確率が高くなり、p は q に比べて相対的に大きくなり、その結果精度が大きくなる。

そのため図2のようになると考えられる。

図2．ROC曲線

【モデル選択】

モデル選択、すなわちトピック数Kの決定はAUCのトピック数依存性（図3）をプロットして決定している。

図3．AUCのトピック数依存性

これは、最も「正しい」オーダセットを予測するトピック数を最適なトピック数とみなす考え方である。確かに現実を最も忠実に再現するモデルが良いモデルに違いないのでこの考え方は合理的である。

【読後感】

非常にシンプルなLDAの医療への適用例である。特徴はトピック数Kの決定法とモデルの評価法にある。モデルが実際に発行されたオーダをどのくらい忠実に再現（予測）しているかによってモデルのパラメタであるトピック数を決定している。これはPerplexityによるトピック数決定の手法やheld-outデータによる尤度最大基準による方法と等価である。
モデルの評価に人手で作成した既存のオーダセットを外部基準として用いている。

果たして既存のオーダセットはモデル評価のための外部基準となりえるのだろうか
オーダセットには様々な種類があると思われるが、それを無視して単にP(Item_i|Patient_k)だけでオーダを予測していいのだろうか
「正しい」オーダセットや既存のオーダセットにはオーダの「発行順序」という概念があるが、トピックモデルによるオーダセットにはない。あるのはその患者に対してそのオーダがどの程度発行され易いかという確率だけである。しかるに精度や適合率を計算するのにオーダの出現準を暗黙に仮定しているが、これはよいのだろうか

といった疑問が残った。

Jonathan H Chen, Mary K Goldstein, Steven M Asch, Lester Mackey, Russ B Altman: Predicting inpatient clinical order patterns with probabilistic topic models vs conventional order sets. Journal of the American Medical Informatics Association, 24(3), 472-480, 2017.

モデル選択

LDAにおけるモデル選択とは最適なトピック数の決定である。最適なトピック数の決定方法として次のようなものがある。

Perplexity
log-likelihood
predictive log-likelihood
HDP (Hierarchical Dirichlet Process)

Perplexityはモデルがどのくらい単語を絞ることができるかという指標。小さいほど良い。これを最小とするトピック数を最適なトピック数と考える。

図1．Perplexityのトピック数依存性

図1はDPCデータに対してPerplexityのトピック数依存性を描いたものである。トピック数が500になってもまだPerplexityは減少傾向にある。

log-liklihoodは訓練データを用いて作成したモデルのパラメタを使って訓練データの対数尤度（の調和平均）を求めたものである。あるブログ記事を読むと、GriffithとSteyvers (2004) によってこの近似法が提案されたとのことである。
DPCデータに対してこの値のトピック数依存性を描いたものを図2に示す。

図2．log-likelihoodのトピック数依存性

図2を見ると、トピック数が350あたりでlog-liklihoodが最大になっている。

predictive log-likelihoodは前回のブログにも書いたように、未知のデータに対して求めた対数尤度である。訓練データで求めた対数尤度はデータにオーバーフィッティングする可能性があるのでこれを避けるのがpredictive log-likelihoodである。
しかし、LDAに対してこれを計算するライブラリが見当たらず計算しあぐねている。

HDPはデータから自動的に最適なトピック数を求める手法である。まだ勉強していないので詳細はわからない。

これら以外に、先に紹介したあるブログ記事にIdtuningというRパッケージで最適なトピック数を求めることができると紹介されていた。
ここでは、次の4つの指標が紹介され、それらを計算してトピック数依存性をグラフ表示する例が掲載されている。

上述したブログ記事によれば最初の3つ（1～3）はトピック間の単語分布の類似性に注目した手法であるとのこと（2はcosine類似性、3はKLダイバージェンスによって類似性を計算し、1も何らかの距離を計算しているらしい）。そして4はlog-likelihoodを用いる方法と同じである。
DPCデータに対してこのIdtuningを使って4つの指標のプロット依存性をプロットしたものを図3に示す。

図3．Idtuningの結果

図3によるとArun2010とCaoJuan2009はトピック数が250～300あたりで最小（最適）になっており、Deveaud2014はトピック数が100で最大（最適）になっている。そしてGriffiths2004はトピック数が450で最大（最適）になっている。

ドキュメントに相当する病院数が1600程度しかないのにトピック数が300～400というのはあまりに多すぎはしないだろうか。

医療分野でトピックモデルを利用した研究（その２）

[4]Identifying Patterns of Associated-Conditions through Topic Models of Electronic Medical Records

【論文概要】

電子カルテ（EMR）に蓄積された患者記録をドキュメント、SNOMED-CTでコーディングされた患者の健康状態（health conditions）を単語とみなしてトピックモデルを適用し、潜在的な健康状態のパターン（latent patterns of associated-conditions：これがトピック）を抽出するという研究。トピック数は20としている。

構築したトピックモデルを評価する方法として以下の2つの評価を行った。

医学的な妥当性の評価（定性的な評価）
弁別性（distinctiveness）と気密性（tightness）による定量的な評価

前者は、同一トピック内に高確率で出現する健康状態が関連しあっているかどうかを医学文献で検証することによって妥当性の評価を行っている。

後者の弁別性（distinctiveness）は、トピックの単語分布間のJensen-Shannon Divergenceを計算してトピック間の違いを定量的に評価している。
気密性（tightness）は1つのトピックがどのくらい少ない健康状態（SNOMED-CTコード）で表されるかを計算して評価している。

代表的な6つのトピックについて、出現確率の大きいものを10個抽出したところ、疾患の共起が医学文献で裏付けられたトピックが得られた。

20個のトピック間のJensen-Shannon Divergenceを計算したところ、平均値は0.666、中央値は0.692、そして最小値は0.483となった。Jensen-Shannon Divergenceは、2つの確率分布間の"距離"を表す量で、2つの分布が同一の場合は0に、まったくオーバーラップしない場合はln2 (～0.693)になるので、この結果は十分な弁別性を示している。
トピック内の出現確率が閾値（0.01）を超える健康状態（SNOMED-CTコード）の数を数えたところ、180あるSNOMED-CTコードのうち10以下のコードを合わせたものが9割を超えており、このことから十分な気密性を示していることがわかる。

【読後感】

LDAのテキストマイニング以外への応用としては興味深いものであった。しかし、医学文献による定性的な評価方法は説得力に欠く印象をぬぐえない。すべてのトピックについて検証したわけでもなさそうだし、類似したSNOMED-CTコードがLDAによってうまくグルーピングされたという客観的な証拠にはならないような気がする。
定量的評価方法については、弁別性（distinctiveness）といい気密性（tightness）といい計算が簡単な割には興味深い指標だと思った。しかし、これらはトピック間は健康状態空間における分布として大きく異なり、各トピックはわずかな健康状態で記述できる（高次元健康状態空間上のトピックが低次元健康状態空間で表現されている）ということを言っているだけで、トピックの解釈については何も触れていない。
最後に、なぜトピック数が20なのか根拠が示されていない（トピック数を14にしても結果はあまり変わらなかったとは書いていたが）。

Moumita Bhattacharya, Claudine Jurkovitz, Hagit Shatkay: Identifying Patterns of Associated-Conditions through Topic Models of Electronic Medical Records. 2016 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), 466-469, 2016.

[5]Interpretable Topic Features for Post-ICU Mortality Prediction.

【概要】

ICU退院後の死亡率（post-discharge ICU mortality）の正確な予測と透明性のある予測モデルの作成を目的に、ICD-9-CMコードをラベルに用いたLabeled-LDA (Latent Dirichret Allocation)を使って患者の診療記録（medical notes）から理解可能なトピック特性表現（topic feature representations）を取り出した。

近年、トピックモデルを用いて死亡率の予測精度を改善しようという試みが行われているが、トピックそのものはフラットな単語の集まりに過ぎず、その臨床的な解釈には専門家による吟味が必要となる。

この研究では解釈可能なトピックを自動的に定義する手法を提案している。その手法とは、ICD-9-CMコードをLabeled-LDAのトピックに用いてモデルの学習を導き、診療記録から理解可能なトピック特徴表現（understandable topic feature representations）を抽出するという方法である。

Labeled-LDAのラベルとしてICD-9-CMを用いる利点には次の2つがある。

ある患者の診療記録はその患者に割り当てたICD-9-CMコードに対応するトピックにのみ寄与する
トピックの解釈はICD-9-CMコードの定義とそのトピックに含まれる単語集合をとおして達成される

モデルを訓練する段階ではICD-9-CMコードを使っている。しかしながら、患者が退院してしまうとICD-9-CMコードは利用できなくなるので、予測モデルの特徴量としてICD-9-CMコードを直接含めることができない。Labeled-LDAを使えばこの問題を回避できる。なぜならば、予測段階ではICD-9-CMコードに対応するトピックの割合を推測するのにICD-9-CMコードを必要としないからである。これは、トピックに高確率で含まれる単語を特徴量に用いて文書のトピック含有率を予測するからである。

【評価方法】

提案手法の評価は、得られたトピック含有率を特徴量に用いてSVMでICU退院後の死亡率（post-discharge ICU mortality）を予測することによって行っている。すなわち、外部基準として予測の精度を用いているわけである。死亡率の予測は、もとをただせば死亡／生存の分類問題であるから本質的にこれまで読んできた文献と同じである。トピックモデルの評価はやはりこれしかないのだろうか。
予測する死亡率は30日後と6か月後の2つで、6か月後死亡率の予測SVMモデルの特徴量に用いるのは以下の3パターンである。

ベースライン：年齢、性別、入院時のSAPS-Ⅱスコア（ICUの重症度スコア）、最小SAPS-Ⅱスコア、最大SAPS-Ⅱスコア、Elixhauser Comorbidity Index（30個のへ依存症を考慮に入れたICU重症度の補正指標）
ベースライン＋通常のLDAから得られた50個のトピック含有率
ベースライン＋Labeled-LDAから得られた50個のトピック含有率（提案手法）

30日死亡率予測モデルでは、最小SAPS-Ⅱスコア、最大SAPS-Ⅱスコア、ICD-9-CM由来のElixhauser Comorbidity Indexを特徴量から取り除いている。これは、ICD-9-CMコードの割り当ては通常退院後2週間まで有効とされているからで、それ以降はICD-9-CM由来のElixhauser Comorbidity Indexを特徴量に用いることは適切ではないからである。

【結果】

最も予測精度が高かったのは、baseline + LDA with 50 topicsで、30日死亡率の予測でAUCが0.860、6か月死亡率は0.842だった。次いで baseline + Labeled-LDA with 111 labelsで、各々0.835、0.829だった。これに対してbaselineのみでは、それぞれ0.736、0.776であった。
この結果からLDAで得たトピック成分を特徴量に用いると予測精度が高くなることが分かったが、Lableled-LDAよりも純粋なLDAの方が精度が高いという結果となった（Table 2）。これは、本文に書かれていることと矛盾する。
しかし、Discussionに書かれているように純粋なLDAは専門家によるトピックの解釈が必要になる。これに対してLabeled-LDAを用いればTable 3に示すようなトピック別（つまりICD-9－CMコード別）死亡率を推定できる。
また、興味深いことにLabeled-LDAモデルは、トピックに関連付けられた単語によって、潜在的に異なる疾患の間の関係を発見する能力を持っている。例えば、"Other metabolic and immunity disorders（他の代謝疾患及び免疫疾患）"に関連付けられたトピックに頻出する単語に、chest（胸部）、artery（動脈）、coronary（冠状）、cabg（冠動脈バイパス手術）などがあるが、これは循環器疾患と代謝疾患に関連性があることを示している。これは、ラベルの頻度や相互依存性を考慮したLabeled-LDAを用いれば異なるラベル間の相関を探索できる可能性があることを示唆している。

【読後感】

Labeled-LDAというものをはっきりと勉強したわけではないのでぼんやりとしか理解できないのだが、この手法は、LDAの大きな欠点である「得られたトピックの解釈困難性」を緩和する手法のように思われる。
モデルを学習するときラベルを与えることにより、トピックを明示的に定義できる。つまり、当該トピックに関連するドキュメントのみを与えることによってトピックの意味をモデルに教えている。LDAと同じように得られたトピックは単語分布になる。これを用いてトピックの予測は単語を使って行うが、トピックにはラベルが対応しているので、それを用いてトピックの意味を解釈できる。つまり、こうだ。モデルに文書を与えると、モデルは文書のトピック分布を返す。各トピックにはラベルがついているので、その文書の解釈が可能となる。

これは文献[2]と手法が似ていると感じた。文献[2]では請求データというドキュメントから疾患というトピックを推測し、その中で最もコストのかかったトピック（疾患）を最大資源投入疾患としている。そこでもLabeled-LDAを用いている（ラベルはDPCコード？）。

この論文を読んで感心したのは、

退院後はICD-9-CMコードが付与されていないにも拘わらず、ICD-9-CM別の死亡率を予測できること
トピックに高確率で含まれる単語を用いて疾患間の関係を調べることができること

である。

Yen-Fu Luo, Anna Rumshisky: Interpretable Topic Features for Post-ICU Mortality Prediction. AMIA 2016 Annual Symposium Proceedings, 827-836, 2016.

[6]Identifying prescription patterns with a topic model of diseases and medications

【論文概要】

DMPM（disease-medicine pattern model）という拡張版トピックモデルを使って大量の保険請求データを解析し、疾患と処方薬の関係から処方パターンを抽出した。
この研究では処方せんをドキュメント、その中に含まれる疾患（ICDコード）や処方薬を単語、処方パターンをトピックとみなしている。
DMPMの生成プロセスを図1に示す。

図1．DMPMの生成プロセス

また、グラフィカルモデルを図2に示す。

図2. DMPMのグラフィカルモデル

DMPMはLDAを拡張して各々の処方せんについて処方薬と疾患を生成するプロセスを加えている。
疾患は、図1のBのiiで生成し、医薬品はBのiiiのaで生成している。
通常のLDAであれば、図2において、ハッチのかかったノードDが単語に相当し、一番外側の四角はドキュメントの繰り返し、その内側の四角は単語の繰り返しを表し、ノードMに相当する確率変数はない。
これに対してDMPMでは、個々の患者pに対してトピック分布（この論文ではトピックという用語ではなくパターンという用語を使っている）θ_pをαをパラメタとするディリクレ分布から生成し、その患者に対して発行された個々の処方せんのトピックz_p,rをθ_pをパラメタとする多項分布から生成する。疾患d（図2中のノードDのインスタンス）は、トピックz_p,r=kに対するβ_kをパラメタとする多項分布から生成する。また、処方薬m（図2中のノードMのインスタンス）は、トピックz_p,r=kに対するγ_kをパラメタとする多項分布からN_pr個生成する。
このモデルで推定するパラメタは

θ_p = {θ_p,k}：患者pにトピックkが含まれる割合（患者pのトピックk含有率）
β_k = {β_k,d}：トピックkに疾患dが含まれる割合（トピックkの疾患d含有率）
γ_k = {γ_k,m}：トピックkに医薬品mが含まれる割合（トピックkの医薬品m含有率）

の3つである。

【評価方法】

得られたモデルの評価方法としては、これまでは最大トピック含有率を与えるトピックでドキュメントを分類したり^[1]、トピック含有率を新たな特徴量に使って分類問題の精度を比較する^[5]といった外部基準を使った評価手法が多かった。
この研究では{β_k,d}や{γ_k,m}でパターン（トピック）間の類似度（cosine類似度、K-L divergence、Jensen-Shannon divergence）を計算し、他の手法（ICDコードによるグルーピングや疾患LDA、医薬品LDA）に比べて年齢、性別、居住地などの患者属性による違いをどのくらい顕著に区別できるかを定量的に評価している。
一例として年齢による患者の多様性をDMPMとICDによるカテゴリ化で比べたのが図3である。

図3. 年齢による多様性の捕捉（上：DMPM、下：ICD）

図4．他のグルーピング手法との比較

【読後感】

著者らはICDのような標準コードを用いた疾患のカテゴリ化を「トップダウン方式」と呼び、それに対して提案手法によるデータに基づいたカテゴリ化を「ボトムアップ方式」と呼んで両者の比較を行っている。
また、それらの類似点や差異をもとに本手法の妥当性（ICDによるカテゴリをある程度再現しているから妥当だ）の評価と有効性（ICDによるカテゴリ化では説明できない事実を捉えているので有効だ）を議論している。
たとえば同じ疾患でも使用する医薬品が異なるとか逆に異なる疾患に同じ医薬品を処方しているかといった知見を得ている。
この研究の興味深い点は疾患と医薬品の両方を単語とみなしてLDAを拡張しているところである。これによってパターン（トピック）を疾患含有率ベクトル{β_k,d}と医薬品含有率ベクトル{γ_k,m}の2通りで表現し疾患と医薬品の相関を議論している。

ドキュメントにラベルが付けられていたら図3のようなグラフを描いて、トピックモデルを使った場合とラベルだけを使った場合でどのような違いが見られるか（ラベルだけでは捕捉できないどのような特徴をトピックモデルで捉えることができるか）を議論するのは興味深い。

また、図5のようなグラフを描いてラベル別のトピック含有率を比較するのも面白いかもしれない。

図5．ICDによるカテゴリ化とDMPMによるパターンとの比較

【モデル選択】

この論文ではモデル選択（最適なトピック数の推定）手法としてPerplexityとpredictive log-likelihoodの2つの方法を用いてトピック数＝15が最適であるとしている。
RのtopicmodelsパッケージにはPerplexityを計算するメソッドはあるが、predictive log-likelihoodの計算メソッドは見当たらない。
predictive log-likelihoodは未知データR_newを使ってP(R_new)で計算する。topicmodelsのLDAメソッドは訓練データからモデルを作成する。LogLiksというメソッドを使えば作成したモデルからlog-likelihoodを計算できるが、これはP(R_train)なので、トピック数Kを増やせばいくらでも大きくすることができる（オーバーフィッティング）。
P(R_train)を計算するメソッドはないのだろうか。"Perplexity To Evaluate Topic Models"にはPerplexityはpredictive log-likelihoodと等価であると書いてあるが・・・。

P(R_new)を計算する方法を考えてみよう。LDAの生成モデルでは、ドキュメントr∈R_newのトピック分布θ_kからトピックkを選び、トピックkの単語分布φ_k,wから単語wを選択した結果、単語wがrに現れると考える。これをR_newだけで計算するにはどうすればよいのだろう？そもそも、単語wがどのトピックから生成されたなどどうしてわかるのだろうか？単語wはトピックkからφ_k,wの確率で生成されるのだから単語だけではトピックは決定できない。
これがpredictive log-likelihoodの計算メソッドが提供されていない理由なのかもしれない。
しかし、この論文では(8)式を使ってpredictive log-likelihoodを計算している。ポイントはwのトピックkなどわからないので、期待値E[θ_k|R_obs,Θ]を使って近似しているところだ。これを使えば最終的にpredictive log-likelihoodはモデルのパラメタ{β_k,d}、{γ_k,m}を使って(9)式で求めることができるとのこと（途中の式の展開はよくわからないのでフォローできていない）。これを用いて図6に示すpredictive log-likelihoodのトピック数K依存性が得られたということである。

図6．predictive log-likelihoodのパターン数（トピック数）依存性

訓練データを用いて計算したlog-likelihoodはトピック数の増加に伴っていくらでも増加する。これはモデルがデータにオーバーフィッティングしているためだと考えられる。したがってトピック数Kの決定にはpredictive log-likelihoodを用いるべきではないだろうか。
論文中の式(8)のE[θ_k|R_obs,Θ]はどのようにして求めるのだろうか？

Park S, Choi D, Kim M, Cha W, Kim C, Moon IC: Identifying prescription patterns with a topic model of diseases and medications. Journal of Biomedical Informatics archive, 75C, 35-47, 2017.

LDAvisの日本語文字化け

LDAvisの文字化けの原因

LDAvisは以下の5つのファイルを作成する。

index.html
lda.css
d3.v3.js
ldavis.js
lda.json

この中で、lda.jsonにLDAモデルから作成されたデータが入っている。
LDAvisパッケージの関数createJSONのvocabに日本語を設定すると文字化けが起きて画面に何も表示されない。原因は文字コード体系にある。

json <- createJSON(
 phi = posterior(model)$terms
 theta = posterior(model)$topics
 vocab = colnames(posterior(model)$terms)
 doc.length = rowSums(as.matrix(dpc_dtm))
 term.frequency = colSums(as.matrix(dpc_dtm))
 mds.method = svd_tsne
)

リスト1．createJSON

LDAvisはcharset="utf-8"でindex.htmlを作成する。したがってlda.jsonもutf-8で作成されなければならない。しかし、Windows版のR Studioの内部コード (option()$encoding) はShift_JIS (CP932)なので、lda.json中の日本語部分はShift_JISで書き出され、index.htmlとの間に齟齬が生じて動かなくなるようだ。

対策

調べた限り、LDAvisにはエンコードを指定するオプションはないようなので、Windowsの内部コードを一時的にUTF-8に変更して処理を行い、終わったら元に戻す（参考サイト）。

tmp.enc <- options()$encoding
options(encoding = "UTF-8") 

dpc_df <- read.table("input.csv", header=T, sep=",", fileEncoding = "CP932")

・・・LDAvisのメイン処理・・・

options(encoding = tmp.enc)

リスト2．LDAvisの処理

R Studioの内部エンコードはoptions()$encodingに格納されているので一時的にこの値をテンポラリ変数tmp.encに退避したのちにoptions(encoding = "UTF-8")でUTF-8を設定する。
読み込むCSVファイルがShift－JISで作成されている場合はread.tableにオプションfileEncoding = "CP932"を付けて読み込む。ファイルのエンコードが"UTF-8"であればこれは不要である。
その後、LDAvisのメイン処理を行って、処理が終わったら退避していた元のエンコードに戻す。

トピックモデルの可視化 - LDAvis

概要

LDAの結果を可視化するツールLDAvisを使ってモデルを俯瞰する。

モデルの構築

モデルを構築するRスクリプトをリスト1に示す。

################################################
#
# モデルの生成
#
################################################
# log likelihood の最大値を与えるトピック数＝350
topic_num = 350
burnin = 500
iter = 1000
keep = 10
model <- LDA(dpc_dtm, k = topic_num, method = "Gibbs", control = list(burnin = burnin, iter = iter, keep = keep) )
model

LL <- data.frame(
  topic_nums = seq(10, 1500, 10),
  logLiks = model@logLiks
)

################################################
#
# 対数尤度の変化
#
################################################
ggplot(LL, aes(x=topic_nums, y=logLiks)) + geom_line() + labs(title="Evolution of log likelihood along with iterations", x="iterations",y="log likelihood") + theme(plot.title=element_text(vjust=0.5, hjust=0.5))

リスト1．モデルの構築

モデルの構築にあたってGibbsサンプリングを用い、トピック数は350とした。これはブログ記事「topicmodelsパッケージでLDA（パーツ作りました）」の図4で、トピック数が350のとき対数尤度の調和平均が最大になったことに基づいている。つまり、モデル選択の基準として尤度を最大にするという基準を用いている。
また、LDAのオプションであるburninやiterの妥当性を確かめるためにiterationが10, 20, ..., 1500に対する対数尤度を求めてプロットしたのが図1である。

図1．対数尤度のiteration依存性

図1から、対数尤度はiterationsが200あたりで急激な増加から緩やかな増加に転じており、500を超えたあたりではほぼ飽和しているとみなせるのでburnin=500, iter=1000という設定は妥当なものと考えられる。

LDAvisの結果

作成したモデルをLDAvisを使って可視化するスクリプトをリスト2に示す。

################################################
#
# LDAvisで使うデータを作成する関数
# 　fit: LDAモデル
# 　doc_term: DocumentTermMatrix
#
################################################
buildVisdata1 <- function(fit, doc_term) {
  phi <- posterior(fit)$terms %>% as.matrix
  theta <- posterior(fit)$topics %>% as.matrix
  vocab <- colnames(phi)
  doc.length <- rowSums(as.matrix(doc_term))
  term.frequency <- colSums(as.matrix(doc_term))
  params <- list(phi = phi,
                 theta = theta,
                 doc.length = doc.length,
                 vocab = vocab,
                 term.frequency = term.frequency)
  return(params)
}

install.packages("tsne")
library("tsne")
svd_tsne <- function(x) tsne(svd(x)$u)

install.packages('servr') 
library(servr)

# LDAvis
install.packages("LDAvis")
library(LDAvis)
VisualizeModel <- function(fit, doc_term, dir) {
  params <- buildVisdata1(fit, doc_term)
  json <- createJSON(phi = params$phi,
                     theta = params$theta,
                     doc.length = params$doc.length,
                     vocab = params$vocab,
                     term.frequency = params$term.frequency,
                     mds.method = svd_tsne
  )
  serVis(json, out.dir = dir, open.browser = TRUE)
}

VisualizeModel(model, dpc_dtm, './')

リスト2．

リスト2を実行した結果を図2に示す。

図2．LDAvisの実行結果（トピック数=350）

左図はトピックの全体像を示す。各円がトピックを表し、その面積はコーパス中におけるそのトピックの相対的な割合（prevalence）に比例して描かれている。各円に付けられた数字は割合の大きい順につけられた番号である。

右図は選択されたトピックに関連のある単語を横棒グラフで示したものである。灰色の棒はコーパス全体における各単語の出現頻度を表しており、赤色の棒は選択されているトピックに固有な出現頻度を表している。

図2で問題なのは

右図で単語（疾患）がコードで表されており、疾患名がわからない
トピックに関連する疾患がたった1つしかないように見える

点である。
まず１については、buildVisdata1でvocabへcolnames(phi)、すなわちposterior(fit)$termsのカラム名を代入しているので、疾患名ではなくDPCの疾患コードが設定され、当然の結果になっている。
ここに疾患名を出すのであれば、buildVisdata1でvocabへ疾患名リストを設定する必要がある。疾患名リストが疾患コード（DPCコード）と対応付けられていればそれを使えばよい。
まず、現在のvocabには

colnames(posterior(model)$terms)

が設定されている。
次に、疾患名と疾患コードの関係は、dpc_tibbleから(DPC2, 疾患名)を抽出してDPC2で圧縮すればよい。最後にDPC2をキーとして結合し、疾患名だけ抜けば疾患名リストが完成する。以上の処理をリスト3に示す。

disease_dic_tibble <- dpc_tibble %>%
  tidyr::unite(疾患手術名, c('疾患名', '手術')) %>%
  dplyr::select(DPC2, 疾患手術名) %>%
  dplyr::distinct(DPC2,.keep_all=TRUE) %>%
  dplyr::arrange(DPC2)
  
terms_tibble <- tbl_df(colnames(posterior(model)$terms))
colnames(terms_tibble) <- c("DPC2")
vocab_tibble <- dplyr::left_join(terms_tibble, disease_dic_tibble, by="DPC2")
disease_names <- vocab_tibble$疾患手術名

リスト3．疾患名リストの作成

こうして求めた疾患名リスト「disease_names」をリスト2のbuildVisdata1内でvocabへ設定して実行してみた。すると何も表示されない。そこで、VisualizeModelが作成したJSONデータ「lda.json」を見てみると、文字コードがShift_JISになっている。
そこで、index.htmlのcharsetをutf-8からShift_JISに変更してみた。すると今度は”d3がない”といった関係のないエラーが出る。どうも、Shift_JISではd3（グラフィックスライブラリ）はだめなようである。ということは、CSVファイルから作り直さなければならないということらしい。
WindowsのR-StudioはデフォルトでShift-JISになっているので、このあたりから設定しなおす必要がある（ああ大変だ・・・Pythonではこんな面倒はなかったのに・・・日本語が絡むと絶対に何か起きる）。

そこで、その場しのぎで「lda.json」の文字コードを手作業でShift-JISからutf-8に変換してみた。その結果を図2bに示す。

図2b．疾患名を表示したLDAvis

これは下記のURLから操作できる。

http://hinfokumw.html.xdomain.jp/lda/LDAvis/350/

トピックに含まれる疾患の数

まず、トピックごとの疾患の出現確率φを抽出する。それを行うのがposterior(model)$termsである。リスト4は各トピックがどのくらいの疾患を含んでいるか調べるスクリプトである。

phi <- posterior(model)$terms
threshold <- 1 / model@wordassignments$ncol # 1/V
terms_per_topic <- unlist(
  lapply(1:model@k, function(k) {
    return(sum(phi[k,] > threshold))
  })
)
df <- data.frame(value = terms_per_topic)
ggplot(df, aes(x = value)) + 
  geom_histogram(binwidth = 1) + 
  annotate("text", x=25,   y=40, label=paste("k =", model@k), hjust=0) +
  annotate("text", x=25,   y=38, label=paste("トピックに含まれる最大単語数 =", round(max(terms_per_topic), digit=1)), hjust=0) +
  annotate("text", x=25,   y=36, label=paste("トピックに含まれる最低単語数 =", round(min(terms_per_topic), digit=1)), hjust=0) +
  annotate("text", x=25,   y=34, label=paste("トピックに含まれる平均単語数 =", round(mean(terms_per_topic), digit=1)), hjust=0) +
  annotate("text", x=25,   y=32, label=paste("threshold = 1 / ", model@wordassignments$ncol), hjust=0) +
  labs(title = "トピックを構成する単語（疾患）の数", x = "単語の数", y = "頻度")

リスト4．トピックに含まれる疾患数の可視化スクリプト

posterior(model)$termsはトピック数×疾患数のmatrixで、(k, v)要素にはトピックkに含まれる疾患vの割合（φ_kv）が格納されている。
φ_kvの閾値thresholdとして1 / model@wordassignments$ncolを設定する。ここで、model@wordassignments$ncolには疾患数として801という値が入っているので、閾値はthreshold=1/801=0.00125になる。
上記のトピック×疾患matrixの要素がこの閾値を超えるものだけをそのトピックの構成疾患とみなして、各トピックがいくつの疾患で構成されているかを調べる。それが、terms_per_topicで、添え字がトピック番号、要素の値がそのトピックを構成する疾患数になっている。こうして求めたterms_per_topicをもとにヒストグラムを描いたのが図3である。

図3．トピックに含まれる疾患数

トピックに含まれる平均疾患数は12疾患で（中央値：7、最頻値：5）図に示すようにロングテールの分布を示す。

病院のトピック分布

同様にして各病院がいくつのトピックから構成されているかを求めてヒストグラムにするスクリプトをリスト5に、その結果を図4に示す。

theta <- posterior(model)$topics
threshold <- 1 / model@k
topics_per_document <- unlist(
  lapply(1:model@wordassignments$nrow, function(d) {
    return(sum(theta[d,] > threshold))
  })
)
df <- data.frame(value = topics_per_document)
ggplot(df, aes(x = value)) + 
  geom_histogram(binwidth = 1) + 
  annotate("text", x=25,   y=40, label=paste("k =", model@k), hjust=0) +
  annotate("text", x=25,   y=38, label=paste("病院が含む最小トピック数 =", round(min(topics_per_document), digit=1)), hjust=0) +
  annotate("text", x=25,   y=36, label=paste("病院が含む最大トピック数 =", round(max(topics_per_document), digit=1)), hjust=0) +
  annotate("text", x=25,   y=34, label=paste("病院が含む平均トピック数 =", round(mean(topics_per_document), digit=1)), hjust=0) +
  annotate("text", x=25,   y=32, label=paste("threshold = 1 / ", model@k), hjust=0) +
  labs(title = "病院に含まれるトピックの数", x = "トピックの数", y = "病院の数")

リスト5．病院のトピック分布を求めてヒストグラムにするスクリプト

閾値は1 / model@kとした。ここで、model@kはトピック数（350）である。

図4．病院のトピック分布

図4に示すように、病院が含むトピック数の最大は115で最小は3、平均は57である（中央値も57）。
図4から次のことがわかる。

図4はほぼ左右対称である
最大トピック数は350であるが、最も多いトピックを含む病院でも115トピックしか含んでいない

トピックを特徴づける疾患

各トピックに含まれる疾患の割合が閾値(1 / model@wordassignments$nrow=1/1664)を超えるものについて、その割合の降順に並べたリストを作成するスクリプトをリスト6に示す。

install.packages("radiant.data")
library(radiant.data)

phi_tibble <- dplyr::as_tibble(posterior(model)$terms, rownames=NA) %>%
  radiant.data::rownames_to_column('topic') %>%
  tidyr::gather(key = DPC2, value = phi, -topic)

phi_tibble$topic <- as.integer(phi_tibble$topic)

phi_with_disease_name <- phi_tibble %>%
  dplyr::filter(phi > 1 / model@wordassignments$nrow) %>%
  dplyr::arrange(topic, desc(phi)) %>%
  dplyr::left_join(disease_dic_tibble, by = "DPC2")

write.csv(phi_with_disease_name, "./phi_with_disease_name.csv")

リスト6．トピックに含まれる疾患をリストアップするスクリプト

トピック×疾患マトリックスposterior(model)$termsをtibble型に変換するためにas_tibbleを使っている。これは、tbl_dfだと行名が失われるからである。as_tibbleを使うと、rownames=NAとすることにより、行名を残すことができる（rownamesを省略すると行名は除去される）。
tibble型に変換したらrownames_to_columnによって行名を列に変換し（列名は引数に指定した'topic'）、gatherによってデータフレームの構造を"列型"から"行型"に変換する。ここで"列型"というのは、データが列方向に伸びているデータ構造のことで、列をバラバラに切り離して行方向にデータを並べたものを"行型"と呼んでいる。ちょうどリレーショナルデータベースにおいて繰り返しを取り除く第一正規化に似た処理である。
gatherの第一引数keyには”列型”の列名を値とする項目に付ける名前を指定する。ここではDPC2としている。第2引数valueには列名に対応する値に付ける名前を指定する。ここではphiとしている。最後の引数-topicは”列型”のデータフレームの列topic以外すべての列をバラバラにすることを意味している（図5）。

図5．dplyr::gather()関数のイメージ

リスト7にはリスト6の実行結果を示す。

> phi_with_disease_name
# A tibble: 5,173 x 4
   topic DPC2          phi 疾患手術名                       
   <int> <chr>       <dbl> <chr>                            
 1     1 05017002 0.631    閉塞性動脈疾患_02                
 2     1 05013099 0.323    心不全_99                        
 3     1 04008099 0.0150   肺炎等_99                        
 4     1 06003501 0.00226  結腸（虫垂を含む。）の悪性腫瘍_01
 5     1 06016001 0.00192  鼠径ヘルニア_01                  
 6     1 05005099 0.00140  狭心症、慢性虚血性心疾患_99      
 7     1 06002097 0.000881 胃の悪性腫瘍_97                  
 8     2 09001002 0.553    乳房の悪性腫瘍_02                
 9     2 09001005 0.333    乳房の悪性腫瘍_05                
10     2 06003501 0.0559   結腸（虫垂を含む。）の悪性腫瘍_01
# ... with 5,163 more rows
>

リスト7．各トピックを特徴づける疾患

医療分野でトピックモデルを利用した研究

はじめに

医療分野に限定してトピックモデル・・・なかでもLDA、およびその派生・・・を利用した文献を漁ってみた。
最初はCiNiiやPubMedで検索していたが、フリーで手に入るものは圧倒的にGoogleの検索結果が充実している。
トピックモデルと言えばテキストマイニングが多いが、医療分野ではカルテをマイニングするだけでなく、問診票や電子レセプトをLDAやその派生形を使って解析する例も見られ、大変興味深い。
今後少しずつ読んでまとめていきたい。

[1] 問診データに対する潜在トピックモデルに基づく健診データ解析

問診票をドキュメント、問診項目を単語、被験者の生活習慣を潜在トピックとしてLDAを用いてトピック分析を行ったという研究。
被験者を出現率が最大となるトピックで分類し、各被験者群の間で検査値データに違いがあるかどうかでモデルの評価を行っている。
一般に、トピックモデルは教師なし学習であるため評価が難しいとされている。さらにトピックを人間が解釈するのも難しいとされている。したがって何をもってモデルが妥当か評価するのは難しい。そこでよく行われるのが外部基準を使って適切に分類されているかを検証するという手法である。この研究では外部基準として検査値を使用している。つまり、トピックモデルで分類された被験者群の間で検査値に有意な違いがあるかどうかでモデルの妥当性を検証している。
しかし、そもそも「トピックモデルで分類された被験者」とは一体何だろう？トピックモデルは潜在トピックを抽出する手法であり、対象（この場合は被験者）を分類するものではない（べつに分類してもいいけど・・・）。しかし、潜在トピックに意味付けができない以上、それらがどのくらい含まれているかを議論しても意味がない。そこで、最も含まれる割合が高いトピックで被験者を分類している。また、そのようにしてもよいことをネットワークグラフでクラスタリングされた群と比較して検証している。ここで、ネットワークグラフを構築する上で被験者間の類似性をJensen-Shannon divergenceによって計算し、それがある閾値（ここでは0.05としている）より小さい場合にノード間（被験者間）にエッジを設けてネットワークグラフを作成している。クラスタリングはNewmannアルゴリズムを使って被験者を群に分け、それが最大トピック成分で分類したものと同等であることを示している。
この文献には様々な手法（Jensen-Shannon divergence, ネットワークグラフ理論, Newmannアルゴリズム, concordance）が用いられており、非常に示唆に富む。何より問診票の問診項目を単語に見立ててトピックモデルを適用しているのが秀逸である。考えてみれば、Recommendation分析も商品アイテムを単語とみなして解析を行っているのだから、それと同じと言えば同じであるが。

畠山豊, 宮野伊知郎, 片岡浩巳, 中島典昭, 渡部輝明, 奥原義保: 問診データに対する潜在トピックモデルに基づく健診データ解析. 医療情報学, 33(5), 2013.

Hatakeyama Y, Miyano I, Kataoka H, Nakajima N, Watabe T, Yasuda N, Okuhara Y: Use of a Latent Topic Model for Characteristic Extraction from Health Checkup Questionnaire Data. Methods Inf Med. 54(6), 515-21, 2015.

[2] A Most Resource-Consuming Disease Estimation Method from Electronic Claim Data Based on Labeled LDA

電子レセプトをドキュメント、診療行為を単語、疾患を潜在トピックとみなして、最も医療資源を使った疾患を推定するといった内容の研究。ただし、ここではただのLDAではなく"Labeled LDA"なるものを使っている。
Labeled LDAというものをよく知らないが、この論文によると「ラベル付きの文書コレクションを生成するためのプロセスを記述する確率的トピックモデル」と書いてある。
図1が論文中に示されたグラフィカルモデルである。

通常のLDAに比べてΦ→λ→θが付け加わっている。λはトピックの有無を示すインディケータで、Φはその事前分布ということである。この図でλとwにハッチが付けられているので、これら2つが観測されるデータなのであろう。
どのようにしてモデルを生成するのかこの論文を読んでもよくわからなかった。特にλの正体がわからない。ノーテーションも λ と思われる箇所が l になっていたりして、間違っているのか、それとも自分の理解不足なのかもわからない。
検証は、モデルがはじき出したトピック（すなわち疾患）が電子レセプトに人手でつけられた疾患（DPC病名と思われる）に一致するかどうかで行っている。ここで、「モデルがはじき出したトピック」は、単に構成比率（つまりθ_{d,k}）が最大の疾患というのではなく、料金（点数）が最大の疾患としている点がポイントである。「最も医療資源を投入した」というのを「最もコストのかかった」と解釈してのことだろう。
しかし、この「コスト（点数）」はどのようにしてモデルに取り込んだのであろうか？通常の文書を対象としたLDAでは単語の出現頻度が重みの役割を果たす。この研究では医療行為の点数を頻度の代わりに使っているのだろうか？そのあたりが読み取れない。
また、盛んに疾患（つまりλ？）が分かっていなくとも問題ないといった記述がみられるが、これがどういうことなのかわからない。ラベル付きLDAなんだからラベル（つまり疾患）がデータとして与えられなかったら機能しないだろうと思うのだが。
モデルの評価方法は文献[1]と同じで「正解」（人手で付けたDPC病名）とモデルが分類した結果を比較している。
評価指標はRecall（再現率）、Precision（精度）、F-value（F値）で、SVMとナイーブベイズと比較し、かなり良い結果が得られたと述べている。
ここでもトピックモデルを分類に使っているが、文献[1]と違って構成比率（つまりθ_{d,k}）が最大のトピックには「最も医療資源を使った疾患」という大義名分があるのでトピックモデルの使い方としては間違っていないと感じた。

Yasutaka Hatakeyama, Takahiro Ogawa, Hironori IKEDA, Miki Haseyama: A Most Resource-Consuming Disease Estimation Method from Electronic Claim Data Based on Labeled LDA. IEICE Transactions on Information and Systems E99.D(3):763-768, 2016.

【補足】

もう一度読み返してみた。前読んだ時分からなかったことがいくらか分かるようになった。通常のLDAではドキュメントがすべてのトピックを含みうるが，Labeled LDAではドキュメントごとに含みうるトピックが異なり，それが「ラベル」ということらしい。つまり，DPCデータ（請求データ）にはいくつかの疾患が割り当てられており（併存症）それが当該DPCデータにおける「ラベル」になる。Labeld LDAは，これらの疾患の構成比率（含有率）θ_{d,k}と疾患ごとの診療行為の割合β_{k,w}を求め，(8)式によって与えられた請求データにおける疾患kが診療行為w_{d,i}を含む確率を求める。さらに，当該請求データの診療行為が疾患kに投入されたものであるか否かのインディケータを(7)式から計算し，それを用いて(6)式から疾患kに費やしたコスト（点数）を計算している。これが最大になる疾患kが最も医療資源を投入した疾患というわけである。
(6)式の右辺は内積であることに注意する。point(w_{d})は請求dの診療行為に対応する点数を要素に持つ行ベクトルである。一方，I(k,w_{d})は，請求dの診療行為が当該疾患kに投入されたものかどうかを示すインディケータを要素とする列ベクトルである。これらの内積を計算することで当該請求において疾患kに投入された医療資源の点数が計算できる。

[3]Adverse Drug Reaction Prediction with Symbolic Latent Dirichlet Allocation

医薬品文書をドキュメント、薬物副作用用語（ADR terms）を単語、医薬品文書にまたがって頻繁に現れる薬物副作用用語の集合をトピックとみなして3種類のLDA類似モデルを構築し、副作用（ADR）を予測するといった内容の研究である。
副作用の予測を行うために、まず対象とする医薬品を特徴づけるADRトピック分布を求め、次に薬物構造の特徴（features）をADRトピック分布と関連付けるために予測モデルを構築する。その後、その医薬品に関連する副作用を、そのトピック分布を通じて予測することができる。
この論文では3つのモデルが検討されている。まず、基本となるモデルのグラフィカルモデルを図1に示す。

図1．ベースモデル（The base model）

ベースモデルはLDAをそのまま利用したものである。

次にドメイン知識を採り入れた正規化モデル（The regularized rodel）のグラフィカルモデルを図2に示す。

図2．正規化モデル（The regularized rodel）

ドメイン知識として階層ADRオントロジーシステムを利用する。これは上位のADRが下位ADRの抽象化バージョンになるように階層化されたオントロジーで、これに違反した場合にペナルティーを与えることによってドメイン知識を反映させるものである。
具体的には同じトピックを持つ子ADRが異なる親ADRを持つ場合にペナルティーを科す。

最後に混合入力モデル（The mixed input model）のグラフィカルモデルを図3に示す。

図3．混合入力モデル（The mixed input model）

これは、入力となるADR用語と薬剤構造の特徴（feature）を別々に学習する代わりにADR用語と薬剤構造の特徴（図3のx）の特徴の両方を医薬品文書の単語として扱う。これによって学習速度が大幅に向上する。

構築したモデルはADReCSという副作用データベースを用いて評価を行っている。そして他の手法（lasso、CC: canonical correlation analysis）とROCのAUC（Area Under Curve）やPR（Precision-Recallのことだろうか？）のAUCを使って比較している。
ということは、やはり外的な基準（教師データ）を用いて分類問題に帰着させているのだろうか・・・。このあたりが読み取れなかった。

Cao Xiao, Ping Zhang, W. Art Chaowalitwongse, Jianying Hu, Fei Wang: Adverse Drug Reaction Prediction with Symbolic Latent Dirichlet Allocation. Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence (AAAI-17), 2017.

topicmodelsパッケージでLDA（パーツ作りました）

topicmodelsパッケージでLDAを行う際のパーツを作ってみた。

DTM

分析対象のDTMです。DPCデータです。

> dpc_dtm
<<DocumentTermMatrix (documents: 1664, terms: 801)>>
Non-/sparse entries: 167793/1165071
Sparsity           : 87%
Maximal term length: 8
Weighting          : term frequency (tf)
> dpc_dtm$nrow
[1] 1664
> dpc_dtm$ncol
[1] 801

リスト1．分析対象のDPCデータのDTM

モデルの作成

トピック数を10から100まで10ずつ増やしながらモデルを作ってみました。リスト2はそのスクリプトです。

topic_nums = c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)

burnin = 500
iter = 1000
keep = 100

models <- lapply(topic_nums, function(topic_num){
  model <- LDA(dpc_dtm, k = topic_num, method = "Gibbs", control = list(burnin = burnin, iter = iter, keep = keep) )
  return(model)
})

リスト2．モデルの作成

モデルの保存

リスト3は作成したモデルをファイルに保存するスクリプトです。

save(models, file = "models-10-100.rda")

リスト3．モデルをファイルに保存

保存したモデルはloadで読み込むことができます。

load("models-10-100.rda")

リスト4．モデルの読み込み

対数尤度のトピック数依存性

対数尤度の調和平均を計算してトピック数依存性をグラフ表示するスクリプトです。

logLiks_harmonicMean <- lapply(models, function(model){
  logLiks <- model@logLiks[-c(1:(burnin/keep))]
  return(harmonicMean(logLiks))
})

LL <- data.frame(
  topic_nums = topic_nums,
  logLiks = unlist(logLiks_harmonicMean)
)
ggplot(LL, aes(x=topic_nums, y=logLiks)) + geom_line()

リスト5．対数尤度のトピック数依存性をグラフ表示

実行結果を図1に示します。対数尤度はトピック数の増加とともに増大します。とても収束しそうにありません。過適合しているのでしょうか。

図1．対数尤度の調和平均のトピック数依存性

対数尤度の変化

モデル作成時に対数尤度がどのように変化するかをトピック数ごとにグラフにするスクリプトです。これを見ればburninはどのくらい必要かがわかります。

logLiks_array <- lapply(models, function(model){
  logLiks <- model@logLiks
  return(logLiks)
})

iter_array <- lapply(topic_nums, function(topic){
  return(seq(keep,burnin+iter,keep))
})

topic_array <- lapply(topic_nums, function(topic){
  return(rep(topic, length(seq(keep,burnin+iter,keep))))
})

LL <- data.frame(
  topic = factor(unlist(topic_array)),
  iter = unlist(iter_array),
  logLiks = unlist(logLiks_array)
)

ggplot(data=LL, aes(x=iter, y=logLiks, colour=topic)) +
  geom_line()

リスト6．対数尤度の変化（トピック数別）

リスト6の実行結果を図2に示します。トピック数の増大に伴い対数尤度は大きくなっていきますが、その勢いは小さくなっているようです。また、収束に必要なburninもトピック数に応じて増えているように見えます。

図2．対数尤度の変化

トピックモデルで最適なトピック数を決定する方法

トピックモデルで最適なトピック数を決定する方法について調べた。

optimal_k

トピック数を変えながらtopicmodelsを使ってモデルを作成し、データから対数尤度の調和平均を求め、最適なトピック数を求める。
この方法を使ったサンプルがTopic Models Learning and R Resourcesにある。

#Control List
control <- list(burnin = 500, iter = 1000, keep = 100, seed = 2500)
(k <- optimal_k(dpc_dtm, 40, control = control))

リスト1. optimai_k

関数optimal_kが最適なkを求める関数である。ここで、dpc_dtmはDocumentTermMatrixである。これを実行すると最適なトピック数 k とともに、対数尤度の調和平均のトピック数依存性のフラフが表示される。
dpc_dtmのプロフィールを以下に示す。

<<DocumentTermMatrix (documents: 1664, terms: 801)>>
Non-/sparse entries: 167793/1165071
Sparsity           : 87%
Maximal term length: 8
Weighting          : term frequency (tf)

リスト2. DocumentTermMatrixの概要

これはDPC参加病院のみ（1664病院、疾患別手術801件）を対象として作成したDocumentTermMatrixである。

これに対してリスト1を実行したところ、全部で10時間25分かかった。

結果を図1に示す。

図1．optimal_kの実行結果

図から、対数尤度の調和平均は飽和していないことがわかる。もう少しトピック数を増やして調べる必要がある。
しかしながら、optimal_kは、トピック数を2～最大トピック数まで1ずつ増やしながら対数尤度の調和平均を計算しているため、非常に時間がかかる。せめてトピック数を間引きながら増やすことができればよいのだが、関数仕様にそのようなパラメータはない。

そこで、リスト1のcontrolにおけるburninとiterを調整できないか検討してみる。具体的には、LDAをいくつかのトピック数で実行して、対数尤度の推移をプロットし、どのくらいの繰り返し (iter) で飽和するか調べ、最適なburninとiterを決定する。リスト1では、burnin=500、iter=1000になっているが、これを少しでも小さくできれば、optimal_kをもっと高速にできるかもしれない。
これを行うために以下のRスクリプトを作成した。

fitted <- LDA(dpc_dtm, k = 40, method = "Gibbs", control = list(burnin = 500, iter = 1000, keep = 50) )
LL <- data.frame(
  topic_nums = seq(1, 1500, 50),
  logLiks = fitted@logLiks
)
ggplot(LL, aes(x=topic_nums, y=logLiks)) + geom_line()

リスト3. 対数尤度の推移

実行結果を図2に示す。

図2．対数尤度の推移（トピック数=40）

図から分かるように最低でもburnin=250, iter=500は必要で、burnin=500, iter=1000くらいはあった方が良いと考えられる。

したがって、burninとiterの最適化は諦めて、トピック数の間引きを考える。そのために、リスト4のようなRスクリプトを作成した。

topic_nums = seq(5, 100, 5)

burnin = 500
iter = 1000
keep = 100

models <- lapply(topic_nums, function(topic_num){
  model <- LDA(dpc_dtm, k = topic_num, method = "Gibbs", control = list(burnin = burnin, iter = iter, keep = keep) )
  return(model)
})
# モデルを保存
save(models, file = "models-b500-i1000-5-100-5.rda")

# 対数尤度の調和平均のトピック数依存性
logLiks_harmonicMean <- lapply(models, function(model){
  logLiks <- model@logLiks[-c(1:(burnin/keep))]
  return(harmonicMean(logLiks))
})

LL <- data.frame(
  topic_nums = topic_nums,
  logLiks = unlist(logLiks_harmonicMean)
)
ggplot(LL, aes(x=topic_nums, y=logLiks)) + geom_line()

リスト4．対数尤度のトピック数依存性（k=5,10,15,...,100）

リスト4はトピック数を5から始めて5ずつ増やしながら100まで対数尤度の調和平均を計算し、対数尤度の調和平均のトピック数依存性をグラフ化している。結果を図3に示す。

図3．対数尤度の調和平均のトピック数依存性（k=5~100）

図3を見るとトピック数が増えるにつれて対数尤度の調和平均は単調に増加しており、飽和しそうにない。
そこで、今度はトピック数を50から始めて50ずつ増やしながら500まで変化させて対数尤度の調和平均を求めてグラフにした。結果を図4に示す。

図4．対数尤度の調和平均のトピック数依存性（k=50~500）

図からトピック数が350のとき対数尤度の調和平均が最大になっていることがわかる。

図5は繰り返し数（iteration）が進むにつれて対数尤度がどのように変化するかをトピック数ごとに示したものである。

図5．対数尤度の変化

トピック数が大きくなるにつれて対数尤度が飽和するのに要するステップ数が増加する傾向が見られるものの、いずれも収束傾向にある。

Perplexity

これらの結果からトピック数を増やせば増やすほど対数尤度の調和平均は増加することが分かった。原因としてモデルがデータに過適合（オーバーフィッティング）していることが考えられる。

そこで、別の指標で最適なトピック数を求めてみる。ここでは、指標としてPerplexityを考える。ここでは、Topic models: cross validation with loglikelihood or perplexityに従ってPerplexityを計算する。

その前に、Perplexityとは何かについて説明する。
「エントロピーとパープレキシティ」によれば、Perplexityとは「単語の平均分岐数を表しており・・・大きいほど，単語の特定が難しく，言語として複雑になる」というものである。１単語あたりのエントロピーをHとした場合、Perplexityは2^Hになる。
そこで、さまざまなトピック数kに対してPerplexityを計算し、最小のPerplexityを与えるトピック数 k を求めてみる。
先述したTopic models: cross validation with loglikelihood or perplexity
には"Perplexity is a measure of how well a probability model fits a new set of data."という記述がある。つまり、Perplexityを計算するにはモデルを構築したデータ（すなわち訓練データ）の他に検証データが必要となる。

そこで、DTM形式のDPCデータの75%をランダムに抽出して訓練データとし、残りを検証データとする。

full_data  <- dpc_dtm
n <- nrow(full_data)

splitter <- sample(1:n, round(n * 0.75))
train_set <- full_data[splitter, ]
valid_set <- full_data[-splitter, ]

リスト5．全データを訓練データと検証データに分割（3:1）

次に、訓練データを用いてあるトピック数（ここでは10としている）に対してモデルを構築する。

topic_num = 10

burnin = 500
iter = 1000
keep = 100

model <- LDA(train_set, k = topic_num, method = "Gibbs", control = list(burnin = burnin, iter = iter, keep = keep) )

リスト6．モデルの構築

最後に訓練データと検証データに対してPerplexityを計算する。

PPL_train <- perplexity(model, newdata = train_set)
PPL_valid <- perplexity(model, newdata = valid_set)

リスト7. Perplexityの計算

これをあるトピック数から始めて少しずつ増やしながらPerplexityを計算し、最小を与えるトピック数が最適なトピック数と判断するという考え方である。
リスト8は、これを行うための"Using perplexity and cross-validation to determine a good number of topics"にあるスクリプトである。

install.packages("doParallel")
library(doParallel)
cluster <- makeCluster(detectCores(logical = TRUE) - 1) # leave one CPU spare...
registerDoParallel(cluster)

# load up the needed R package on all the parallel sessions
clusterEvalQ(cluster, {
  library(topicmodels)
})

folds <- 5
splitfolds <- sample(1:folds, n, replace = TRUE)
candidate_k <- c(10, 20, 30, 40, 50, 100, 200, 300, 400, 500) # candidates for how many topics

# export all the needed R objects to the parallel sessions
clusterExport(cluster, c("full_data", "burnin", "iter", "keep", "splitfolds", "folds", "candidate_k"))

# we parallelize by the different number of topics.  A processor is allocated a value
# of k, and does the cross-validation serially.  This is because it is assumed there
# are more candidate values of k than there are cross-validation folds, hence it
# will be more efficient to parallelise
system.time({
  results <- foreach(j = 1:length(candidate_k), .combine = rbind) %dopar%{
    k <- candidate_k[j]
    results_1k <- matrix(0, nrow = folds, ncol = 2)
    colnames(results_1k) <- c("k", "perplexity")
    for(i in 1:folds){
      train_set <- full_data[splitfolds != i , ]
      valid_set <- full_data[splitfolds == i, ]
      
      fitted <- LDA(train_set, k = k, method = "Gibbs",
                    control = list(burnin = burnin, iter = iter, keep = keep) )
      results_1k[i,] <- c(k, perplexity(fitted, newdata = valid_set))
    }
    return(results_1k)
  }
})
stopCluster(cluster)

results_df <- as.data.frame(results)

ggplot(results_df, aes(x = k, y = perplexity)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  ggtitle("5-fold cross-validation of topic modelling with the 'DPC' dataset",
          "(ie five different models fit for each candidate number of topics)") +
  labs(x = "Candidate number of topics", y = "Perplexity when fitting the trained model to the hold-out set")

リスト8．Perplexityから最適なトピック数を決定する

昨日から動かしているけどなかなか終わらない。やっと終わったところで経過時間を見るとなんと所用時間は241541.88（秒）・・・日にち換算で約2.8日。とてもintensiveな計算だ。結果を図6に示す。

図6．Perplexity（リスト8の結果）

トピック数が100あたりまでは急激に減少し、100を超えたあたりからはなだらかに減少し、500に到達しても最小には達せず、引き続き減少傾向にある。

これを図4（対数尤度の調和平均）と比較すると、必ずしも極値は一致しないように思われる。

調べているとこんな記事があった。なんでもldatuningというパッケージがあって、モデルからトピック数を選択するいくつかの指標をはじき出してくれるらしい。その中に「Griffiths2004」というのがあって「要するに、トピック数をT、コーパス全体の単語をwとして、対数尤度logP(w|T)が最大になるTにすればいいじゃん、というアイデアである。」ということらしい。これは、もしかしてリスト4でやっていることなんだろうか？グラフを見ていると似ているなぁと思う。
リスト9にトピック数を50から50ずつ増やしながら500まで変化させたときの各種指標を計算するldatuningのスクリプトを示す。

install.packages("ldatuning")
library("ldatuning")

result <- FindTopicsNumber(
  dpc_dtm,
  topics = seq(from = 50, to = 500, by = 50),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 77),
  mc.cores = 2L,
  verbose = TRUE
)

# 可視化
FindTopicsNumber_plot(result)

リスト9．ldatuningによる4つの指標の計算と可視化（k=50~500）

図7にリスト9の実行結果を示す。

図9．4つの指標のトピック数依存性

図9を見ると、Arun2010とCaoJuan2009はトピック数が250～350あたりが最適で、Griffiths2004はトピック数が450でピーク、そしてDeveaud2014はトピック数が100でピークになっている。

LDAの場合、HDP（Hierarchical Dirichlet Process）というのがあって、これを使うとトピック数の自動決定が可能だそうである。しかし、RのHDP実装は無いみたいで、その代わり、GensimがPythonで実装している。加えてRからGensimを呼び出して利用できるようである。

※その後調べているとRにも"R pkg for Hierarchical Dirichlet Process"というのがあるようだ。しかし、"Works on MacOS and Linux, but may not install on Windows."と書かれている。

【補足】

関数optimal_kが何をしているかを知りたいならばソースがGitHubに公開されている。

また、出典となった論文はFindibg Scientific Topicsである。
トピック数の決定については"The input parameters for using latent Dirichlet allocation"や"Topic models: cross validation with loglikelihood or perplexity"で議論されている。

登録: 投稿 (Atom)