廿TT

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

対応のある標本の分析をわざわざStanでやる

今日の川柳

sleep データ

モデルに個体差のパラメータを入れるかどうかのポイントは、「個体差や場所差が識別できてしまうようなデータのとりかたをしているか」だという話があります(久保 (2012)、p.161)。

「個体差が識別できてしまうデータ」にはたとえば R に最初から入ってる sleep データがあります。

library(tidyverse)
p1<-ggplot(sleep,aes(x=group,y=extra,group=ID))+
  geom_point()+
  geom_line()+
  geom_text(data=sleep[sleep$group==2,],aes(label=ID),hjust=-1)+
  theme_bw(14)
print(p1)

f:id:abrahamcow:20180223211103p:plain

これはおなじ被験者に偽薬(1)と睡眠薬(2)を与えて、睡眠時間ののびを調べたデータです。

上図を見ると、睡眠薬(2)を飲んだほうが睡眠時間が伸びそうです。

しかし、平均の差の検定(ウェルチの t 検定)をやると5%水準では有意になりません。

> with(sleep,t.test(extra[group == 2],extra[group == 1]))

	Welch Two Sample t-test

data:  extra[group == 2] and extra[group == 1]
t = 1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2054832  3.3654832
sample estimates:
mean of x mean of y 
     2.33      0.75 

そこで、個体差を入れたモデルを考えます。

y_i \sim N(\beta_0+\beta_1 x_i +w_{ID_i},\sigma^2)

w_j \sim N(0,\tau^2)

y_i は睡眠時間の伸び、x_i睡眠薬(2)を飲んだか否か、 ID_i は被験者のIDを表す変数です。

Stan で書くとこうなります。

//mod_sleep.stan
data{
  int<lower=1> n;
  int<lower=1,upper=n> ID[2*n];
  real y[2*n];
  int<lower=0,upper=1> x[2*n];
}
parameters{
  real beta[2];
  real w[n];
  real<lower=0> sigma;
  real<lower=0> tau;
}
model{
  for(i in 1:(2*n)){
    y[i] ~ normal(beta[1]+beta[2]*x[i]+w[ID[i]],sigma);
  }
  w ~ normal(0,tau);
}
library(rstan)
rstan_options(auto_write = TRUE)
mod_sleep <-stan_model("mod_sleep.stan")
datlist <- list(n=n_distinct(sleep$ID),y=sleep$extra,
                x=abs(as.integer(sleep$group)-2),
                ID=as.integer(sleep$ID))
fit_sleep <-sampling(mod_sleep,datlist)
traceplot(fit_sleep,pars=c("sigma","tau","beta","w"))

無事収束しました。

f:id:abrahamcow:20180223212753p:plain

> round(summary(fit_sleep,pars="beta[2]")$summary[,c("mean","2.5%","97.5%")],2)
 mean  2.5% 97.5% 
 1.58  0.63  2.54 

係数 \beta_1 の95%信用区間は 0 から離れており、薬には効果があると言えそうです。

賢明なあなたはすでにお気づきだと思いますが、こんな面倒なことをしなくても、対応のある t 検定をやれば薬の効果は見れます。

> with(sleep,t.test(extra[group == 2],extra[group == 1], paired = TRUE))

	Paired t-test

data:  extra[group == 2] and extra[group == 1]
t = 4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.7001142 2.4598858
sample estimates:
mean of the differences 
                   1.58 

だいたい似たような結果になりました。

でも、対応のある検定、繰り返し測定の数が増えてくると面倒じゃないですか?

次の例にいきます。

Orthodont データ

繰り返し測定の回数が3回以上の例として、Orthodont データを取り上げます。

data(Orthodont,package = "nlme")
p2 <-ggplot(Orthodont,aes(x=age))+
  geom_line(aes(y=distance,group=Subject,colour=Sex))+
  theme_bw(14)
print(p2)

f:id:abrahamcow:20180223214920p:plain

これは歯の矯正かなんかのデータです。2年おきに4回、おなじ被験者から測定されています。

混合モデルとか使うのが王道かもしれませんが、混合モデルは固定効果とかランダム効果とか難しい言葉が出てきて混乱する。

以下のモデルを当てはめます。

y_i \sim N(\alpha_{age_i}+\beta x_i +w_{ID_i},\sigma^2)

w_j \sim N(0,\tau^2)

x_i は性別(0:男性 or 1:女性)、age_i は年齢を示す変数です。

Stan で書くとこうなります。

//modrepmes.stan
data{
  int<lower=1> l;
  int<lower=1> n;
  int<lower=1> m;
  int<lower=0,upper=l> age[n];
  int<lower=1,upper=n> ID[n];
  int<lower=0,upper=1> x[n];
  real y[n];
}
parameters{
  real w[m];
  real alpha[l];
  real beta;
  real<lower=0> sigma;
  real<lower=0> tau;
}
model{
  w ~ normal(0,tau);
  for(i in 1:n){
   y[i] ~ normal(alpha[age[i]]+beta*x[i]+w[ID[i]],sigma); 
  }
}
modrepmes <-stan_model("modrepmes.stan")
datlist <- list(n=nrow(Orthodont),m=n_distinct(Orthodont$Subject),
                l=n_distinct(Orthodont$age),
                age=as.integer(factor(Orthodont$age)),
                ID=as.integer(Orthodont$Subject),
                y=Orthodont$distance,
                x=as.integer(Orthodont$Sex)-1)
fitOrth <- sampling(modrepmes,datlist)
traceplot(fitOrth,pars=c("beta","alpha","sigma","tau","w"))

無事収束しました。

f:id:abrahamcow:20180223221822p:plain

性別の効果 \beta の事後分布を見てみます。

ex <- rstan::extract(fitOrth)
hist(ex$beta,xlab = expression(beta),main="",breaks ="FD")

f:id:abrahamcow:20180223221832p:plain

はい。

年齢ごとのベースラインとなる  \alpha_j の事後期待値と95%信用区間を図に示します。

CI_alpha <-t(apply(ex$alpha,2,quantile,prob=c(0.025,0.975))) %>% 
  as_data_frame() %>% 
  mutate(age=seq(8,14,by=2),
         mean =get_posterior_mean(fitOrth,"alpha")[,"mean-all chains"])


p2_2 <-p2+geom_line(data=CI_alpha,aes(y=mean),colour="royalblue",size=1,linetype=2)+
  geom_ribbon(data=CI_alpha,aes(ymin=CI_alpha[,1],ymax=CI_alpha[,2]),
              fill="royalblue",alpha=0.3)
print(p2_2)

f:id:abrahamcow:20180223222122p:plain

はい。

つづくかも。