[R][初心者の質問] mice で多重代入法の結果を統合できない時の対処法

概要

結論から言うとヘルプドキュメントを読めばどうすべきか分かる問題であるが, 初心者が初心者である所以はそこが分からないことなので (誰でも最初は初心者だったのでそれが悪だとは言わない) 書いておくことにもある程度価値はあるだろう.

R-wakalang にて最近 miceパッケージを使った多重代入法 (multiple imputation) がうまくいかないという質問を何度か見かけた. mice を使った多重代入法について日本語で書かれた教科書として高橋 and 渡辺 (2017)によるものがあるが, 今回のようなケースには詳しくないためここに解決法を書いておく (パッケージの使い方はヘルプドキュメントに書いてあるためこの教科書に問題があるわけではない. ちゃんとこの教科書などを読んで多重代入法の使い方を理解してから使ってほしい).

注意: この記事のパッケージ解説は公式ドキュメントではありません. パッケージは常に改良されているため, なるべく最新の公式ドキュメントを参考にしてください.

概要
問題
解説と解決方法
まとめ
補足: Rubin のルール
参考文献

問題

類似する質問を投稿した人間が複数いるが, 共通するのは多重代入法で補完したデータにモデルや検定を当てはめて, pool 関数で結果を統合しようとしたらエラーが出たというもの. 例えばこんなコードを実行したとする.

require(mice)
imp <- mice(nhanes)
fit <- with(imp, t.test(bmi ~ hyp))
pool(fit)

するとこのようなエラーが出るはず.

エラー: Problem with `summarise()` input `ubar`.
x Column `std.error` not found in `.data`
i Input `ubar` is `mean(.data$std.error^2)`.
Run `rlang::last_error()` to see where the error occurred.

今回は基本パッケージの t.test を適用した結果だが, 生存時間分析の timereg::aalen を使ってうまくいかなかったという投稿もあった. こちらの場合は以下のようなエラー.

 エラー: No tidy method for objects of class aalen
 追加情報:  警告メッセージ:
  get.dfcom(object, dfcom) で:  Infinite sample size assumed.

解説と解決方法

確かにこれらのエラーメッセージだけでは原因がどこにあるのか少し分かりにくいかもしれない. まずは mice::pool の内部の挙動を解説する必要がある. この関数は「Rubin のルール」に基づいて結果の統合をする関数であるが, 内部では broom::tidy 関数を使ってモデルの推定結果を取り出すという処理がある. この tidy 関数は, 従来テキストで表示されていた回帰分析などの結果を表形式に変換してくれるものである. たとえば lm 関数の結果は従来なら以下のようにして表示できる.

fit2 <- lm(mpg ~ wt + qsec, data = mtcars)
summary(fit2)

Call: lm(formula = mpg ~ wt + qsec, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max
 -4.3962 -2.1431 -0.2129  1.4915  5.7486

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
    (Intercept)  19.7462     5.2521   3.760 0.000765 ***
wt           -5.0480     0.4840 -10.430 2.52e-11 ***
qsec          0.9292     0.2650   3.506 0.001500 **
 ---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.596 on 29 degrees of freedom
Multiple R-squared:  0.8264,    Adjusted R-squared:  0.8144
F-statistic: 69.03 on 2 and 29 DF,  p-value: 9.395e-12

これを tidy 関数で変換するとこうなる.

broom::tidy(fit2)

表 1: `broom::tidy` 関数による出力
term	estimate	std.error	statistic	p.value
(Intercept)	19.746	5.252	3.759	0.000
wt	-5.047	0.483	-10.429	0.000
qsec	0.929	0.265	3.506	0.001

モデルオブジェクトのデータ構造は実装によってバラバラなので, 出力がこの形式で統一されることで pool 関数は様々なモデルに対して結果の統合処理をすることができる. summary 関数の出力をコピーしてレポートや卒論にそのまま貼り付ける人はいないだろうから, tidy 関数は単体でも結果のレポーティングに便利である (もちろん, stargazer とかを使っても良い).

それはさておき, tidy 関数はあらゆる統計モデルの結果を表形式に変換できるわけではない. これは tidy のヘルプドキュメントに書いてあるとおりで, 対応しているモデルの一覧に, Aalen モデル (timereg::aalen) は含まれていない. つまり aalen のエラーは tidy に対応していないモデルを使ったのが原因である*1. 一方で t.test の場合は表形式にすることができる(表 2).

tidy(t.test(bmi ~ hyp, data = nhanes))

表 2: `t.test` を `broom::tidy` 関数で変換した場合
estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method
-0.483	26.441	26.925	-0.309	0.7616	13.961	-3.835	2.868	Welch Two Sample t-test

しかし表 1 と比較すると項目が異なる. エラー内容が Column ‘std.error‘ not found in ‘.data‘ だったことからも, これは標準誤差の列が存在しなかったのがエラーの原因だとわかる.

では tidy 関数や pool 関数が対応していない場合, どうやって結果を統合するか, いくつか方法がある.

pool.scalar を使う
pool 関数が対応している等価な別のモデルに置き換える
自分で処理を書く (そして broom パッケージにコミットする)

たぶん (1 - 2) が簡単で, (3) はそれなりに複雑なモデル, 例えば標準誤差を代数的に計算できない場合とかだろう.

(1) の方法

pool が対応していない場合, 1変数バージョンの pool.scalar 関数を使うのが最も簡単な方法だろう. 名前の通り, 単一の統計量に対して結果の統合を行う関数である. 必要なのは (1) 多重代入ごとの推定値, (2) それらに対応する分散 (標準誤差を直接与えてはいけない), (3)元のデータのサンプルサイズ, (4) モデルの推定するパラメータの数である.

試しに lm に対して結果が一致するか見てみる. lm は pool でも統合可能だが, ならば pool.scalar 関数でもできるはずである. 公式のサンプルコードを参考に, 以下のような回帰モデルの, age のパラメータについてのみ比較してみる.

imp <- mice(nhanes, maxit = 2, m = 2, print = FALSE, seed = 42)
fit3 <- with(data = imp, lm(bmi ~ age))
pool(fit3)[2, ]

Class: mipo    m = NA
   term m  estimate      ubar          b         t dfcom       df        riv     lambda       fmi
2  age 2 -1.875604 0.8647848 0.01895132 0.8932118    23 20.13587 0.03287174 0.03182557 0.1155202

次に, pool.scalar で age の結果を統合する. Rubinのルールによる計算に対応する処理は pool.scalar 関数でのみなされており, それ以降の処理は結果を表形式に整形するためのものである.

pool_age <- pool.scalar(
  Q = unlist(lapply(fit3$analyses, function(x) coef(x)["age"])),
  U = unlist(lapply(fit3$analyses, function(x) summary(x)$coefficients["age", "Std. Error"]))^2,
  n = NROW(nhanes),
  k = 2
  )
pool_age <- append(list(term = "age"), pool_age)
out <- as.data.frame(pool_age)[1, -(3:4)]
colnames(out)[c(3, 8)] <- c("estimate", "riv")
out

fit3$analyses はリストオブジェクトで, その各要素は多重代入ごとに適用したモデルの結果である. そこで上記では lappy 関数を使って必要な値を取りだす処理を並列して行っている. lm オブジェクトならば coef 関数で回帰係数を取り出せる. そして標準誤差は suumary 関数の結果から取りだし, 忘れずに2乗する. この辺は何通りかやり方があり, モデルオブジェクトによって全く方法が変わってくるはずなので, 推定値と分散 (or 標準誤差) をどうやって取りだすかは自分で考えねばならない.

表 3 は pool.scalar 関数による結果. pool 関数と同様に, summary 関数を使えば p 値を出してくれる. 6桁まで一致しているので同じ結果が出ていると見て良いだろう.

表 3: `pool.scalar` による計算結果
term	m	estimate	ubar	b	t	df	riv	fmi
age	2	-1.875604	0.8647848	0.0189513	0.8932118	20.13587	0.0328717	0.1155202

(1) の補足: その他の拡張パッケージについて

mice 以外のパッケージも探してみる. miceadds パッケージは名前の通り mice にない機能を補うパッケージで, mitools パッケージも同様のモチベーションで作られている. mitools::MIcombine 関数と miceadds::pool_mi 関数は mice::pool.scalar の多変量バージョンと言える. 補完データごとのパラメータ推定値をリストオブジェクトにまとめて結果を与えることで計算できる. 以下, 比較に使ったコード. より一般的に, パラメータが多くても動作することを確認するため先ほどよりモデルに使用する変数を増やしてみる.

require(mitools)  # ver.2.4
require(miceadds) # ver. 3.10-28
imp <- mice(nhanes, maxit = 2, m = 10, seed = 42)
models <- with(data = imp, exp = lm(bmi ~ hyp + chl))
summary(pool(models), conf.int = T)
summary(pool_mi(
  qhat = lapply(models$analyses, function(x) summary(x)$coefficients[, "Estimate"]),
  se = lapply(models$analyses, function(x) summary(x)$coefficients[, "Std. Error"]),
  dfcom = 22
))
summary(MIcombine(
  results = lapply(models$analyses, function(x) summary(x)$coefficients[, "Estimate"]),
  variances = lapply(models$analyses, function(x) summary(x)$coefficients[, "Std. Error"]^2),
  df.complete = 22
))

MIcombine 関数はt統計量とp値を表示してくれないようだ. これも推定値と標準誤差, そして検定統計量は6桁まで一致するが, MIcombine のみ区間推定の誤差が大きい. デフォルトではどれも95%信頼区間を表示するので, それ以外の点で違いがあることになる*2.

また, 細かい解説をしておくと, 上記の書き方はやや不正確である. pool.scalar 同様に1つ1つが独立したパラメータなら問題ないが, 実際にはそうとは限らない. 重回帰モデルであっても共分散行列を求めるのが正確である. 今回の例では結果が変わらないが, 以下のように書くのが適切である.

pool_mi(
  ...,
  u = lapply(models$analyses, vcov)
)
MIcombine(
  ...,
  variances = lapply(models$analyses, vconv),   
)

さらに, これら2つの関数は, 「完全データの自由度」がデフォルトで無限大になっているため, 手動で指定しないと pool 関数とp値や区間推定の結果が一致しない*3 *4.

というわけで, 現状では単変量ならpool.scalar関数を使い, 2つ以上のパラメータがあるのならmiceadds::pool_mi関数を使うのが簡単なやり方だろう.

なお, 多重代入法に関連するパッケージには Amelia というのもあるが, これは mice で実装されていない新しい多重代入法を提供するだけで結果の統合機能はない (Amelia のサポートする手法の解説は高橋 and 渡辺 (2017) を参考に).

(2) 別のモデルに置き換える

これは適用できるモデルが限定されるやり方である. 例えば Heymans and Eekhout (2019) の 5章に書いてあるように, t.test 関数のデフォルトは Welch の検定で, 単回帰モデルの係数の有意検定と同じ値になる. よって 1元配置分散分析 (1-way ANOVA) の結果は lm 関数を使って統合できる. それ以外でも変わったモデルのように見えて実質的に線形回帰モデルということはよくあるので, lm 関数で同じことができないか考えてみるべきだろう.

(3) 自分で計算処理を書く

(1), (2) でうまくいかない, あるいは (1) で計算できない, モデル全体の統計量も正確に求めたいなら, Rubin のルールを自分で実装する必要がある. もしやりたいのなら, それは次の節が参考になる.

まとめ

miceパッケージで pool関数を使った場合に起こる冒頭のようなエラーは, miceパッケージがあなたの使用している統計モデルを提供しているパッケージに対応していないのが原因
モデルのパラメータの推定値と標準誤差がわかっているならば, モデルのパラメータが1つならpool.scalar関数を, 2つ以上のパラメータがあるなら miceadds パッケージの pool_mi 関数を使う
そうでないなら, Rubin のルールを理解して自分で計算する

補足: Rubin のルール

Rubin のルールに関しては 1987年の Donald Rubin の教科書が毎度引用されるが, 古いし高いので大学図書館に自由にアクセスできる身分でない私には確認する気がおきない. mice パッケージなども紹介している高橋 and 渡辺 (2017) でも, Rubin のルールについて大まかなアイディアは説明しているが, 細かい手順までは書いていない. なので調べるのが少し面倒なのだが, mice パッケージの開発者の一人である Buuren の教科書 (mice パッケージのリポジトリからもリンクされている) van Buuren (2018) “Flexible Imputation of Missing Data” の2章と5章や Heymans and Eekhout (2019) “Applied Missing Data Analysis With SPSS and RStudio” の9-10章では, 最近の修正された方法も含め, それなりに詳しく書かれている. これらはWeb上で公開されておりこれとヘルプドキュメントを読んでねで済む話で, これをコピペして書く意味もないのだが, まあ今回は「初心者の質問」への回答なので書く.

(これを書いてから気づいたのだが, 高井, 星野, and 野間 (2016) は計算できるソフトウェアの解説はないものの, 4章で Rubin のルールの詳しい話に言及している)

まず, Rubin のルールは多重代入ごとに求めたパラメータ $\{\theta_{i}\}_{i=1}^{m}$ が正規分布するという前提である*5 (よってモデルによっては Rubin ルールで統合するのが不適切な場合もありえる). ここから, 統合後 (pooled) のパラメータは単純にこれらの標本平均になる.

$\begin{align} \bar{\theta}_{\mathit{pool}}:= & \frac{1}{m}\sum_{i=1}^{m}\theta_{i}\end{align}$

各標準誤差 $\{\mathit{se}_{i}\}$ も, 群間 (Between)・郡内 (Within) 分散を考慮して次のように統合した標準誤差を求められる.

$\begin{align} \mathit{se}_{\mathit{pool}}:= & \sqrt{V_{T}},\\ V_{\mathit{T}}:= & V_{W}+V_{B}+\frac{V_{B}}{m},\\ V_{W}:= & \frac{1}{m}\sum_{i=1}^{m}\mathit{se}_{i}^{2},\\ V_{B}:= & \frac{\sum_{i=1}^{m}(\theta_{i}-\bar{\theta}_{\mathit{pool}})^{2}}{m-1}\end{align}$

さらに, $\mathrm{H}_{0}:\theta=\theta_{0}$ の帰無仮説に対して, 以下のような検定統計量が帰無仮説のもとで近似的にt分布に従う.

$\begin{align} t:= & \frac{\bar{\theta}_{\mathit{pool}}-\theta_{0}}{\sqrt{V_{T}}}\sim t(\mathit{df})\end{align}$

信頼区間の計算もここからできるので, あとはt分布の自由度 $\mathit{df}$ の計算だけだが, これが少し複雑である. 当初は Rubin によって以下の $\mathit{df}_{\mathit{old}}$ を丸めた値を自由度としていた.

$\begin{align} \mathit{df}_{\mathit{old}}:= & \frac{m-1}{\lambda^{2}}\\ = & (m-1)\left(1+\frac{1}{r^{2}}\right)\end{align}$

しかしこの方法では, 特にサンプルサイズの小さい時に多重代入法の個別の補完結果に対して過大になる (欠損値が一切ない場合よりも大きくなることすらありうる) 傾向があることから, pool 関数では後に発表された Barnard and Rubin (1999) による修正方法を採用している. よって, この関数が表示する df は以下で計算される.

$\begin{align} \mathit{df}:= & \mathrm{round}\left(\frac{\mathit{df}_{\mathit{old}}\mathit{df}_{\mathit{observed}}}{\mathit{df}_{\mathit{old}}+\mathit{df}_{\mathit{observed}}}\right)\end{align}$

$\lambda$ は Proportion of total variance due to missingness で, $r$ は RIV (relative increase in variance) を表し, 上記 $\mathit{df}$ の丸める前の値とともにそれぞれ pool 関数で lambda と riv と df して出力される. $\lambda$ , $\mathit{RIV}$ , $\mathit{df}_{\mathit{observed}}$ はそれぞれ以下で与えられる.

$\begin{align} \lambda:= & \frac{V_{B}+V_{B}/m}{V_{T}},\\ \mathit{RIV}:= & \frac{V_{B}+V_{V}/m}{V_{W}},\\ \mathit{df}_{\mathit{observed}}:= & \frac{\mathit{df}_{\mathit{com}}+1}{\mathit{df}_{\mathit{com}}+3}\mathit{df}_{\mathit{com}}(1-\lambda),\\ \mathit{df}_{\mathit{com}}:= & n-k\end{align}$

$n,k$ は既に出てきたサンプルサイズとパラメータの数である. つまり $\mathit{df}_{\mathit{com}}$ は, サンプルが仮に欠損のない完全データ (complete data) だった場合の自由度を意味する (これも pool 関数で dfcom として表示される). これで計算ができるようになった.

この定義では, $\mathit{df}_{\mathit{com}}\to\infty$ のとき $\mathit{df}\to\mathit{df}_{\mathit{old}}$ となり, サンプルサイズの増加に従い漸近的に古いルールに一致することが分かる. 一方で, $\lambda=1$ の場合は $\mathit{df}=0$ となる. 自由度がゼロの分布というのは意味不明である. これは $\lambda$ の定義から, 欠損値補完のばらつきが非常に大きい場合に起こりうる. つまり欠損値によって失われた情報が多すぎるため, この補完方法では信頼できる推定ができない可能性があるということになる.

これが多変量の場合も自然に拡張できる. 補完ごとのパラメータベクトルを $\boldsymbol{\theta}_{i}$ , 対応する共分散行列を $\boldsymbol{\Omega}_{i}$ とすると,

$\begin{aligned} \bar{\boldsymbol{\theta}}_{\mathit{pool}}:= & \frac{1}{m}\sum_{i=1}^{m}\boldsymbol{\theta}_{i},\\ \mathit{se}(\bar{\boldsymbol{\theta}}_{\mathit{pool}}):= & \mathrm{trace}(\mathbf{V}_{T}^{1/2}),\\ \mathbf{V}_{T}:= & \mathbf{V}_{W}+\left(1+\frac{1}{m}\right)\mathbf{V}_{B},\\ \mathbf{V}_{w}:= & \frac{1}{m}\sum_{i=1}^{m}\boldsymbol{\Omega}_{i},\\ \mathbf{V}_{B}:= & \frac{1}{m-1}\sum_{i=1}^{m}(\boldsymbol{\Omega}_{i}-\bar{\boldsymbol{\Omega}})(\boldsymbol{\Omega}_{i}-\bar{\boldsymbol{\Omega}})^{\top},\\ \bar{\boldsymbol{\Omega}}:= & \frac{1}{m}\sum_{i=1}^{m}\boldsymbol{\Omega}_{i}\end{aligned}$

以降は $\lambda$ を以下の $\bar{\lambda}$ に置き換えた以外は同じである.

$\begin{align} \bar{\lambda}:= & \left(1+\frac{1}{m}\right)\frac{1}{k}\mathrm{trace}(\mathbf{V}_{B}\mathbf{V}_{T}^{-1})\end{align}$

なお RIV は以下のようになる.

$\begin{align} \bar{r}:= & \left(1+\frac{1}{m}\right)\frac{1}{k}\mathrm{trace}(\mathbf{V}_{B}\mathbf{V}_{W}^{-1})\end{align}$

なお, SPSS で使用できる統合処理はさらにその後提案された別の方法に従っているため, これとも少し違う値になるらしい. また, 正規分布の仮定が適切でないケースの解説は Buuren の5章で言及がある.

参考文献

Barnard, J, and Donald B. Rubin. 1999. “Small-Sample Degrees of Freedom with Multiple Imputation.” Biometrika 86 (4): 948–55. DOI: 10.1093/biomet/86.4.948.

Heymans, Martijn W, and Iris Eekhout. 2019. “Applied Missing Data Analysis with SPSS and RStudio.” https://bookdown.org/mwheymans/bookmi/.

van Buuren, Stef. 2018. Flexible Imputation of Missing Data. Second. Chapman and Hall/CRC. https://stefvanbuuren.name/fimd/.

高井啓二・星野崇宏・野間久史. 2016. 『欠測データの統計科学』. 岩波書店.

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際第1巻)

発売日: 2016/04/20
メディア: 単行本（ソフトカバー）

高橋将宜・渡辺美智子. 2017. 『欠測データ処理 — R による単一代入法と多重代入法—』. 統計学 One Point 5. 共立出版.

欠測データ処理: Rによる単一代入法と多重代入法 (統計学One Point)

作者:将宜, 高橋,美智子, 渡辺
発売日: 2017/12/09
メディア: 単行本

*1:より正確には, tidy または glance 関数に対応していないとエラーになる.

*2:p値が表示されないので, 自由度の計算が違うのではないかと疑っている. しかしヘルプドキュメントが書きかけだったり, どの理論に準拠した実装なのかとかを書いてないのでどういう仕様なのかわからない. 他の方法もあることだしソースコードを確認するまでの手間はかけたくない.

*3:ヘルプによれば pool 関数の場合は, サンプルサイズ - パラメータ数で計算する. この数値は (1) glance 関数で取りだした値を使う, (2) df.residual 関数で取りだした値を使う (3) いずれも不可能な場合, その旨を表示し, 自由度を999999 とする. また, この値は ‘dfcom‘として保存されるし, 引数で任意に指定もできる.

*4:サンプルサイズが大きければこのような近似で問題ないのだが, デフォルトが無限大なのは不親切な気がする. 必須引数にしてくれたほうがミスが少なくなると思うのだが...

*5:仮説検定が不要なら, この仮定は緩和できる

ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません