Metric 予測：SmartDrug は IQ を高めるか？

本投稿は John Kruschke 著 Doing Bayesian Data Analysis の第 16 章 "Metric-Predicted Variable on One or Two Groups" をもとにした。

Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan

John Kruschke

Academic Press
Amazon.co.jpで詳細を見る

本章までは、コインの表裏などの「0 か 1」の分析で、predicted y（被説明変数 y）が

y ~ bernoulli(μ)

であったが、ここでは predicted が「metric（数値）」

y ~ normal(μ, σ)

分析内容は「薬 "SmartDrug" は IQ 値を高めるのか？」を調べる。データは以下の通り。

> summary(myDataFrame)
     Score               Group   
 Min.   : 50.00   Placebo   :57  
 1st Qu.: 92.75   Smart Drug:63  
 Median :102.00                  
 Mean   :104.13                  
 3rd Qu.:112.00                  
 Max.   :208.00                  

左のヒストグラムで上段がプラシーボ（偽薬）、下段が SmartDrug 。右の Boxplot で左からプラシーボ、SmartDrug 。

> p1<-ggplot(d=myDataFrame,aes(x=Score)) + geom_histogram(fill="white",colour="black") + facet_grid(Group ~ .)
> p2<-ggplot(d=myDataFrame,aes(y=Score,x=Group)) + geom_boxplot()
> library(grid)
> grid.newpage()
> pushViewport(viewport(layout = grid.layout(1,2)))
> print(p1,vp=viewport(layout.pos.row=1,layout.pos.col = 1))
> print(p2,vp=viewport(layout.pos.row=1,layout.pos.col = 2))

グラフから SmartDrug の効果は、多少は見られるが、大きな違いとは言い難い。偽薬と "SmartDrug" の平均と標準偏差は

# Placebo
> d1<-(myDataFrame$Score[myDataFrame$Group=="Placebo"]); mean(d1);sd(d1)
[1] 100.0351
[1] 17.8945
# Smart Drug
> d2<-(myDataFrame$Score[myDataFrame$Group=="Smart Drug"]); mean(d2);sd(d2)
[1] 107.8413
[1] 25.4452


MCMC 近似

このデータを以下のような「階層モデル」（図16.2）で考えることができる。

実行スクリプトは、Jags-Ymet-Xnom1grp-Mnormal.R に定義されている以下の「uniform版」（図16.2 の右のモデル）を使用。

  model {
    for ( i in 1:Ntotal ) {
      y[i] ~ dnorm( mu , 1/sigma^2 )
    }
    mu ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma ~ dunif( sdY/1000 , sdY*1000 )
  }

次は、Jags-Ymget-Xnom1grp-Mnormal-Example.R 実行結果例（本書図 16.3 に相当）。

これは、データで Group が "Smart Drug" を MCMC サンプリングした結果。IQ の平均 100、標準偏差 15 を ROPE に使用して、Smart Drug の効果を測る。

上段左図：95%HDI[101, 114] は ROPE[99,101] を僅かに外している程度で、効果を断定しにくい。

下段右図：Effect Size とは the "standardized" change（標準化された変化）。ここでは(μ-100)/σ で求める。95%HDI[0.0492, 0.551] は ROPE[-0.1, 0.1]（心理学調査において、通常「小さい」effect size は 0.2 ）を外していない。つまり「小さい effect size ですら無い」ということ。

下段左図：標準偏差、すなわちデータのばらつきは SmartDrug では明確な違いが出ている。一般的な IQ の標準偏差は 15 なので ROPE[14, 16] と設定すると、95%HDI[21.5, 30.8] は明らかに ROPE を外している。SmartDrug の服用は、IQ を良い方にも悪い方にも働いたと見ることができる。

下段左図：これは posterior-predictive check 。データのヒストグラムに、MCMC のサンプリングで得た正規分布を重ねたもの。今回の観測データは標本数が少ないため評価は難しいが、データの「外れ値」をモデルが捉えようとしている、などのヒントは得られる。

この「外れ値の対応」について、次回「外れ値と Robust 推定（ t 分布）」で行う。