TJO氏にふられた質問で、反復測定分散分析（repeated measures anova）をベイズ推定に置き換えられないかと言う話があったので、質問文から想定される状況を整理し、欠損値補定をFully Bayesian Imputation Method *1で、モデル比較をBayes Factorで行なうコードを描いて試してみました。
結論はRのlme4パッケージあたりで時点ダミー×介入群ダミーを入れた線形混合モデルを推定した方がコード量も計算時間も1/100で済むといった感じなのですが *2、各時点各群の平均と分散をそれぞれ推定できること、欠損値の補定方法のパラメーターも同時推定できる利点もあるので、ファンがうぃーんと言うのが好きな人には良いかもです。うぃーん。
真のモデルは以下の表の通りです。
| 群 |
パラメーター |
t=1 |
t=2 |
t=3 |
| C |
 |
0.0 |
0.0 |
0.0 |
| C |
 |
1.0 |
1.0 |
1.0 |
| T |
 |
0.0 |
1.0 |
2.0 |
| T |
 |
1.0 |
1.1 |
1.2 |
t=2時点からT群とC群に差があるとしています。
この設定を元に、架空のデータを生成します。練習感あふれてしまいますが、真のパラメーターが分かることから、推定結果が妥当なものか検討できるので、手法の確認には実データを使うより良いとも言えます。
set.seed(604)
rm(list=ls())
n <- 120
t <- 1:3
g <- c("C", "T")
dimnames <- list(g, sprintf("t=%d", t))
mu <- matrix(c(0+0*t, 0+1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
sigma <- matrix(c(1+0*t, 1+0.1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
beta <- 1
x <- rep(runif(n), length(t))
df01 <- data.frame(
id=rep(1:n, length(t)),
time=rep(t, each=n),
group=rep(rep(g, each=n/length(g)), length(t)),
x=x,
z=x/2 + rnorm(n, sd=1),
score=numeric(n*length(t))
)
rm(x)
for(j in t){
for(i in 1:length(g)){
s <- n*(j-1)+n/length(g)*(i-1) + 1
e <- n*(j-1)+n/length(g)*i
df01$score[s:e] <- rnorm(n/length(g), mean=mu[i, j], sd=sigma[i, j]) + beta*df01$x[s:e]
}
}
df01$x[df01$id %in% sample(1:n, n/2)] <- NA
こんな感じのデータセットが出来上がります。
| id |
time |
group |
x |
z |
score |
| 1 |
1 |
C |
0.03264839 |
-0.39328775 |
0.1545386 |
| 2 |
1 |
C |
0.07093954 |
-0.07395381 |
-0.3282276 |
| 3 |
1 |
C |
0.96846472 |
0.70290536 |
1.6737583 |
| 4 |
1 |
C |
NA |
-0.79364909 |
1.0485097 |
| 5 |
1 |
C |
0.40517678 |
-0.50145817 |
0.7247011 |
| 6 |
1 |
C |
0.09162524 |
-1.43459105 |
0.5523630 |
| … |
… |
… |
… |
… |
… |
3. 推定コード
ベイズファクターの計算の都合上、推定はこれでもかと言うぐらいループが繰り返されます。
library(MCMCpack)
get_density <- function(m, samp, n=200){
max <- max(samp)
min <- min(samp)
range <- (max - min)/n
sum(samp > (m - range) & samp < (m + range))/length(samp)/2/range
}
get_posterior_density <- function(objf, init_theta, post.samp, estimated, log=TRUE){
pdensity <- numeric(length(init_theta))
len <- length(estimated)
pdensity[len] <- get_density(estimated[len], post.samp[,len])
for(i in 1:(len - 1)){
l <- len - i
post.samp.cnd <- MCMCmetrop1R(objf, theta.init=init_theta[1:l], force.samp=TRUE, estimated=estimated)
pdensity[l] <- get_density(estimated[l], post.samp.cnd[,l])
}
ifelse(log, sum(log(pdensity)), prod(pdensity))
}
theta2param <- function(theta){
theta <- theta[-(1:4)]
mlen <- length(t)*length(g)*2
tlen <- length(theta)
slen <- (mlen - tlen)/2
if(0 < slen){
len <- tlen/2
m <- c(theta[1:length(t)], theta[1:slen])
s <- c(theta[(1+len):(length(t)+len)], theta[(1+len):(slen+len)])
if((length(t)+1)<=len){
m <- c(m, theta[(length(t)+1):len])
s <- c(s, theta[(length(t)+1+len):(2*len)])
}
theta <- c(m, s)
}
len <- length(theta)/2
mu <- matrix(theta[1:len], ncol=length(t), byrow=TRUE, dimnames=dimnames)
sigma <- matrix(theta[(len+1):(2*len)], ncol=length(t), byrow=TRUE, dimname=dimnames)
list(mu=mu, sigma=sigma)
}
llf.nd <- function(theta){
p <- theta2param(theta)
beta <- theta[1]
alpha <- theta[2:4]
theta <- theta[-(1:4)]
len <- length(theta)/2
xi <- df01$x
s <- sum(dnorm(df01$x[!is.na(xi)] - alpha[1] - alpha[2]*df01$z[!is.na(xi)], sd=alpha[3], log=TRUE))
xi[is.na(xi)] <- alpha[1] + alpha[2]*df01$z[is.na(xi)]
for(j in t){
for(i in g){
s <- s + sum(dnorm(with(df01, { score[group==i & time==j] - beta*xi[group==i & time==j] })
, mean = p$mu[i, j]
, sd = p$sigma[i, j]
, log = TRUE))
}
}
s
}
prior.nd <- function(theta){
beta <- theta[1]
alpha <- theta[2:4]
theta <- theta[-(1:4)]
len <- length(theta)/2
mu <- 1:len
sigma <- (1+len):(2*len)
sum(dnorm(mu, mean=1, sd=1, log=TRUE)
,dgamma(sigma, shape=1.5, scale=2, log=TRUE)
,dnorm(beta, mean=1, sd=1, log=TRUE)
,dnorm(alpha[1], mean=0, sd=1, log=TRUE)
,dnorm(alpha[2], mean=0, sd=1, log=TRUE)
,dgamma(alpha[3], shape=1.5, scale=2, log=TRUE))
}
objf.nd <- function(theta, estimated=theta){
len <- length(estimated)
d <- len - length(theta)
if(0 < d){
p <- (len-d+1):len
theta[p] <- estimated[p]
}
llf.nd(theta) + prior.nd(theta)
}
calcPosterior <- function(l=16){
MCMCmetrop1R(objf.nd, theta.init=rep(1, l))
}
getMerginalLikelihood <- function(post.nd, l){
mean.nd <- summary(post.nd)$statistics[, "Mean"]
pdensity.nd <- get_posterior_density(objf.nd, rep(1, l), post.nd, mean.nd)
llf.nd(mean.nd) - pdensity.nd
}
ptrn <- seq(2*max(t), 0, -2) + 10
library(parallel)
library(doParallel)
library(foreach)
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)
r <- foreach(i=ptrn, .packages = c("MCMCpack")) %dopar% {
post.name <- "posterior"
mll.name <- "marginal likelihood"
r <- list(i, calcPosterior(i))
names(r) <- c("nop", post.name)
r[mll.name] <- getMerginalLikelihood(r[[post.name]], i)
return(r)
}
stopCluster(cl)
GetBayesFactor <- function(mll.a, mll.b){
exp(mll.a - mll.b)
}
printMeans <- function(r){
cat(sprintf("パラメーターの数：%d\n", r[[1]]))
cat(sprintf("パラメーターの点推定値（mean）：\n"))
print(theta2param(summary(r[[2]])$statistics[,"Mean"]))
}
len <- length(ptrn)
bf <- matrix(NA, len, len, dimnames=list(c(ptrn), c(ptrn)))
for(i in 1:len){
rownames(bf)[i] <- sprintf("numerator:t=%d", (max(ptrn) + 2 - r[[i]][[1]])/2)
for(j in 1:len){
colnames(bf)[j] <- sprintf("denominator:t=%d", (max(ptrn) + 2 - r[[j]][[1]])/2)
bf[i, j] <- GetBayesFactor(r[[i]][[3]], r[[j]][[3]])
}
}
rownames(bf)[len] <- colnames(bf)[len] <- "no change"
print(bf)
choiced_model_no <- which(bf[,1]==max(bf[,1]))
choiced_model <- r[[choiced_model_no]]
printMeans(choiced_model)
summary(choiced_model[[2]])
plot(choiced_model[[2]])
長い？ — 4つあるモデルを1つの処理で片付けようとしたために、ほんの少し複雑になっていますが、大したことありません。時点ダミー×介入群ダミーを入れた線形混合モデルだと10行以内に終わるとしてでもです。なお、ベイズファクターを用いているので、事前分布は絶対に省略しないでください。ところで記事を書いてから気付いたのですが、推定式に
を含めるのを忘れていますが、まぁ、含めても結果は変わらないと思うので気にしないでください。
4. 推定結果
各時点120とそれなりのサンプルサイズなのですが、ベイズファクターで選択されたモデルの推定結果はt=3で変化となりま･･･と最初書いていたのですが、コードにミスがあったので直したらt=1で変化となりました。真のモデルはt=2で変化のはずなので妙な気も済ますが、t=1とt=2のモデルのベイズ・ファクターを比較すると1.03とほとんど1であり、両者に差が無い事が分かります。t=2のモデルの方が簡素なので、第2時点で変化があったと解釈して良いでしょう。
ベイズファクターで選択されたモデル以外、今回の場合はt=2もチェックした方が良さそうです。今回のコードではprintMeans(r[[2]])で、もしくはMCMCPackを使っているので、summary(r[[2]][[2]])で数字を見ることができます。
なお、set.seed(611)にしておくと、t=2で変化と言う結果になります。例に使うシード値を間違えた感。
