廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ノンパラメトリック検定からの増減比較という二段構えの方法について

今日の川柳

就職してから一年ちょっとたちました。就職して以来、一応バイオインフォマティクスっぽいことをやっています。

それで多少は医学っぽい論文もよむようになりました。

今日は腸内細菌の分析手法について書きます。

最近(2000年代以降)の腸内細菌の研究のブレークスルーはメタゲノム解析というやつです。

メタゲノム解析とは試料(糞便など)に含まれるDNAの断片をまとめて回収してきてつなぎあわせ、データベースに照会することで菌種を特定する手法のことを指します。

これにより、培養が難しい細菌の情報も手に入れられるようになりました。

腸内細菌は流行のトピックで、さまざまな疾病と腸内細菌の関連を示唆する論文がたくさん出ている一方、「この病気は決定的に腸内細菌が原因だ」「この細菌で病気が治る」というようなキラー分析はまだ少ない、というかほとんどない印象です。

病気と細菌の関係を探すのに、いまのところたぶん最も幅広く使われてる統計手法がノンパラメトリック検定です。

つまり、疾病を持つグループと、持たないグループに分けて2群を比較し、ある細菌が多いとか少ないとかやります。

実にシンプルでオーソドックスな分析に見えるかもしれません。

でもぼくはこれがあんまり好きじゃないです。

ノンパラメトリック検定をして→次に増減比較、という二段構えのやり方が素直じゃないと思うからです。

ノンパラメトリック検定の多くは

F_1(x) = F_2(x)

という帰無仮説を置きます。グループ1(疾病を持つグループ)とグループ2(持たないグループ)の分布は等しい、という仮定です。

ぼくは増減を比較したいなら、

\mu_1 = \mu_2

グループ1(疾病を持つグループ)とグループ2(持たないグループ)の平均は等しい、という仮定を置いたほうがいいと思います。

ちょっとシミュレーションしてみます。

細菌のデータは人ごとの細菌の構成比として扱われることが多いです。

構成比(足して1になる正の実数)を生成する分布としてディリクレ分布というのがあるので、これを使うことにします。

グループ1のディリクレ分布のパラメータは全部1、グループ2のディリクレ分布のパラメータは全部10とします。

細菌は全部で10種類いるとし、このうちのある菌種Aのみを比較します。有意水準は5%としておきましょう。

この場合、菌種Aの構成比の平均はどちらのディリクレ分布でも0.1です。

まず上記の設定でコルモゴロフ・スミルノフ検定を1万回やった場合のp値の分布です。

帰無仮説は、

F_1(x) = F_2(x)

なので「正しく」棄却されることが多いです。

f:id:abrahamcow:20190730205314p:plain

一方t検定の場合、帰無仮説は2群(を正規分布としたとき)の平均が等しいことなので、「正しく」p値が一様分布に近い分布になっています。

f:id:abrahamcow:20190730205639p:plain

t検定の代替としてよく用いられるウィルコクソン検定の場合、こんな感じです。

f:id:abrahamcow:20190730205806p:plain

帰無仮説は棄却されることが多めです。

つまりなにが言いたいかというと、ノンパラの検定とパラメトリックな平均の差の検定では置いてる前提が違うということです。

でも例えばこんな図や、

f:id:abrahamcow:20190730210201j:plain
Potential of fecal microbiota for early-stage detection of colorectal cancer より。アスタリスクはウィルコクソン検定の結果が有意だったことを示す。)

こんな表
f:id:abrahamcow:20190730210423p:plain
Variation in Microbiome LPS Immunogenicity Contributes to Autoimmunity in Humans より。アスタリスクはコルモゴロフ・スミルノフ検定の結果が有意だったことを示す。)

を見ると、いつの間にか2群の分布が等しいという仮定の下で行われた検定が、2群の平均なり中央値なりが異なるという仮説を検討したかのように扱われてしまう。

「ノンパラメトリック検定をして→次に増減比較、という二段構えのやり方」は本当に本来知りたかったことについて考えたことになるのか、という点は一度きちんと考えるべきだと思います。

じゃあどうすればいいのか、という話はまたいずれ書きます。

R のコードを貼ります。

library(gtools)
set.seed(1)
wp <- numeric(10000)
ksp <- numeric(10000)
tp <- numeric(10000)
for(i in 1:10000){
  y1 <- rdirichlet(50,rep(1,10))
  y2 <- rdirichlet(50,rep(10,10))
  wp[i] <- wilcox.test(y1[,1],y2[,1])$p.value
  ksp[i] <- ks.test(y1[,1],y2[,1])$p.value
  tp[i] <- t.test(y1[,1],y2[,1])$p.value
}

hist(ksp,main="ks.test",xlab = "p-value")
hist(tp,main="t.test",xlab = "p-value")
hist(wp,main="wilcox.test",xlab = "p-value")

実験医学増刊 Vol.37 No.2 腸内細菌叢 健康と疾患を制御するエコシステム

実験医学増刊 Vol.37 No.2 腸内細菌叢 健康と疾患を制御するエコシステム