線形混合モデル欠損値の多重代入補完: multiple imputation
井口豊（生物科学研究所，長野県岡谷市）
最終更新：2025 年 10 月 23 日
1. はじめに
通常の線形回帰モデル（単回帰または重回帰）でも，混合効果モデルでも，実際のデータには，いわゆる欠損値が生じる場合が，しばしばある。一般的には，欠損値を含むデータ行全体を除外してしまう。それに対して，ここでは，線形混合モデル（Linear Mixed Model）を取り上げ，統計ソフト R による多重代入法（multiple imputation）による欠損補完を取り上げる。
2. R による混合モデルで欠損値の取り扱い
線形混合モデルの計算として，ここでは R の lme4 パッケージにある lmer 関数を利用するが， p 値を計算するために，実際には lmerTest パッケージを利用した。
通常の線形モデル関数 lm 同様に， lmer 関数でも，欠損値があると，そのデータセット（データフレーム行）は，丸ごと削除になる。いわゆる「リストワイズ削除」である（Wikipedia, Listwise deletion 参照）ただし警告もなく自動的に削除されるので，注意が必要だ。
例えば，以下のような 12 行のデータセットで，そのうち 4 行に欠損値 NA を含むデータを想定する。欠損値を含んだままで lmer 関数に投入した場合と，欠損値を含む 4 行のデータセットを予め除外して lmer 関数に投入した計算を比較してみる。
###################################
# データ例
# 12 行データ中， 4 行に欠損値あり
Input<- ("
id x y
1 2 17
2 4 NA
3 NA 50
4 8 65
1 2 16
2 4 33
3 6 56
4 8 69
1 NA 20
2 4 41
3 6 NA
4 8 63
")
dat.1<- read.table(
textConnection(Input), header=T
)
library(lmerTest)
# 欠損値を含んだまま計算した線形混合モデル
mod<- lmer(y ~ x + (1 | id), data = dat.1)
summary(mod)
# 欠損値を予め除外して計算した線形混合モデル
dat.2<- na.omit(dat.1)
mod<- lmer(y ~ x + (1 | id), data = dat.2)
summary(mod)
#######################
結果を抜粋すると，以下のとおりである。
# 欠損値を含んだまま計算した線形混合モデル
Number of obs: 8, groups: id, 4
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.3706 4.5156 1.4973 0.525 0.6668
x 8.1934 0.8097 1.3540 10.120 0.0306 *
# 欠損値を予め除外して計算した線形混合モデル
Number of obs: 8, groups: id, 4
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.3706 4.5156 1.4973 0.525 0.6668
x 8.1934 0.8097 1.3540 10.120 0.0306 *
どちらも，計算に使われたデータ行数は，
Number of obs: 8
であり，全く同じ結果になることがわかる。つまり，欠損値があるデータ行は，その全体が削除されている。
3. 多重代入による欠損補完
実際に線形混合モデルを適用した多重代入による欠損補完のデータセットとして， R の nlme パッケージにあるデータセット Orthodont を利用する。
もとのデータにおける変数 age は 8, 10, 12, 14 という定期的な整数値だが， ここでは age に区間 (0, 1) 一様乱数を加えて， 不定期な連続変数とした。
さらに， age と Sex 変数の交互作用項を投入するので，変数の中心化（centering）を行う。これを忘れがちだが，最近は推奨される変数変換であり， R の多重代入法パッケージ mice の解説である Van Buuren and Groothuis-Oudshoorn (2011), p.29 でも次のように書かれている。
Interactions between two continuous variables are often defined by subtracting the mean and taking the product.
Van Buuren, S. and Groothuis-Oudshoorn, K. (2011) mice: Multivariate imputation by chained equations in R. Journal of statistical software, 45: 1-67.
二値のカテゴリ変数に対しても，しばしば使われる 01 ダミー変数ではなく， +0.5, -0.5 などのコード変換した変数にしたほうが良い。コーディングに関しては，例えば，以下の文献参照。
Alday, P. M. (n.d.) Coding schemes for categorical variables
今回は，目的変数 distance と説明変数 age に欠損値 NA を代入する。
###################################
library(lmerTest)
data(Orthodont,package="nlme")
dat<- Orthodont
set.seed(4222)
# age を不定期な連続変数化
# age 中心化と Sex コード変換
dat$age<- as.numeric(scale(
dat$age + runif(length(dat$age), 0, 1),
center = TRUE, scale = FALSE
))
dat$Sex<- ifelse(
dat$Sex == "Male",-0.5, +0.5
)
## 欠損なしの線形混合モデル
mod.0<- lmer(
distance ~ age*Sex + (1 | Subject),
data = dat
)
summary(mod.0)
## 欠損代入
dat.1<- dat
md <- as.logical(rbinom(nrow(dat.1), 1, 0.2))
ma <- as.logical(rbinom(nrow(dat.1), 1, 0.2))
dat.1[md, "distance"] <- NA
dat.1[ma, "age"] <- NA
## 欠損除外の線形混合モデル
# 欠損値を含むデータ行の除外
dat.2<- na.omit(dat.1)
dat.2$age<- dat.2$age - mean(dat.2$age)
mean(dat.2$age)
mod.2 <- lmer(
distance ~ age*Sex + (1 | Subject), dat.2
)
summary(mod.2)
#######################
結果は以下のとおりである。
まず，欠損値が無い場合の結果として，固定効果を示す。
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 23.80789 0.38110 24.99881 62.472 < 2e-16 ***
age 0.61685 0.06129 79.08897 10.064 8.21e-16 ***
Sex -2.30219 0.76220 24.99881 -3.020 0.00575 **
age:Sex -0.27226 0.12259 79.08897 -2.221 0.02921 *
次に，欠損値がある場合の結果，つまり，欠損行が削除された結果として，固定効果を示す。
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 23.77247 0.41250 23.96447 57.630 < 2e-16 ***
age 0.58764 0.07487 41.76462 7.849 9.45e-10 ***
Sex -2.33673 0.82501 23.96447 -2.832 0.00922 **
age:Sex -0.24321 0.14974 41.76462 -1.624 0.11185
欠損値がある状態で計算すると，もとの欠損値が無い場合と比べて，交互作用が有意でなくなっていることに気づく。
3. 多重代入による欠損補完
多重代入法では，補完を繰り返して生成された複数のデータセットの基づく計算をプールして使うことが多い。その場合，以前は， coef 関数と vcov 関数が使われたが，最近は，その方法は推奨されなくなった。例えば， pool 関数のヘルプ（あるいは ?pool と打ち込む）を見ると，以下のような注意書きがある。
In versions prior to mice 3.0 pooling required that coef() and vcov() methods were available for fitted objects. This feature is no longer supported. The reason is that vcov() methods are inconsistent across packages, leading to buggy behaviour of the pool() function.
多重代入法を用いて，今回のような線形混合モデルの欠損補完するためには， mice パッケージにある mice.impute.panImpute 関数を利用するか，または，そのオリジナルな関数である mitml パッケージの panImpute 関数を利用すると便利である。
#############
# 多重代入補完
library(mitml)
dat.1$age<- dat.1$age - mean(dat.1$age, na.rm =T)
mean(dat.1$age, na.rm =T)
dat.1$age.Sex <- dat.1$age * dat.1$Sex
## 欠損推定の式
# ～左辺が欠損あり変数，～右辺が欠損なし変数
fm<- distance + age + age.Sex ~ Sex + (1 | Subject)
imp <- panImpute(
dat.1,
formula = fm,
n.iter = 20000,
m = 500,
silent = T
)
implist <- mitmlComplete(imp)
# fit multilevel model using lme4
fit.lmer <- with(
implist,
lmer(distance ~ age + Sex + age.Sex + (1 | Subject))
)
testEstimates(fit.lmer)
################
結果の主要部を抜粋すると，以下の通り。
Final parameter estimates and inferences obtained from 500 imputed data sets.
Estimate Std.Error t.value df P(>|t|)
(Intercept) 23.772 0.407 58.447 1.593e+05 0.000
age 0.601 0.064 9.411 2.120e+03 0.000
Sex -2.204 0.812 -2.715 1.815e+05 0.007
age.Sex -0.281 0.132 -2.120 1.853e+03 0.034
欠損補完することによって，もとの欠損なしの結果と同様に，交互作用が有意になることがわかる。
いつでも，このように上手くいくわけではないが，欠損がある結果と，それを補完した結果を比較検討することは，感度分析（sensitivity analysis）の観点からも重要だろう。