1 はじめに

データセットに欠損値（Missing Data）が含まれることは、統計解析において一般的である。
欠損値を適切に処理しないと分析結果にバイアスが生じる可能性がある。
主にRのmiceパッケージを使用して、欠損値を処理する方法を記載する。
miceパッケージのサイトhttps://amices.org/mice/index.htmlを参考にしている。

2 欠損の種類

MCAR (Missing Completely at Random): 欠損が完全にランダムで、他の変数やデータに依存しない。
MAR (Missing at Random): 欠損が観測された他の変数に依存するが、欠損値そのものには依存しない。
MNAR (Missing Not at Random): 欠損が欠損値そのものに依存する（最も処理が難しい）。

3 欠損データの特定と可視化

3.1 欠損値の特定

Rの関数is.na()やsummary()を使用して、データセット内の欠損値を確認できる。
- 例: sum(is.na(data))で欠損値の総数を計算。

3.2 可視化

パッケージmiceやVIMを使用して、欠損パターンを視覚的に把握。
- 例: VIM::aggr()で欠損値の割合やパターンをプロット。
これにより、どの変数に欠損が多いか、欠損がランダムか系統的かを評価可能

4 欠損データ処理の基本手法

4.1 完全ケース解析（Listwise Deletion）

欠損値を含む行を削除（na.omit()）。
メリット: 簡単で迅速。
デメリット: データ量が減少し、バイアスが生じる可能性（特にMCARでない場合）。

4.2 単一代入（Single Imputation）

欠損値を平均値、中央値、最頻値などで埋める。
例: impute.mean()（Hmiscパッケージ）で平均値代入。
デメリット: 欠損値の不確実性を考慮せず、標準誤差が過小評価される。

4.3 多重代入（Multiple Imputation, MI）

多くの場合、推奨される手法。
欠損値を複数回異なる値で埋め、複数データセットを生成。
不確実性を考慮し、ルビンの規則（Rubin’s Rules）で結果を統合。
例: miceパッケージを使用。

5 ルビンの規則

ルビンの規則（Rubin’s Rules）は、多重代入（Multiple Imputation, MI）において、複数の代入データセットから得られた分析結果を統合（Pooling）するための統計的手法。
この手法は、欠損値の不確実性を適切に考慮し、点推定値（例: 回帰係数）とその標準誤差を計算するために使用。

基本的な考え方

多重代入では、欠損値を複数回（例: ( m ) 個のデータセット）異なる値で埋めて、複数の完全データセットを生成。
各データセットで個別に分析（例: 回帰分析）を行い、結果（例: 係数や標準誤差）を統合する必要がある。
ルビンの規則は以下の統合を行う。

点推定値の統合: 各データセットの推定値（例: 回帰係数）の平均を取る。
標準誤差の計算: データセット内での分散（Within-Imputation Variance）とデータセット間の分散（Between-Imputation Variance）を組み合わせて、欠損値の不確実性を反映。
統計的推論: 統合された推定値に基づく信頼区間やp値を計算。

6 Rでの多重代入

6.1 `mice`パッケージの概要

mice（Multivariate Imputation by Chained Equations）パッケージは、連鎖方程式による多重代入（MICE）を実装。
各変数の条件付き分布をモデル化し、反復的に欠損値を代入。

6.2 パッケージのインストールとロード

# install.packages("mice")
library(mice)

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

6.3 データセットの準備と確認

欠損値がある観測は、mice::md.pattern()で確認できる。
nhanesデータ
- age（年齢）：欠損なし
- bmi（BMI）：欠損あり
- hyp（高血圧）：欠損あり
- chl（コレステロール）：欠損あり

一部の変数に欠損値が存在。

# データセットの読み込み
data(nhanes, package="mice")

# データセットの確認
mice::md.pattern(nhanes)

##    age hyp bmi chl   
## 13   1   1   1   1  0
## 3    1   1   1   0  1
## 1    1   1   0   1  1
## 1    1   0   0   1  2
## 7    1   0   0   0  3
##      0   8   9  10 27

例えば、以下のようなことがわかる：

欠損値が0の観測が、13観測値ある
欠損値が1の観測が、3観測値ある
欠損値が3の観測が、3観測値ある

6.4 代入の実行

mice()関数で多重代入を実行
3つのデータセットを生成。
mice()関数の主な引数は以下の通り：
- m: 生成するデータセット数。
- maxit: 反復回数。
- seed: 再現性のための乱数シード。
method: 代入方法。例: “pmm”は予測平均マッチング（Predictive Mean Matching）。PMMは、連続変数に適している。

imp <- mice::mice(nhanes, m=3, maxit=50, seed=500, method="pmm")

## 
##  iter imp variable
##   1   1  bmi  hyp  chl
##   1   2  bmi  hyp  chl
##   1   3  bmi  hyp  chl
##   2   1  bmi  hyp  chl
##   2   2  bmi  hyp  chl
##   2   3  bmi  hyp  chl
##   3   1  bmi  hyp  chl
##   3   2  bmi  hyp  chl
##   3   3  bmi  hyp  chl
##   4   1  bmi  hyp  chl
##   4   2  bmi  hyp  chl
##   4   3  bmi  hyp  chl
##   5   1  bmi  hyp  chl
##   5   2  bmi  hyp  chl
##   5   3  bmi  hyp  chl
##   6   1  bmi  hyp  chl
##   6   2  bmi  hyp  chl
##   6   3  bmi  hyp  chl
##   7   1  bmi  hyp  chl
##   7   2  bmi  hyp  chl
##   7   3  bmi  hyp  chl
##   8   1  bmi  hyp  chl
##   8   2  bmi  hyp  chl
##   8   3  bmi  hyp  chl
##   9   1  bmi  hyp  chl
##   9   2  bmi  hyp  chl
##   9   3  bmi  hyp  chl
##   10   1  bmi  hyp  chl
##   10   2  bmi  hyp  chl
##   10   3  bmi  hyp  chl
##   11   1  bmi  hyp  chl
##   11   2  bmi  hyp  chl
##   11   3  bmi  hyp  chl
##   12   1  bmi  hyp  chl
##   12   2  bmi  hyp  chl
##   12   3  bmi  hyp  chl
##   13   1  bmi  hyp  chl
##   13   2  bmi  hyp  chl
##   13   3  bmi  hyp  chl
##   14   1  bmi  hyp  chl
##   14   2  bmi  hyp  chl
##   14   3  bmi  hyp  chl
##   15   1  bmi  hyp  chl
##   15   2  bmi  hyp  chl
##   15   3  bmi  hyp  chl
##   16   1  bmi  hyp  chl
##   16   2  bmi  hyp  chl
##   16   3  bmi  hyp  chl
##   17   1  bmi  hyp  chl
##   17   2  bmi  hyp  chl
##   17   3  bmi  hyp  chl
##   18   1  bmi  hyp  chl
##   18   2  bmi  hyp  chl
##   18   3  bmi  hyp  chl
##   19   1  bmi  hyp  chl
##   19   2  bmi  hyp  chl
##   19   3  bmi  hyp  chl
##   20   1  bmi  hyp  chl
##   20   2  bmi  hyp  chl
##   20   3  bmi  hyp  chl
##   21   1  bmi  hyp  chl
##   21   2  bmi  hyp  chl
##   21   3  bmi  hyp  chl
##   22   1  bmi  hyp  chl
##   22   2  bmi  hyp  chl
##   22   3  bmi  hyp  chl
##   23   1  bmi  hyp  chl
##   23   2  bmi  hyp  chl
##   23   3  bmi  hyp  chl
##   24   1  bmi  hyp  chl
##   24   2  bmi  hyp  chl
##   24   3  bmi  hyp  chl
##   25   1  bmi  hyp  chl
##   25   2  bmi  hyp  chl
##   25   3  bmi  hyp  chl
##   26   1  bmi  hyp  chl
##   26   2  bmi  hyp  chl
##   26   3  bmi  hyp  chl
##   27   1  bmi  hyp  chl
##   27   2  bmi  hyp  chl
##   27   3  bmi  hyp  chl
##   28   1  bmi  hyp  chl
##   28   2  bmi  hyp  chl
##   28   3  bmi  hyp  chl
##   29   1  bmi  hyp  chl
##   29   2  bmi  hyp  chl
##   29   3  bmi  hyp  chl
##   30   1  bmi  hyp  chl
##   30   2  bmi  hyp  chl
##   30   3  bmi  hyp  chl
##   31   1  bmi  hyp  chl
##   31   2  bmi  hyp  chl
##   31   3  bmi  hyp  chl
##   32   1  bmi  hyp  chl
##   32   2  bmi  hyp  chl
##   32   3  bmi  hyp  chl
##   33   1  bmi  hyp  chl
##   33   2  bmi  hyp  chl
##   33   3  bmi  hyp  chl
##   34   1  bmi  hyp  chl
##   34   2  bmi  hyp  chl
##   34   3  bmi  hyp  chl
##   35   1  bmi  hyp  chl
##   35   2  bmi  hyp  chl
##   35   3  bmi  hyp  chl
##   36   1  bmi  hyp  chl
##   36   2  bmi  hyp  chl
##   36   3  bmi  hyp  chl
##   37   1  bmi  hyp  chl
##   37   2  bmi  hyp  chl
##   37   3  bmi  hyp  chl
##   38   1  bmi  hyp  chl
##   38   2  bmi  hyp  chl
##   38   3  bmi  hyp  chl
##   39   1  bmi  hyp  chl
##   39   2  bmi  hyp  chl
##   39   3  bmi  hyp  chl
##   40   1  bmi  hyp  chl
##   40   2  bmi  hyp  chl
##   40   3  bmi  hyp  chl
##   41   1  bmi  hyp  chl
##   41   2  bmi  hyp  chl
##   41   3  bmi  hyp  chl
##   42   1  bmi  hyp  chl
##   42   2  bmi  hyp  chl
##   42   3  bmi  hyp  chl
##   43   1  bmi  hyp  chl
##   43   2  bmi  hyp  chl
##   43   3  bmi  hyp  chl
##   44   1  bmi  hyp  chl
##   44   2  bmi  hyp  chl
##   44   3  bmi  hyp  chl
##   45   1  bmi  hyp  chl
##   45   2  bmi  hyp  chl
##   45   3  bmi  hyp  chl
##   46   1  bmi  hyp  chl
##   46   2  bmi  hyp  chl
##   46   3  bmi  hyp  chl
##   47   1  bmi  hyp  chl
##   47   2  bmi  hyp  chl
##   47   3  bmi  hyp  chl
##   48   1  bmi  hyp  chl
##   48   2  bmi  hyp  chl
##   48   3  bmi  hyp  chl
##   49   1  bmi  hyp  chl
##   49   2  bmi  hyp  chl
##   49   3  bmi  hyp  chl
##   50   1  bmi  hyp  chl
##   50   2  bmi  hyp  chl
##   50   3  bmi  hyp  chl

代入されたbmiの値を確認する。
観測ごとに、3個の代入値がある。

imp$imp$bmi

##       1    2    3
## 1  22.0 30.1 28.7
## 3  27.2 27.2 28.7
## 4  24.9 20.4 20.4
## 6  22.5 21.7 24.9
## 10 20.4 28.7 33.2
## 11 30.1 33.2 30.1
## 12 26.3 22.5 27.4
## 16 22.7 22.5 30.1
## 21 20.4 33.2 22.5

6.5 代入データの出力

complete()関数で、代入されたデータセットを取得。
complete(imp, 1)で1つ目の代入データセットを取得。
complete(imp, "long")で全ての代入データセットを長形式で取得。
complete(imp, "long", include = TRUE)で、元のデータセットと代入データセットを結合。
include = TRUEを指定すると、元のデータセットと代入データセットを結合したデータフレームが得られる。

# 代入データの出力
complete1 <- complete(imp, 1)
complete2 <- complete(imp, 2)
complete3 <- complete(imp, 3)

completeAll <- complete(imp, "long")
completeAll2 <- complete(imp, "long", include = TRUE)

# 代入データの確認
head(complete1)

##   age  bmi hyp chl
## 1   1 22.0   1 131
## 2   2 22.7   1 187
## 3   1 27.2   1 187
## 4   3 24.9   2 204
## 5   1 20.4   1 113
## 6   3 22.5   1 184

head(completeAll)

##   age  bmi hyp chl .imp .id
## 1   1 22.0   1 131    1   1
## 2   2 22.7   1 187    1   2
## 3   1 27.2   1 187    1   3
## 4   3 24.9   2 204    1   4
## 5   1 20.4   1 113    1   5
## 6   3 22.5   1 184    1   6

6.6 代入データの確認

summary(imp)で代入結果を確認。
stripplot(imp)で代入値の分布を可視化。

# 代入結果の確認
summary(imp)

## Class: mids
## Number of multiple imputations:  3 
## Imputation methods:
##   age   bmi   hyp   chl 
##    "" "pmm" "pmm" "pmm" 
## PredictorMatrix:
##     age bmi hyp chl
## age   0   1   1   1
## bmi   1   0   1   1
## hyp   1   1   0   1
## chl   1   1   1   0

ここで、PredictorMatrixの値は以下のように解釈される：

1: その変数が他の変数の代入に使用された（例: ageはbmi, hyp, chlの代入に使用）。
0: 使用されなかった（例: bmiは自身の代入には使用しない）。
stripplot()は、代入値の分布を可視化し、元の観測値と比較。
これにより、代入値が妥当かどうかを評価。
赤い点: 代入値（Imputed values）。
青い点: 観測値（Observed values）。

stripplot(imp)

6.7 分析と統合

with()で各データセットにモデルを適用（例: 線形回帰）。
各データセットごとに推定結果を得る。
with()関数は、各データセットに対して指定されたモデルを適用し、結果をリストとして返す
with()関数の主な引数は以下の通り：
- data: データセット。
- expr: 適用するモデル（例: 線形回帰）。。

# 線形回帰モデルの適用
fit <- with(imp, lm(bmi ~ hyp + chl))
summary(fit)

## # A tibble: 9 × 7
##   term        estimate std.error statistic   p.value  nobs df.residual
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl> <int>       <dbl>
## 1 (Intercept)  19.2       3.87      4.95   0.0000590    25          22
## 2 hyp          -1.38      1.86     -0.739  0.468        25          22
## 3 chl           0.0425    0.0194    2.19   0.0392       25          22
## 4 (Intercept)  22.3       4.28      5.20   0.0000326    25          22
## 5 hyp          -2.81      2.08     -1.35   0.190        25          22
## 6 chl           0.0423    0.0240    1.76   0.0916       25          22
## 7 (Intercept)  23.7       4.52      5.25   0.0000287    25          22
## 8 hyp          -0.211     2.18     -0.0967 0.924        25          22
## 9 chl           0.0170    0.0219    0.777  0.446        25          22

6.8 統合結果の取得

pool()関数でルビンの規則に基づき、各データセットの結果を統合。
summary()で統合結果を確認。

# pool()でルビンの規則に基づき結果を統合。

pooled <- mice::pool(fit)
summary(pooled)

##          term    estimate  std.error  statistic       df     p.value
## 1 (Intercept) 21.72261156 5.01569854  4.3309245 9.019006 0.001893146
## 2         hyp -1.46467162 2.53528740 -0.5777142 7.271823 0.580897760
## 3         chl  0.03392592 0.02761348  1.2285996 6.699584 0.260624683

7 カテゴリー変数の多重代入

nhanesのhyp（バイナリ変数）をlogregで代入する場合
method引数にpmm（連続変数）とlogreg（バイナリ変数）を指定。

# hypをカテゴリー変数に変換
data(nhanes, package="mice")
nhanes2 <- nhanes
nhanes2$hyp <- as.factor(nhanes2$hyp)

imp2 <- mice(nhanes2, m=3, 
            method=c("", "pmm", "logreg", "pmm"), 
            seed=500)

## 
##  iter imp variable
##   1   1  bmi  hyp  chl
##   1   2  bmi  hyp  chl
##   1   3  bmi  hyp  chl
##   2   1  bmi  hyp  chl
##   2   2  bmi  hyp  chl
##   2   3  bmi  hyp  chl
##   3   1  bmi  hyp  chl
##   3   2  bmi  hyp  chl
##   3   3  bmi  hyp  chl
##   4   1  bmi  hyp  chl
##   4   2  bmi  hyp  chl
##   4   3  bmi  hyp  chl
##   5   1  bmi  hyp  chl
##   5   2  bmi  hyp  chl
##   5   3  bmi  hyp  chl

stripplot(imp2, hyp ~ .imp)

8 Stataでの多重代入

Royston (2009)の論文では、Stataのiceパッケージを用いたMICEの実装が解説されている。

参考文献

Princeton University Library. Missing Data Imputation in R: Missing data R tutorial: a brief guide to handle missing data in R. https://libguides.princeton.edu/c.php?g=1391216&p=10290937#s-lg-box-wrapper-38134556
Hayati Rezvan, P., Lee, K. J., & Simpson, J. A. (2015). The rise of multiple imputation: A review of the reporting and implementation of the method in medical research. BMC Medical Research Methodology, 15(1), 30. https://doi.org/10.1186/s12874-015-0022-1
Little, R. J., & Rubin, D. B. (2019). Statistical analysis with missing data. John Wiley & Sons.
Royston, P. (2009). Multiple Imputation of Missing Values: Further Update of Ice, with an Emphasis on Categorical Variables. The Stata Journal, 9(3), 466–477. https://doi.org/10.1177/1536867X0900900308
Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581-592. https://doi.org/10.1093/biomet/63.3.581
Sterne, J. A. C., White, I. R., Carlin, J. B., Spratt, M., Royston, P., Kenward, M. G., Wood, A. M., & Carpenter, J. R. (2009). Multiple imputation for missing data in epidemiological and clinical research: Potential and pitfalls. BMJ, 338(jun29 1), b2393–b2393. https://doi.org/10.1136/bmj.b2393

欠損値の処置

田中鮎夢

2025-04-22

1 はじめに

2 欠損の種類

3 欠損データの特定と可視化

3.1 欠損値の特定

3.2 可視化

4 欠損データ処理の基本手法

4.1 完全ケース解析（Listwise Deletion）

4.2 単一代入（Single Imputation）

4.3 多重代入（Multiple Imputation, MI）

5 ルビンの規則

6 Rでの多重代入

6.1 `mice`パッケージの概要

6.2 パッケージのインストールとロード

6.3 データセットの準備と確認

6.4 代入の実行

6.5 代入データの出力

6.6 代入データの確認

6.7 分析と統合

6.8 統合結果の取得

7 カテゴリー変数の多重代入

8 Stataでの多重代入

参考文献

欠損値の処置

田中 鮎夢

2025-04-22

1 はじめに

2 欠損の種類

3 欠損データの特定と可視化

3.1 欠損値の特定

3.2 可視化

4 欠損データ処理の基本手法

4.1 完全ケース解析（Listwise Deletion）

4.2 単一代入（Single Imputation）

4.3 多重代入（Multiple Imputation, MI）

5 ルビンの規則

6 Rでの多重代入

6.1 miceパッケージの概要

6.2 パッケージのインストールとロード

6.3 データセットの準備と確認

6.4 代入の実行

6.5 代入データの出力

6.6 代入データの確認

6.7 分析と統合

6.8 統合結果の取得

7 カテゴリー変数の多重代入

8 Stataでの多重代入

参考文献

田中鮎夢

6.1 `mice`パッケージの概要