JAGSで2群の差の推定

JAGSで2群の差の推定を行ったので、備忘録を残しておきます。

参考にしたサイトは以下です。
例によって、コードは間違っていたらご助言ください^^;。
model {#prior 
mu.y ~ dnorm(0, 1.0e-8) 
tau.y <- 1/(sigma.y*sigma.y) 
sigma.y ~ dunif(0,1.0E+8) 
mu.x ~ dnorm(0, 1.0e-8) 
tau.x <- 1/(sigma.x*sigma.x) 
sigma.x ~ dunif(0,1.0E+8) 

#likelihood 
for(i in 1:N.y){ 
y[i] ~ dnorm(mu.y, tau.y)
} 
for(j in 1:N.x){ 
x[j] ~ dnorm(mu.x, tau.x) 
} 

#generated quantities 
delta <- mu.y - mu.x 
cohen_d <- delta /((sigma.y^2*(N.y-1)+sigma.x^2*(N.x-1)) /((N.y-1)+(N.x-1)))^0.5 }
このモデルはtTest1.txtと名づけて保存します。

差の推定はgenerated quantitiesで表し、また効果量も推定できるようにしています。

このモデルは欠損値があってもなくても、同じように使えました。

Stanの場合、欠損値ありなしで以下のようにモデルが変わります。

あっ、上記のJAGSコードは一様分布ですけど、下記のStanコードはコーシー分布ですw。
#欠損値なし
data{
 int N1;
 int N2;
 real x[N1];
 real y[N1]; 
} 
parameters{
 real x_mu;
 real<lower=0> x_sigma;
 real y_mu;
 real<lower=0> y_sigma; 
} 
model{ 
 x_mu ~ normal(0, 100); 
 x_sigma ~ cauchy(0, 5);
 y_mu ~ normal(0, 100); 
 y_sigma ~ cauchy(0, 5); 
 x ~ normal(x_mu, x_sigma); 
 y ~ normal(y_mu, y_sigma); 
} 
generated quantities{
 real diff;
 diff = y_mu - x_mu; 
}
#欠損値あり
data{
 int x_N;
 int y_N;
 int x_Nmis;
 int y_Nmis;
 real x[x_N];
 real y[y_N]; 
} 
parameters{
 real x_mu;
 real<lower=0> x_sigma;
 real y_mu;
 real<lower=0> y_sigma;
 real x_mis[x_Nmis];
 real y_mis[y_Nmis]; 
} 
model{
 x_mu ~ normal(0, 100);
 x_sigma ~ cauchy(0, 5);
 y_mu ~ normal(0, 100);
 y_sigma ~ cauchy(0, 5);
 x ~ normal(x_mu, x_sigma);
 y ~ normal(y_mu, y_sigma);
 x_mis ~ normal(x_mu, x_sigma);
 y_mis ~ normal(y_mu, y_sigma); 
} 
generated quantities{
 real diff;
 diff = y_mu - x_mu; 
} 
JAGSなら欠損値があっても同じモデルでそのまま扱えるので良いですね。

秋期は、大学院修士課程の講義でベイズ統計モデリングを教えるのですが、欠損値の扱いやすさを考えると、StanよりもJAGSのほうがよいかもなぁと思いはじめました。




1.欠損値なしデータの場合

まずシミュレーションデータを作成し、普通にt検定です。
set.seed(123) 
x <- rnorm(60,10,4) 
y <- rnorm(40,15,7) 
summary(x)
t.test(x,y)
## 
## Welch Two Sample t-test 
## 
## data: x and y 
## t = -4.99, df = 55.486, p-value = 6.314e-06 
## alternative hypothesis: true difference in means is not equal to 0 
## 95 percent confidence interval: 
## -7.891531 -3.369775 
## sample estimates: 
## mean of x mean of y 
## 10.26247 15.89312

結果はxとyに有意差ありで、95%信頼区間も0をまたいでいません。

にJAGSで差を推定し、ggmcmcでプロットの確認です。

結論を先にいっておくと、Rhatは収束を示し、差に関してはJAGSのモデルでも同様の結果であり、加えて効果量が確認できます。

library(R2jags)
data <- list(y=y, N.y=40, x=x, N.x=60) fit <- jags(data=data, model.file="tTest1.txt", parameters.to.save=c("mu.y", "mu.x", "delta", "cohen_d"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit)
## Inference for Bugs model at "tTest1.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved 
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## cohen_d 1.102 0.243 0.616 0.944 1.100 1.260 1.591 1.001 
## delta 5.614 1.157 3.343 4.838 5.599 6.370 7.875 1.001 
## mu.x 10.263 0.473 9.345 9.940 10.270 10.580 11.212 1.001 
## mu.y 15.876 1.066 13.801 15.153 15.869 16.592 18.017 1.001 
## deviance 590.542 2.923 586.918 588.404 589.852 591.953 597.903 1.001 
## n.eff 
## cohen_d 2700 
## delta 3200 
## mu.x 4000 
## mu.y 4000 
## deviance 4000 
## 
## For each parameter, n.eff is a crude measure of effective sample size, 
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 4.3 and DIC = 594.8 
## DIC is an estimate of expected predictive error (lower deviance is better).
library(ggmcmc)
mcmc <- as.mcmc(fit) 
ggs.result <- ggs(mcmc) 
ggs_histogram(ggs.result)
ggs_density(ggs.result)ggs_traceplot(ggs.result)ggs_running(ggs.result) 
NewImage
NewImage
NewImage
NewImage
2.欠損値ありデータの場合
続いて、欠損値ありデータの例も示しておきます。
まずシミュレーションデータの生成です。
基本的な設定は上記と同じで、x1で10、y1で7の欠損値をつくっています。
set.seed(123) 
x1 <- rnorm(60,10,4) 
y1 <- rnorm(40,15,7) 
x1[sample(60, 10)] <- NA 
y1[sample(40, 7)] <- NA 
summary(x1)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 
## 2.134 8.260 10.360 10.290 12.310 18.680 10
summary(y1)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 
## 6.455 11.560 16.270 15.660 18.840 25.730 7
t.test(x1,y1)
## 
## Welch Two Sample t-test 
## 
## data: x1 and y1 
## t = -5.1148, df = 50.73, p-value = 4.845e-06 
## alternative hypothesis: true difference in means is not equal to 0 
## 95 percent confidence interval: 
## -7.480277 -3.262959 
## sample estimates: 
## mean of x mean of y 
## 10.28908 15.66070
結果は先と同じで、xとyに有意差ありで95%信頼区間も0をまたいでいません。
次にJAGSで差を推定します。
plotは割愛します。
使用するJAGSコードは、上記のtTest1.txtです。
library(R2jags) data1 <- list(y=y1, N.y=40, x=x1, N.x=60) fit1 <- jags(data=data, model.file="tTest1.txt", parameters.to.save=c("mu.y", "mu.x", "delta", "cohen_d"), n.chains=4, n.iter=2000, n.burnin=1000, DIC=TRUE)
print(fit1)
## Inference for Bugs model at "tTest1.txt", fit using jags, 
## 4 chains, each with 2000 iterations (first 1000 discarded) 
## n.sims = 4000 iterations saved 
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat 
## cohen_d 1.106 0.243 0.628 0.941 1.105 1.269 1.581 1.001 
## delta 5.645 1.174 3.327 4.842 5.643 6.433 8.025 1.001 
## mu.x 10.262 0.488 9.314 9.934 10.267 10.580 11.217 1.001 
## mu.y 15.907 1.076 13.774 15.196 15.907 16.616 18.013 1.001 
## deviance 590.656 3.048 586.925 588.398 589.951 592.204 598.202 1.001 
## n.eff 
## cohen_d 4000 
## delta 4000 
## mu.x 4000 
## mu.y 4000 
## deviance 3400 
## 
## For each parameter, n.eff is a crude measure of effective sample size, 
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## 
## DIC info (using the rule, pD = var(deviance)/2) 
## pD = 4.6 and DIC = 595.3 
## DIC is an estimate of expected predictive error (lower deviance is better).
欠損値なしデータの場合とほとんど同じ結果になりました(めでたしめでたし。