Rで線形混合モデル×多重代入法

反復測定分散分析をlme4パッケージとmiceパッケージを使った線形混合モデル*1×多重代入法で行なう作業を確認したので公開します*2。以前のエントリー「反復測定分散分析を主観ベイズ推定に置き換えよう」で、こっちの方が楽と書いたのですが、欠損値補定がやや煩雑だったので。

1. データセット

前回とほとんど同じモノにしますが、比較のために欠損値なしのデータセットも作っておきます。細かい説明は以前のエントリーを参照してください。

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に欠損値をつくる
df00 <- df01
df01$x[df01$id %in% sample(1:n, n/2)] <- NA

2. 補定と推定

miceパッケージのお手軽関数は、lme4パッケージが出力するオブジェクトを現時点で処理できません。miceが補定*3したデータセットに逐次推定をまわし、その結果となる係数と分散共分散行列の推定量をリストにまとめて、miceaddsパッケージ、mitoolsパッケージに渡して推定結果を統合します。こう書くと煩雑で大変そうですが、偉い人が情報をまとめておいてくれました*4。

# 利用パッケージを呼び出す
library(lme4)
library(mice)
library(miceadds)
library(mitools)

# miceで補定したm個のデータセットをつくる
# 欠損値があるのは変数xのみ
mice.out <- mice(df01, blocks=c("x"), print=FALSE, m=100)

# 補定したm個のデータセットをそれぞれ推定する
mods <- lapply(1:mice.out$m, function(i){
    # m個の補定データセットの中からi番目を取り出す
    dfm <- complete(mice.out, action=i)
    # i番目で推定をかけて結果を戻すとリストに格納される
    lmer(score ~ x + factor(time)*factor(group) + (1|group), data = dfm)
})

# 推定結果の集合から係数を抜き出す
cmod <- MIextract(mods, fun=function(r_lmer){
    # lme4の癖でcoef(r_lmer)が複数行返ってくるため
    summary(r_lmer)$coefficients[,1]
})

# 推定結果の集合から分散共分散行列を抜き出す
vmod <- MIextract(mods, fun=function(r_lmer){
    # lme4の癖でvcov(r_lmer)が行列にならないため
    as.matrix(vcov(r_lmer))
})

# 自由度を計算
df <- with(summary(mods[[1]]), length(residuals) - length(coefficients[,1]))

# 推定結果を統合
r_pool <- pool_mi(qhat=cmod, u=vmod, dfcom=df)

# coef(r_pool)とvcov(r_pool)ができるので、後はまとめる

欠損値補定後にクロス項をつくっていますが、今回はクロス項に使う変数は補定されていないので問題ないです。unable to evaluate scaled gradientや Model failed to converge: degenerate Hessian with 1 negative eigenvaluesと言う警告が出るのですが、どうも内部で生成している乱数が原因で、推定結果を大きく左右するものでは無さそうなので無視します。

3. 結果の比較

欠損値補定、劇的に推定結果が代わるというものではないのですが、欠損値なしの場合と、listwise法の場合と、miceの場合で比較してみましょう。

# 欠損値なしデータでの推定（比較用）
mixed.lmer.c <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df00)

# listwise法（比較用）
mixed.lmer.m <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df01)

cmp_beta <- cbind(summary(mixed.lmer.c)$coefficients[,1], summary(mixed.lmer.m)$coefficients[,1], coef(r_pool))
cmp_se <- cbind(summary(mixed.lmer.c)$coefficients[,2], summary(mixed.lmer.m)$coefficients[,2], sqrt(diag(vcov(r_pool))))
colnames(cmp_beta) <- colnames(cmp_se) <- c("No Missing", "listwise", "MI/mice")

listwise法よりも欠損値なしに近づけば改善ですが、係数の比較は改善しているか、悪化しているのか良く分かりません。

標準誤差は改善していると言えるでしょう。

A. Ameliaを使う

lme4パッケージとmiceパッケージを併用するにはどうしたらいいかと言う質問に、Ameliaパッケージを使えと回答がついていたのを見かけて試してみたのですが、微妙でした。Amelia向きのデータ生成のはずですが。

library(Amelia)
library(merTools)
library(plyr)

# Ameliaで補定したm個のデータセットをつくる
amelia.out <- amelia(df01, idvars="id", ts="time", cs="group", m=100)

# 補定したm個のデータセットをそれぞれ推定する
mods <- lapply(amelia.out$imputations, function(dfa){
    lmer(score ~ x + factor(time)*factor(group) + (1|group), data=dfa)
})

# 推定結果の集合から係数を抜き出す
beta <- sapply(mods, function(r_lmer){
    summary(r_lmer)$coefficients[, 1]
})

# 推定結果の集合から標準誤差を抜き出す
se <- sapply(mods, function(r_lmer){
    summary(r_lmer)$coefficients[, 2]
})

# 合成する（分散共分散行列でなくて良いのか感があるが、良いらしい）
r_mi <- mi.meld(q=beta, se=se, byrow=FALSE)

# 比較表に追加
cmp_beta <- cbind(cmp_beta, c(r_mi$q.mi))
cmp_se <- cbind(cmp_se, c(r_mi$se.mi))
colnames(cmp_se)[4] <- colnames(cmp_beta)[4] <- "Amelia"