TJO氏にふられた質問で、反復測定分散分析（repeated measures anova）をベイズ推定に置き換えられないかと言う話があったので、質問文から想定される状況を整理し、欠損値補定をFully Bayesian Imputation Method *1で、モデル比較をBayes Factorで行なうコードを描いて試してみました。
結論はRのlme4パッケージあたりで時点ダミー×介入群ダミーを入れた線形混合モデルを推定した方がコード量も計算時間も1/100で済むといった感じなのですが *2、各時点各群の平均と分散をそれぞれ推定できること、欠損値の補定方法のパラメーターも同時推定できる利点もあるので、ファンがうぃーんと言うのが好きな人には良いかもです。うぃーん。
1. 分析する状況の想定
反復測定分散分析が使われているところなので、医療系の実験で、薬剤を与える介入群Tと、偽薬を与える操作群Cの2群を、第1時点、第2時点、第3時点それぞれ調べ、介入群と操作群の違いが出る時点を探すことを考えます。
数式で書くとこんなモデルになります。
は検査指標、
は群の各時点の状態、
は個体ごとの欠損値もある変数、
は介入群と操作群を現し、
は時点、
は個体、
は誤差項を示します。
こう書くと簡素なモデルが1つしかないと思うかも知れませんが*3、と
で評価して、
（T群とG群にずっと差異なし）
（T群とG群は第3時点で差異が出た）
（T群とG群は第2時点から差異）
（T群とG群は第1時点から差異）
の4パターンのモデルを評価することができます。平均だけではなく分散も異なること、平均が同じでも分散が異なればとなるとしましょう。また、どちらの群でも個体の状態が一定でないこと、
を仮定します。
欠損値の入り方は想像がつかなかったのですが、MCARとしておきます。原理的にはMARでも大丈夫ですし、欠損値が入るパターンが分かっていれば、推定モデルを工夫すればMNARにも対応できます。
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, sigma, betaは真の値になるので、本当は分からない 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] } } # xに欠損値をつくる 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) # ある点の確率密度を得る（Parzen Window） 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 } # Chib (1995)に習って、MCMCを繰り返して確率密度を計算する 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)) } # xの係数と欠損値補定用の計4個のパラメーター # t*g*2=12個のパラメーター theta2param <- function(theta){ theta <- theta[-(1:4)] # xと欠損値ダミーの係数のパラメーターは不要 mlen <- length(t)*length(g)*2 # 制約なしの長さ tlen <- length(theta) # 引数の長さ slen <- (mlen - tlen)/2 if(0 < slen){ # 制約がついている # slen=kならば、t=k時点まではCとTは同じ 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) } # theta2param(c(beta, 0, 0, 1, t(mu), t(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){ # i群j時点のscoreを説明するモデルの方の尤度を計算して加算 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){ # thetaの長さが短ければ、周辺尤度の計算過程なので、estimatedの値で補完する 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 # パラメーターの数は16（t=1に変化）, 14（t=2に変化）, 12（t=3に変化）, 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) # Bayes Factorsからモデル選択 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で変化と言う結果になります。例に使うシード値を間違えた感。
*1:"Bayesian methods for missing data Bayes Pharma"と言う文書を参考にしました。
*2:欠損値処理をしなければ、library(lme4); mixed.lmer <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df01); print(summary(mixed.lmer))で終わります。多重代入法による欠損値処理も容易です（ Rで一般化混合線形モデル×多重代入法 - 餡子付゛録゛ ）。
*3:ごく普通に期待値と分散を線形モデルで説明する事も出来ます。
