データを差し替えて再投稿．

リストワイズは「すべての質問に回答した人」で条件付けることになるので，何でもかんでもリストワイズするのは危険かもしれない．

そこで多重代入法（Multivariate Imputation by Chained Equations; MICE）を使って欠損値の補完をやってみる．

miceパッケージを利用する．

データの成型

irisデータを使う．中身を見てみよう．

ついでに，今回使うパッケージを読み込んでおく．

pacman::p_load(tidyverse, mice, estimatr, modelsummary)
tibble(iris)
## # A tibble: 150 × 5
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           <dbl>       <dbl>        <dbl>       <dbl> <fct>  
##  1          5.1         3.5          1.4         0.2 setosa 
##  2          4.9         3            1.4         0.2 setosa 
##  3          4.7         3.2          1.3         0.2 setosa 
##  4          4.6         3.1          1.5         0.2 setosa 
##  5          5           3.6          1.4         0.2 setosa 
##  6          5.4         3.9          1.7         0.4 setosa 
##  7          4.6         3.4          1.4         0.3 setosa 
##  8          5           3.4          1.5         0.2 setosa 
##  9          4.4         2.9          1.4         0.2 setosa 
## 10          4.9         3.1          1.5         0.1 setosa 
## # ℹ 140 more rows

今回はSepal.Lengthをターゲットにしてやってみる．

Sepal.Lengthの値をランダムにNAに変換して，データフレームを作る．

NAの割合は，ひとまず20%としておこう．

data <- iris

set.seed(1231)

num_na <- floor(nrow(data) * 0.2)
na_indices <- sample(1:nrow(data), num_na)

data$Sepal.Length[na_indices] <- NA

NAの具合は以下の通り．

colSums(is.na(data))
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
##           30            0            0            0            0

miceで多重代入

では実際に，NAを補完していく．シード値の固定を忘れずに．

また，サンプルサイズが大きいと計算に時間がかかり，最悪の場合Rが爆発してそれまでの苦労が吹っ飛んでしまうらしい．

そういう場合は小さめのデータでテストしたり，変数を減らして負荷を軽減するなどの工夫をしたりすることが肝要らしい．

今回は全然大丈夫．

mice(data, m = 5, maxit = 100, seed = 12345) -> data_mice

mは補完データセット数（今回は5個）で，maxitは各欠測に対する補完の回数（今回は100回）となっている．

補完された値が初期値の影響を受けていないか確認する．

plot(data_mice)

いい感じに鎖が混ざっている．こちらについては，特に問題なさそう．

念のため，補完後のデータの分布に問題がないか確認しておく．

densityplot(data_mice)

青線が元データで，赤線が補完後データ（5つ）の分布となっている．

あんまり良くなさげだが，まあよしとしよう．

多重代入法で補完したデータでOLSをやる

補完後のデータを使って，Sepal.LengthについてOLSをしてみる．

with(data_mice, lm_robust(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species)) -> lm_mice

結果をみるときはこれ．ぶっちゃけ見にくい．

summary(pool(lm_mice))
##                term   estimate  std.error statistic       df      p.value
## 1       (Intercept)  2.1141131 0.29750022  7.106257 49.54150 4.291445e-09
## 2       Sepal.Width  0.5241021 0.10144322  5.166458 31.70339 1.255149e-05
## 3      Petal.Length  0.8435038 0.09677312  8.716303 17.67296 8.190570e-08
## 4       Petal.Width -0.4631615 0.18064955 -2.563867 62.68414 1.276167e-02
## 5 Speciesversicolor -0.6222502 0.28296995 -2.198997 52.21001 3.233130e-02
## 6  Speciesvirginica -0.8479610 0.40675129 -2.084716 57.80613 4.152276e-02

無理矢理キレイに結果をまとめる

上記のコードで結果を確認することはできるが，目がチカチカする．

そこでmsummaryを使って，無理矢理キレイにまとめてみようと思う（警告メッセージでmodelsummaryが文句を言ってくるが，コードは通るのでひとまず良しとする）．

ついでに，元のirisデータを使ったOLSの結果も載せて，多重代入版と比較する．

lm_robust(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, data = iris) -> lm_normal

summary <- list("普通のOLS" = lm_normal,
                "多重代入のOLS" = lm_mice)

msummary(summary, stars = c("*" = .05, "**" = .01, "***" = .001),
         notes = list('括弧内は標準誤差'))
普通のOLS 多重代入のOLS
* p < 0.05, ** p < 0.01, *** p < 0.001
括弧内は標準誤差
(Intercept) 2.171*** 2.114***
(0.249) (0.298)
Sepal.Width 0.496*** 0.524***
(0.081) (0.101)
Petal.Length 0.829*** 0.844***
(0.073) (0.097)
Petal.Width -0.315 -0.463*
(0.166) (0.181)
Speciesversicolor -0.724** -0.622*
(0.240) (0.283)
Speciesvirginica -1.023** -0.848*
(0.351) (0.407)
Num.Obs. 150
R2 0.867
R2 Adj. 0.863
AIC 79.1
BIC 100.2
RMSE 0.30

表示されていない箇所があるのはご愛嬌．

多重代入法と少し仲良くなれたかもしれない．