統計解析の実用において、欠損値補定は21世紀になってから大きく進歩しました。それまでは、とくに統計的仮説検定や区間推定を用いるのに不適切な補定方法が用いられる事も多かったのですが、現在では多重代入法により従来より適切な推定量が計算されるようになっています。

欠損値が少ない場合はlist/pair-wise法を用いて、欠損値が多い場合は多重代入法を用いるのが現在の規範な分析です。Rでの定番パッケージはmiceで、概ねmiceを用いておけば無難となります。ただし、miceの挙動やオプションを理解するのには、単一代入法の知識も必要となります。欠損値削除、単一代入法、多重代入法の概要とRによる実践方法を、順番に確認していきましょう。

1 欠損値の入り方の種類

Rubin (1987)は欠損値の入り方を3種類に分類しました。例えばデータセットに変数$X$と変数$Y$,$Z$,$A$があり、$X$に欠損値があるとします。このとき$X$の欠損値の発生形態は、

$X$と独立に生じているMCAR
$X|Y,Z,A$と独立に生じているMAR
$X|Y,Z,A$と相関して生じているMNAR

に分類できます。統計解析を実践する上では、MCARは欠損値補定が不要、MARは欠損値補定をすべき、MNARは欠損値補定はしたいができないと言う違いがあります。

1.1 データ生成プロセス

あとでコード例を動かすために、欠損値の入り方を変えたデータセットをつくる関数を用意します。

dgp <- function(n = 30, missing.type = "MCAR", var.names = c("x"), x.type = "double"){
    S <- matrix(c(1, -0.3, 0.3, -0.3, 1, 0.3, 0.3, 0.3, 1), 3, 3)
    C <- chol(S)
    V <- t(C) %*% matrix(runif(3*n, 0, 2), 3, n)
    x0 <- x <- drop(V[1, ])
    z <- drop(V[2, ])
    a <- drop(V[3, ])
    if("binominal" == x.type){
        p <- 1/(1 + exp(mean(x) - x))
        x <- x0 <- 1*(runif(n) > p)
    }
    y <- 1 + 2*x - 3*z + rnorm(n)
    if("MAR" == missing.type){
        v <- y + x + z + a
        for(i in var.names){
            v <- v - get(i)
        }
        q <- -3 + 2*v
        p <- 1/(1 + exp(-q))
    } else if("MNAR" == missing.type){
        frml <- as.numeric(x) ~ y + x + z + a
        for(i in var.names){
            frml <- update(frml, as.formula(sprintf(". ~ . - %s", i)))
        }
        r_lm <- lm(frml)
        q <- residuals(r_lm)
        p <- 1/(1 + exp(mean(q) - q))
    } else if("MCAR" == missing.type){
        p <- 0.5
    } else if("NONE" != missing.type) stop("type is incorrect!")
    if("NONE" != missing.type){
        for(i in var.names){
            v <- get(i)
            v[runif(n) > p] <- NA
            assign(i, v)
        }
    }
    if("binominal" == x.type){
        x0 <- as.factor(x0)
        x <- as.factor(x)
    }
    data.frame(y, x, z, a, x0)
}

xとzとaの分散共分散行列Sを用意してデータを生成したあと、欠損値の入り方の指定に応じてxに欠損値を入れています。x0は実データでは観測不能な欠損値が入る前のxで、欠損値補定の精度の比較に用います。

2 欠損値削除

MCARもしくは、欠損値の割合が低いMARの場合は、欠損値補定を行わずに欠損値を削除してしまう方が実践的とされています。

$Y$を$X$と$Z$で説明することを考えているとしましょう。$A$は使いません。$X$と$A$に欠損値があるとします。

Y	X	Z	A
$y_1$	$x_1$	$z_1$	$a_1$
$y_2$	`NA`	$z_2$	$a_2$
$y_3$	$x_3$	$z_3$	`NA`
$y_4$	`NA`	$z_4$	`NA`
$\vdots$	$\vdots$	$\vdots$	$\vdots$

NAはnot availableの略で欠損値です。このとき、$Y$、$X$、$Z$、$A$の全データが揃っている行以外を除去するのがlist-wise法で、分析に使う$Y$、$X$、$Z$が揃っている行以外を除去するのがpair-wise法です。

2.1 データセット

変数xとaに欠損値があるデータセットを作成します。

set.seed(848)
df_mcar <- dgp(30, "MCAR")
# list-wiseとpair-wiseの差をつくるために、変数aにも欠損値を入れる
df_mcar$a[runif(nrow(df_mcar))>0.5] <- NA

2.2 list-wise法

欠損データがある行を除去したサブセットをつくるだけです。

df_listwise <- subset(df_mcar, complete.cases(df_mcar))
nrow(df_listwise)

[1] 8

2.3 pair-wise法

lmコマンドあたりは、欠損データがある行を除いて推定してくれますが、明示的にサブセットをつくれます。

df_pairwise <- subset(df_mcar, complete.cases(df_mcar[, c("y", "x", "z")]))
nrow(df_pairwise)

[1] 17

自明ですが、pair-wise法の方が推定に使えるデータが大きくなります。

3 単一代入法

現在では推奨されなくなっていますが、点推定だけが欲しい場合は手軽に計算できる単一代入法は有用です。多重代入法でも内部的に用いており、多重代入法のライブラリの機能を使うのにも必要な知識になるので、概要は把握しておくだけでも有益です。

3.1 データセット

df_mar <- dgp(100, "MAR")

3.2 平均代入法

欠損していないデータの平均値を、欠損値に入れるだけです。昔はよく使われていたのですが、今はむしろ推定結果が悪くなるので止めるように言われています。使ってはいけません。

df_imp_mean <- df_mar
df_imp_mean$x[is.na(df_imp_mean$x)] <- mean(df_mar$x, na.rm = TRUE)

属性ごとに平均をとり欠損値に代入すると、層化平均代入法になります。

3.3 回帰代入法

非欠損データで予測モデルを作成し、欠損データの行を予測値で埋めます。

r_lm_x <- lm(x ~ z + a, df_mar) # pair-wiseでxを推定 
x_imp <- predict(r_lm_x, df_mar) # 予測値を補定値とする
x_imp[!is.na(df_mar$x)] <- df_mar$x[!is.na(df_mar$x)] # 非欠損データは元の値を補定値に入れる

補定モデルの説明変数は、従属変数を入れても、回帰分析に用いない値を入れても良いです。補定モデルはロジスティック回帰のような非線形モデルを使う事もできます。glmで推定して、type = "response"をつけてpredictした値を補定値として使えます。クラスターを表すダミー変数だけを説明変数とすると、層化平均代入法と同様になります。

3.3.1 確率的回帰代入法

補定モデルの確率的な要素を考慮し、乱数を加味して補定値を定める方法です。そのまま使われる事はまずなく、多重代入法の一部分として使われます。

補定モデルが一般線形回帰の場合は、係数の推定量の標準誤差と誤差項の分散を考慮します。

r_lm_x <- lm(x ~ z + a, df_mar) # 補定モデルを作成する
S <- vcov(r_lm_x) # 分散共分散行列
C <- chol(S) # S == t(C) %*% C になる
b <- coef(r_lm_x) + t(C) %*% matrix(rnorm(ncol(C)), ncol(C), 1) # 乱数で係数を決定
s2 <- deviance(r_lm_x)/df.residual(r_lm_x) # 誤差校の分散
za <- model.matrix(~ z + a, df_mar) # 補定に使う説明変数
x_imp <- za %*% b + rnorm(nrow(za), sd = sqrt(s2)) # 非欠損値を含めた全行を補定する
x_imp[!is.na(df_mar$x)] <- df_mar$x[!is.na(df_mar$x)] # 非欠損データは元の値を補定値に入れる

二項モデルなど誤差項がなく予測値が確率になる場合は、予測値の確率から乱数で補定値を定めます。

3.4 最尤法

最尤法で欠損値を推定し、補定する方法です。回帰代入法では1変数しか扱えませんが、最尤法では補定モデルを多変量分布にすれば、複数の欠損が同時に生じた場合も補定できます。

多変量正規分布を仮定して、最尤法で欠損値補定を行ってみましょう。

# 制約なしなので、変数の平均値とそこから計算した分散共分散行列をパラメーターとする
mu <- apply(df_mar, 2, mean, na.rm = TRUE)[1:4]
frml <- ~ 0 + y + x + z + a
X <- model.matrix(frml, df_mar)
X0 <- apply(X, 1, \(x) x - mu)
Sigma <- X0 %*% t(X0) / (ncol(X0) - 1)
invSigma <- chol2inv(chol(Sigma))
# 欠損値ありの行列をつくる
df_a <- model.frame(frml, data = df_mar, na.action = na.pass)
X_a <- model.matrix(frml, df_a)
# 欠損値ありだけの行列をつくる
X_m <- X_a[!complete.cases(X_a), ]
init.p <- rep(0, sum(is.na(X_a)))
library(mvtnorm)
r_optim <- optim(init.p, function(p){
    X_a[is.na(X_a)] <- p
    sum(dmvnorm(X_a, mu, Sigma, log = TRUE))
}, function(p){ # グラディエントは無くても動く
    na <- is.na(X_m)
    index <- 1:length(X_m)
    X_m[na] <- p
    Y <- X_m - rep(mu, each = nrow(X_m))
    g <- t(-invSigma %*% t(Y))
    g[index[na]]
}, method = "BFGS", control = list(fnscale = -1))
X_a[is.na(X_a)] <- r_optim$par

多変量正規分布に従わない変数の場合も、補定した値を用いて推定を行う場合は、これで問題が無い事も多いです。

補定する値を工夫することもできます。ポアソン分布のパラメーターの平方根が正規分布に従うと仮定し、パラメーターの平方根が多変量正規分布に従うとするように組み合わせると、非負の離散分布が扱えます。二値変数や比率のデータは、アフィン変換で$[0, 1]$の値を$[\epsilon, 1-\epsilon]$に変換したあと、逆ロジット変換で上限下限の無い連続値に変換することで、多変量正規分布に従うと見なすこともできます。スケーリングや平行移動で工夫すれば、区間の異なるデータも取り扱えます。

より厳密に行いたい場合は、多変量正規分布ではなくその分布の尤度関数のパラメーターを推定し、多変量正規分布のときと同様に尤度が最大になるように補定値を求めます。ひとつのレコードに複数の欠損値が入る場合は、

3.5 最尤法による回帰分析

回帰分析で扱う場合は、二種類のアプローチがあります。

3.5.1 完全情報最尤法（FIML）

最尤法が複数のモデルを同時推定できることを活かして、完全データ用のモデルと、欠損データ用のモデルを同時推定します。2つのモデルで共通するパラメーターなどは任意に設定することができます。また、欠損データのモデルは何種類あってもよいので、複数の変数に欠損値が含まれる場合にも拡張できます。

上の説明変数$x$に欠損値が入る例を推定するときは、

\[\begin{align} y_i &= \beta_0 + \beta_1 x_i + \beta_2 z_i + \epsilon_i \quad \mbox{（完全データ用モデル）} \\ y_i &= \gamma_0 + \gamma_1 a_i + \beta_2 z_i + \epsilon_i \quad \mbox{（欠損データ用モデル）} \end{align}\]

と言う2つのモデルを同時推定することになります。

3.5.2 完全ベイズ補定法（FBIM）

欠損値補定モデルと、回帰分析モデルを同時推定する方法です。欠損値補定モデルの予測値で欠損値を補定し、それを回帰分析モデルの尤度計算に用います。欠損値補定モデルの尤度と回帰分析モデルの尤度の合計が最大になるように、2つのモデルのパラメーターを決定します。分類すれば単一代入法ですが、多重代入法に近い標準誤差になり、過小推定バイアスは解消されています。最尤法でも計算できますが、主観的事前分布を目的関数に追加し、マルコフ連鎖モンテカルロ法（MCMC）でも計算できます。

3.6 予測平均マッチング（PMM）

非欠損データで予測モデルを作成し、欠損データの予測値と近い予測値を持つ完全データを1つもしくは複数選び、選んだ完全データの中から無作為に1つ観測値を選び、欠損データを補定する値とします。候補になる完全データは、予測値と予測値の距離が一定以下の場合も、候補数が一定数の場合もあります。候補数が適応的に定まる場合もあります。無作為抽出は予測値と予測値の距離に基づくウェイト付きの場合もあります。

r_lm_x <- lm(x ~ z + a, df_mar) # pair-wiseでxを推定 
x_p <- predict(r_lm_x, df_mar) # 予測値を計算
x_imp <- with(df_mar, {
    x_i <- x
    for(i in seq_len(length(x))){
        if(is.na(x[i])){
            # 予測値が近い3つのうち無作為に1つを選んで補定
            d <- sqrt((x_p[i] - x_p)^2)
            c_x <- x[order(d)]
            n <- 3
            c_x <- c_x[!is.na(c_x)][1:n]
            j <- round(runif(1, 0.5, 0.5 + n)) # 1からnになる
            x_i[i] <- c_x[j]
        }
    }
    x_i
})

回帰代入法や完全ベイズ補定法の場合はモデル定式化に失敗して推定がおかしくなることがあるのですが、予測平均マッチングはその可能性を低くすることができます。一方、欠損データの予測値に近い予測値の観測データが無い場合、とくに外挿になる場合はうまく補定できません。

4 多重代入法

回帰代入法や予測平均マッチングは良い点推定量が得られる一方で、標準誤差を過小評価します。統計的仮説検定を行うと偽陽性が出やすくなり、信頼区間が狭くなると言うことです。補定モデルの標準誤差と誤差項を考慮しないためです。そのため、近年は多重代入法を用いることが一般化しています。

多重代入法では、補定モデルの誤差項と係数の推定量の標準誤差を加味した予測値を用いる単一代入法を何度も繰り返し、多数のサンプルをつくります。そして、それぞれのサンプルごとに推定を行い、すべての推定結果を統合することで、標準誤差の過小推定バイアスを解消します¹。

4.1 Amelia

欠損値補定モデルを多変量正規分布とし、そのパラメーターを、パラメーター⁽ⁿ⁾を用いて補定したデータセットから推定したパラメーター⁽ⁿ⁺¹⁾がパラメーター⁽ⁿ⁾と一致するように、EMBアルゴリズムを用いて算出する方法です。欠損値補定モデルのパラメーターの推定後、多数のサンプルをつくってそれぞれ母集団や回帰モデルのパラメーターの推定を行い、推定結果を統合します。Rではameliaパッケージを用いて計算することができます。

4.2 mice

実践ではmiceと呼ばれる手法（と同名のパッケージ）がよく用いられます。Ameliaと似ていますが、より柔軟な構造になっています。

miceでは、以下のような手順で多数のデータセットを作成します。

すべての変数の欠損値に何らかの単一代入法による補定がされたデータセット$X_0$を作成します。
$k=0$と置きます。
データセット$X_k$を用いて、補定モデル（のパラメーター）を推定します。補定モデルは線形回帰や予測平均マッチングなどを、変数ごとに選択できます。
補定モデルから、確率的誤差を加味した単一代入法による新たな補定値を含むデータセット$X_{k+1}$を作成します。補定されたデータセットから推定するので、複数の変数に欠損があってもそれぞれ推定できます。
$X_k$が収束したと見做せれば、その最初の$k$を$s$として保存します
$k-s+1$が十分な大きさになり、必要な数のサンプルが集まったのであれば終了します。
$k$を$k+1$に置き直します（e.g. $k$が$0$であれば$1$、$1$であれば$2$、$2$であれば$3$）
(3)に戻ります。

これで多数のデータセット$X_0, X_1, X_2, \cdots, X_s, \cdots, X_k$ができます。そして、それらの中で補定モデルのパラメーターが収束したと見做せた後のデータセット（e.g. $X_{s}, X_{s+1}, \cdots, X_k$）を用いて、回帰モデルなどの推定量を出します

miceは数理的な正当化は完全に行われていないようですが、シミュレーションでは良好な結果を示すことが知られています。

4.3 `mice`の利用例

DGPから複数の列に欠損値があるデータセットを作成し、miceで推定しましょう。

set.seed(1251)
df01 <- dgp(100, missing.type = "MAR", var.names=c("x", "z"), x.type = "binominal")
summary(df01)

FALSE        y             x            z                  a          x0    
FALSE  Min.   :-4.9619   0  :24   Min.   :-0.54594   Min.   :0.4522   0:49  
FALSE  1st Qu.:-1.4052   1  :32   1st Qu.:-0.06202   1st Qu.:1.1620   1:51  
FALSE  Median : 0.4248   NAs:44   Median : 0.22697   Median :1.6782         
FALSE  Mean   : 0.2390            Mean   : 0.25125   Mean   :1.6308         
FALSE  3rd Qu.: 1.9452            3rd Qu.: 0.56785   3rd Qu.:2.1565         
FALSE  Max.   : 4.5888            Max.   : 1.33155   Max.   :2.7769         
FALSE                             NAs    :48

欠損値があるのは変数xと変数zで、xは2値の要素型になっています。二値や多値の列はfactor、順序の場合はさらにorderedにするのを忘れないようにしてください。

推定の前に、引数を調整します。しなくてもそれなり動きますが、今回は二値データであることを明示します。

library(mice)
# ループ回数ゼロで実行し、雛形をつくる
template <- mice(df01, maxit = 0)
# 変数xを二項モデルとして補定する
method <- template$method
method["x"] <- "logreg"
method

       y        x        z        a       x0 
      "" "logreg"    "pmm"       ""       ""

pmmは予測平均マッチングの指定です。様々な手法が指定できるのですが、

パラメーター	メソッド	型
`norm`	線形回帰モデル	`numeric`
`pmm`	予測平均マッチング	`numeric`/`factor`/`ordered`
`midastouch`	加重予測平均マッチング	`numeric`/`factor`/`ordered`
`logreg`	二項ロジットモデル	`factor`
`polyreg`	多項ロジットモデル	`factor`
`polr`	順序ロジットモデル	`ordered`

あたりが説明しやすく実用的だと思います。

欠損値の予測に、どの列を用いるかの指定を行列で行います。今回のデータフレームのx0は応用では絶対に分からない真のデータなので、予測には使いません。

# 変数xの補定に使う変数をzとaに限定
predictorMatrix <- template$predictorMatrix
# 変数x0は補定に用いない
predictorMatrix["x0", ] <- predictorMatrix[, "x0"] <- 0
# 変数xはy, z, aで補定
predictorMatrix["x", ] <- predictorMatrix[, "x"] <- c(1, 0, 1, 1, 0)
predictorMatrix

   y x z a x0
y  0 1 1 1  0
x  1 0 1 1  0
z  1 1 0 1  0
a  1 1 1 0  0
x0 0 0 0 0  0

引数を整理したので、実際に多重代入法を行い、推定してみましょう。

# 欠損値補定を行った複数のサンプルを作成
# m: 生成する補定されたサンプルの数
imputed_df01 <- mice(df01, method = method, predictorMatrix, m = 10, maxit = 20, seed = 123, printFlag = FALSE)

# サンプルごとに回帰分析を行う
# withはビルトインではなくmiceパッケージ内の関数
# 変数`a`は補定のための補助変数とし、推定では用いない
imputed_fits <- with(imputed_df01, lm(y ~ x + z))

# 複数の回帰分析の結果を一つに統合する
r_mice <- pool(imputed_fits)

# 推定結果
summary(r_mice)

         term   estimate std.error statistic       df      p.value
1 (Intercept)  0.5334467 0.2493631  2.139237 49.30775 3.739276e-02
2          x1  2.2384396 0.4553241  4.916146 15.37459 1.736873e-04
3           z -2.7210761 0.3551051 -7.662735 46.21744 9.059145e-10

簡単にできましたね。自由度（df）が、完全データにOLSをかけた場合と異なり説明変数ごとにまちまちの数字ですが、これはBarnard and Rubin (1999)の方法に沿って計算されています。

係数を補定しない場合と比較してみましょう。

coef(summary(lm(y ~ x + z, df01)))

             Estimate Std. Error   t value     Pr(>|t|)
(Intercept)  1.569873  0.2055483  7.637489 1.806028e-09
x1           1.469840  0.2866327  5.127956 7.019988e-06
z           -2.082863  0.3137485 -6.638639 4.789631e-08

DGPから真の係数は$1, 2, -3$になるのですが、補定した場合の方が近い結果となっています。一方、標準誤差は補定した場合の方が大きく、偽陰性が出づらい事がわかります。

4.3.1 分散共分散行列

miceは標準誤差（もしくは分散共分散行列の対角成分）に関しては複数推定の結果を統合してくれるのですが、バージョン3から共分散については処理しなくなりました²。Dong and Peng (2013)の説明に準じて計算します。

# within imputation variance
vw <- Reduce("+", lapply(imputed_fits$analyses, vcov)) / r_mice$m
# between imputation variance
qbar <- getqbar(r_mice)
qhats <- sapply(imputed_fits$analyses, coef)
vb <- (1 / (r_mice$m - 1)) * (qhats - qbar) %*% t(qhats - qbar)
# weighted sum of two variances
r_mice_vcov <- vw + (1 + 1 / r_mice$m) * vb

miceが計算する分散共分散行列の対角成分と合致するかチェックしておきましょう。

all(diag(r_mice_vcov) == r_mice$pooled$t)

[1] TRUE

点推定量がr_mice$pooled$estimateに、自由度がr_mice$pooled$dfに入り、分散共分散行列をr_mice_vcovが入りました。

あとはBarnard-Rubinの修正自由度を求めれば、複数の変数を結合した場合の信頼区間の計算と統計的仮説検定を、OLSと同様に行う事ができます。

Barnard-Rubinの修正自由度は変数の値によって変わってくるので、計算がちょっと込み入っています。とは言え、定形パターンなので難しくはありません。例えば、帰無仮説が$x$の係数と$z$の係数の合計が$1$のときのp値を求める場合は、以下のようになります。

mask <- matrix(c(0, 1, 1), 1) # x=1, z=1とおく
se <- sqrt(mask %*% r_mice_vcov %*% t(mask)) # 標準誤差
m <- sum(mask %*% qbar) # \beta_1 * x + \beta_2 * z = \beta_1 + \beta_2
# 自由度の計算（Barnard-Rubinの修正自由度）
vw <- mean(sapply(imputed_fits$analyses, \(r){
        mask %*% vcov(r) %*% t(mask)
    }))
tau <- mask %*% ((1 + 1 / r_mice$m) * vb) %*% t(mask) / vw
df <- (r_mice$m - 1)*(1 + 1 / tau)^2
# p値の計算
pt((m - 1)/se, df)

            [,1]
[1,] 0.002587575

調整自由度計算のために分散共分散行列を計算しなおしているのが、ポイントです。

4.3.2 従属変数の信頼区間と予測区間

lmを用いた線形回帰の場合、mice::predict_mi関数を用いると簡単に出ます。

# 予測区間を出す点
df02 <- data.frame(x = as.factor(c(0, 1, 0, 1)), z = c(-1, -1, 1, 1), a = NA)
predict_mi(imputed_fits, df02, interval = "confidence")

          fit        lwr        upr
1  3.25452279  2.2548667  4.2541789
2  5.49296235  4.1245754  6.8613493
3 -2.18762933 -2.8783156 -1.4969431
4  0.05081024 -0.8091023  0.9107228

predict_mi(imputed_fits, df02, interval = "prediction")

          fit        lwr       upr
1  3.25452279  0.3312326 6.1778129
2  5.49296235  2.4355910 8.5503337
3 -2.18762933 -5.0180872 0.6428286
4  0.05081024 -2.8183212 2.9199417

ソースコード中のコメントを見ると、将来的にはlm以外もサポートする方針のようですが、現在はlmしかサポートしません。

非線形回帰などを用いるときの練習のために、predict_miを使わないで計算をしてみましょう。

まずは信頼区間ですが、p値と同様に計算できます。

evars <- model.matrix(~ x + z, df02)
se <- sqrt(apply(evars, 1, \(e){
    t(e) %*% r_mice_vcov %*% e
})) # 標準誤差
# 自由度の計算（Barnard-Rubinの修正自由度）
df <- apply(evars, 1, \(e){
    vw <- mean(sapply(imputed_fits$analyses, \(r){
            t(e) %*% vcov(r) %*% e
    }))
    tau <- t(e) %*% ((1 + 1 / r_mice$m) * vb) %*% e / vw
    (r_mice$m - 1)*(1 + 1 / tau)^2
})
# 95%信頼区間
a <- 1 - 0.95
fit <- evars %*% qbar
margin <- qt(1 - a/2, df) * se
matrix(c(fit, fit - margin, fit + margin), ncol = 3, dimnames = list(NULL, c("fit", "lwr", "upr")))

             fit        lwr        upr
[1,]  3.25452279  2.2548667  4.2541789
[2,]  5.49296235  4.1245754  6.8613493
[3,] -2.18762933 -2.8783156 -1.4969431
[4,]  0.05081024 -0.8091023  0.9107228

predict_miを使ったときと同じ値になりました。

予測区間は、分散共分散行列から計算される分散に、さらに誤差項の分散を足したものになります。

evars <- model.matrix(~ x + z, df02)
# error term variance
sigma2 <- mean(sapply(imputed_fits$analyses, \(r_lm){
    deviance(r_lm) / df.residual(r_lm)
}))
var <- apply(evars, 1, \(e){
    e %*% r_mice_vcov %*% e
}) # 分散
se <- sqrt(var) # 標準誤差
# 自由度の計算（Barnard-Rubinの修正自由度）
df <- apply(evars, 1, \(e){
    vw <- mean(sapply(imputed_fits$analyses, \(r){
            t(e) %*% vcov(r) %*% e
    }))
    tau <- t(e) %*% ((1 + 1 / r_mice$m) * vb) %*% e / (vw + sigma2) # 信頼区間の計算と異なる
    (r_mice$m - 1)*(1 + 1 / tau)^2
})
# 95%予測区間
m <- evars %*% qbar
a <- 1 - 0.95
margin <- qt(1 - a/2, df) * sqrt(var + sigma2)
matrix(c(m, m - margin, m + margin), ncol = 3, dimnames = list(NULL, c("fit", "lwr", "upr")))

             fit        lwr       upr
[1,]  3.25452279  0.3312326 6.1778129
[2,]  5.49296235  2.4355910 8.5503337
[3,] -2.18762933 -5.0180872 0.6428286
[4,]  0.05081024 -2.8183212 2.9199417

4.4 `mice`がサポートしていない関数の分析結果の統合

mice::poolが対応していない結果オブジェクトの場合は、パッケージmitools::MIextractと組み合わせて、補定したデータセットをそれぞれ推定し、それぞれ係数と分散共分散行列を抜き出し、それぞれ自由度を計算した結果をmiceadds::pool_miを用いて統合する必要があります。

以下でlme4パッケージを使った例を示します。

4.4.1 データセット

ここでは以前作成した別のDGPを用います。

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))

df03 <- 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
        df03$score[s:e] <- rnorm(n/length(g), mean=mu[i, j], sd=sigma[i, j]) + beta*df03$x[s:e]
    }
}

# xに欠損値をつくる
df00 <- df03
df03$x[df03$id %in% sample(1:n, n/2)] <- NA

4.4.2 補定と推定

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

# miceで補定したm個のデータセットをつくる
# 欠損値があるのは変数xのみ
mice.out <- mice(df03, 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)ができるので、後はまとめる

4.4.3 結果の比較

細かく見ていくと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 = df03)

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")
cmp_beta

                              No Missing   listwise    MI/mice
(Intercept)                  -0.03775438 0.07265223 0.02566604
x                             0.96478478 0.90188105 0.85281808
factor(time)2                 0.18491648 0.01659969 0.18434017
factor(time)3                 0.08955894 0.03421298 0.10024424
factor(group)T                0.18728565 0.18978138 0.22427210
factor(time)2:factor(group)T  0.62094379 0.85077402 0.61591506
factor(time)3:factor(group)T  1.88640199 1.94195526 1.84862789

5 まとめ

説明用に色々なコードを示しましたが、lmやglmなど典型的な分析におけるmiceの実践に絞ると必要なコードはさほど無いです。半分ぐらい欠損値があるようなデータセットに遭遇したら、臆せず多重代入法を用いていきましょう。

詳しい説明は高橋・伊藤 (2014)「様々な多重代入法アルゴリズムの比較～大規模経済系データを用いた分析～」などを参照してください。↩︎
バージョン2.xではr_mice$tに分散共分散行列が入りました。多様な結果オブジェクトの操作をbroomパッケージで用いて処理するようになり、計算されなくなったそうです。↩︎

Y	X	Z	A
\(y_1\)	\(x_1\)	\(z_1\)	\(a_1\)
\(y_2\)	`NA`	\(z_2\)	\(a_2\)
\(y_3\)	\(x_3\)	\(z_3\)	`NA`
\(y_4\)	`NA`	\(z_4\)	`NA`
\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)

統計学の鉛筆なめなめ欠損値補定