murawaki の雑記

はてなグループから移転してきました

STRUCTURE と ADMIXTURE の混合分布モデル

久しぶりに NLP タグをつけたが、生物系の混合分布モデルの話。ゲノムを使う系統関係の論文では、PCA と並んで、よくこういう図が出てくる。

参考までに caption も引用。

(A) Representative estimate of population structure for 1,384 individuals from worldwide populations, including 432 individuals from India. The plot represents the highest-likelihood run among ten STRUCTURE runs with K = 7 clusters. Eight of the other nine runs identified a cluster largely corresponding to India, and five of these eight produced plots nearly identical to the one shown.

Figure 2. Population Structure Inferred from Microsatellite and Insertion/Deletion Polymorphisms

縦の列 (column) が各個体。各個体はゲノムの列 (sequence) で表現されている。これが K=7 個の潜在クラスで色分けされている。要は混合分布モデル、というか NLP 業界的にいうとトピックモデル。となると、具体的にどういうモデルなのか気になるところ。しかし、論文を読むと、生物系の人が生物の言葉で語っていて何度かくじけた。今回 ADMIXTURE の論文 (2009) を見たところ、最初から統計の言葉で説明されていて、ようやく糸口がつかめた。NLP 的な説明に翻訳してみる。

まずソフトウェアの確認から。STRUCTURE という検索泣かせな名前のソフトが昔からあった。最近、ADMIXTURE というこれまた嫌がらせのようなソフトが出てきた。新しい論文では ADMIXTURE を使っていることが多い。他に frappe というソフトもあるが、それほど見かけない。まずは新しい方の ADMIXTURE を見て、次に STRUCTURE に移る。

ADMIXTURE の混合分布モデルのグラフィカルモデルは以下の通り。

  • 事前分布が設定されておらず、pLSI 的。
  • 3 重の plate になっている。外側の I が個体のループ。次の J が DNA の列のループ。言語のトピックモデルだとこの 2 つ (I: 文書, J: 文書内の単語)。A は染色体の数。最近の genome-wide SNP の話だと、diploid といって、両親から 1 個ずつ受け継ぐため、A = 2 らしい。
  • \thetaは個体ごとの混合比。要素数は K。結果の図で色分けされているのはこれ。
  • \varphi が K と J の 2 重ループになっているのも特徴的。言語のトピックモデルだと K ごとにサイズ V の語彙分布を持っている。DNA の場合は列の場所ごとに別の分布を持っているので K x J 個の変数が必要。SNP の場合はベルヌーイ分布。
  • 記号は言語のトピックモデル風に変更している。また、元の説明だとカウントの分布 (多項分布) を考えているが、ここでは列の分布 (categorical 分布) を示している。

Z で周辺化して、W の確率にすると以下の通り。
\begin{eqnarray} p(W | \Theta,\Phi) &=& \prod_i \prod_j \prod_a \sum_k p(z_{i,j,a}=k | \theta_i) p(w_{i,j,a} | z_{i,j,a}=k, \Phi)\\ &=& \prod_i \prod_j \prod_a \sum_k \theta_{i, k} \,\times\, \varphi_{j,k,w_{i,j,a}} \end{eqnarray}
推論は、論文ではまず EM を導入する。しかし EM は遅いからと、別の手法を提案する。EM で遅いと言われると、サンプリング脳なのでつらい。

次。STRUCTURE のグラフィカルモデルは、Pritchard et al. (2000) によると以下の通り。

ADMIXTURE のモデルとの違いは、事前分布が追加されていること。\alpha\eta はいずれも Dirichlet 分布のパラメータ。symmetric なパラメータを一つ与えるか、経験ベイズ的にデータから推定するかでモデルに変種がある。ほぼ LDA。

推論。\theta\varphi は共役性を利用して積分消去したいところだが、元論文はそのままにしている。\theta\varphi と z を (実は\alphaも) MCMC でサンプリングする。

欠損値は、ADMIXTURE の場合、あらかじめ補完するという。STRUCTURE のような MCMC であれば、補完を sampling に組み込むのは簡単そう。

トピック数 K はあらかじめ指定する。Pritchard et al. (2000) では K を自動推定する怪しげなモデルが説明されている。実際に使われているのだろうか。AIC などを使ってモデル選択をするという手もある。論文でよく見かけるのは、K = 2 ... 5 くらいの結果を並べてお茶を濁すもの。

新しい ADMIXTURE の方がモデルが退化しているのが妙なところ。STRUCTURE はサンプリングの遅さが嫌われて ADMIXTURE への移行が進んでいるみたい。規模感としては、I が千ぐらい、J が数十万。確かに小規模とはいえない。でも、Wikipedia の記事 3M ページに対するトピック推定などと比べると、特別大きいわけでもない。

似た研究を別々に進めるのは不健全。LDA を提案した Blei et al. (NIPS2002) が 2002 年だから、実は STRUCTURE の Pritchard et al. (2000) の方が先行している。NIPS 2002 でも、2003 年の JMLR 版でも、Pritchard et al. (2000) への言及がない。2004 年の Blei の博論では引用されているので、このあたりで生物系の研究に気付いたらしい。というか、Blei の論文リストを眺めていると、2015 年になって Posterior predictive checks to quantify lack-of-fit in admixture models of latent population structure という論文を出しているのに気付いた。

ADMIXTURE の論文は 2009 年に出ているが、トピックモデルへの言及がない。ここ 10 年ぐらいで発展したトピックモデルの手法が DNA データにもそのまま使えそう。例えば、階層 Dirichlet 過程を使ってトピック数 K をデータに決めさせるとか、高速化の手法とか。需要はないのだろうか。