読者です 読者をやめる 読者になる 読者になる

疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R 欠損値の対応(missing value treatment)

今回は、欠損値の発生する理由を紹介し、その後に書籍を参考に欠損値に対応するRのコードを紹介します。
医学研究だけでなく、様々な調査をしていると欠損値(missing value)に出会います。欠損値が発生しているメカニズムによっては、結果を大きく変える可能性もあります。そのため、まずは欠損値のメカニズムを把握すること、その後欠損値に対応することになります。

欠損値が発生する理由

  1. MCAR (Missing completely at random)「完全にランダムな欠損値」
  2. MAR (Missing at random)「欠損値と他の変数には関連があるが、その変数の値自体との関連はない」
  3. MNAR (Missing not at random)「欠損値とその変数自体に関連がある」

欠損値のパターンとその対象法など詳しく書かれています。
norimune.net

欠損値への対処方法

「Useful R 2 データ分析プロセス」中で紹介されている欠損値への対応法は下記のものがある。

方法 概要
リストワイズ法 欠損値のあるサンプルを除外
ペアワイズ法 相関係数や分散の産出で、2変数のどちらかに欠損値のあるサンプルを除外
平均値代入法 平均値を補完する
回帰代入法 欠損値を除き、回帰分析をして推定値を補完
確率的回帰代入法 回帰代入法により推定した値に誤差をランダムに加えて補完
完全情報最尤推定 欠損値パターンに応じた個別の尤度関数を仮定した最尤推定により補完
多重代入法 欠損値に代入したデータセットを複数作成し、各データセットに対して分析を実行し、その結果を統合することにより欠損値を補完

Moving to R

実際にRで欠損値のパターンの認識と、それに対応するコードを紹介する。

サンプルデータ

10人の人でIDを含めて6つの変数を設定した。その中でいくつかにNA(欠損値)を導入している。

ID <- 1:10
height <- c(123, 134, 156, 147, 120, 135, 129, 119, 131, 125) 
weight <- c(35, 39, 43, 40, 29, NA, 29, 22, 25, NA)
age <- c(11, 10, 9, 8, 9, 9, 10, 11, 8, 8)
bp <-c(100, NA, 121, NA, 98, 111, 119, 101, 93, 130)
bs <- c(99, 100, 101, NA, 106, 80, 67, 88, 92, 79)
sample <- data.frame(ID=ID, Hei=height, Wei=weight, Age=age, BP=bp, BS=bs)
sample

##    ID Hei Wei Age  BP  BS
## 1   1 123  35  11 100  99
## 2   2 134  39  10  NA 100
## 3   3 156  43   9 121 101
## 4   4 147  40   8  NA  NA
## 5   5 120  29   9  98 106
## 6   6 135  NA   9 111  80
## 7   7 129  29  10 119  67
## 8   8 119  22  11 101  88
## 9   9 131  25   8  93  92
## 10 10 125  NA   8 130  79

1. 欠損値のパターンを理解する

欠損値のパターンの確認には、miceパッケージのmd.pattern関数やmd.pairs関数、VIMパッケージのmarginplot関数を用いる

1-1. miceパッケージ

md.pattern関数で欠損値の組み合わせを確認する。

md.pattern(sample)

##   ID Hei Age BS Wei BP  
## 6  1   1   1  1   1  1 0
## 2  1   1   1  1   0  1 1
## 1  1   1   1  1   1  0 1
## 1  1   1   1  0   1  0 2
##    0   0   0  1   2  2 5

この結果は、全ての変数に対して、それぞれの欠損値の有無の組み合わせ(1:欠損なし、0:欠損あり)を表示している。全ての項目で欠損なしが6名、weightだけ欠損ありが2名、BPだけ欠損ありが1名、BSとBPの両者で欠損ありが1名という読み方ができる。

md.pairs関数で2項目の欠損有無の組み合わせごとに件数を表示する

md.pairs(sample)

## $rr
##     ID Hei Wei Age BP BS
## ID  10  10   8  10  8  9
## Hei 10  10   8  10  8  9
## Wei  8   8   8   8  6  7
## Age 10  10   8  10  8  9
## BP   8   8   6   8  8  8
## BS   9   9   7   9  8  9
## 
## $rm
##     ID Hei Wei Age BP BS
## ID   0   0   2   0  2  1
## Hei  0   0   2   0  2  1
## Wei  0   0   0   0  2  1
## Age  0   0   2   0  2  1
## BP   0   0   2   0  0  0
## BS   0   0   2   0  1  0
## 
## $mr
##     ID Hei Wei Age BP BS
## ID   0   0   0   0  0  0
## Hei  0   0   0   0  0  0
## Wei  2   2   0   2  2  2
## Age  0   0   0   0  0  0
## BP   2   2   2   2  0  1
## BS   1   1   1   1  0  0
## 
## $mm
##     ID Hei Wei Age BP BS
## ID   0   0   0   0  0  0
## Hei  0   0   0   0  0  0
## Wei  0   0   2   0  0  0
## Age  0   0   0   0  0  0
## BP   0   0   0   0  2  1
## BS   0   0   0   0  1  1

rは非欠損値、mは欠損値、rmは行が非欠損値、列が欠損値を示している。weightとBSを例にすると、ともに観察したのは7名(rr)、weightのみ欠損しているのは1名(rm)、BSのみ欠損しているのは2名(mr)、両方とも欠損しているのは0名(mm)であった。

1-2. VIMパッケージ

aggr(sample, prop=FALSE, number=TRUE)

図の左は単一のデータ項目に対する欠損値の棒グラフ(個数)を、図の右は複数のデータ項目に対する欠損値のパターンを可視化している。md.pattern関数によるものと同等だが、視覚的に理解しやすいことがわかる。
f:id:ryosukefujii0320:20170203164713p:plain

#pbox関数で特定の項目とその他の項目間の関係について、欠損値の有無ごとに箱ひげ図をプロットする
#nhanesの項目の欠損有無によるheightの値の変化
pbox(sample, pos=2, int=FALSE, cex=1.2)

f:id:ryosukefujii0320:20170203164745p:plain

2. 欠損値に対処する

2-1. リストワイズ法(list wise)

少なくとも一つの変数が欠損しているデータを除去する方法である。この方法は欠損のメカニズムとしてMCAR(完全にランダムな欠損)を仮定している。MCAR以外の欠損値パターンであると、結果を大きく歪める可能性がある。Rではna.omit関数(popularな方法)を使用することで、実施できる。

#リストワイズ法前のデータ行列数
dim(sample)
 ## [1] 10  6

#リストワイズ法の実施
sample.lw <- na.omit(sample)
sample.lw
##   ID Hei Wei Age  BP  BS
## 1  1 123  35  11 100  99
## 3  3 156  43   9 121 101
## 5  5 120  29   9  98 106
## 7  7 129  29  10 119  67
## 8  8 119  22  11 101  88
## 9  9 131  25   8  93  92

#リストワイズ法後のデータ行列数
dim(sample.lw)
## [1] 6 6

2-2. ペアワイズ法(pair wise)

前述のリストワイズでは多くのデータの損失を伴うため、それを軽減するために開発された方法である。2つのデータ項目間の相関係数や共分散を求める場合に、いずれかが欠損値を持つサンプルを使用せずに計算を実施する。

cor(sample$Age, sample$BS) #NAを含んだ状態で計算しようとする
## [1] NA

cor(sample$Age, sample$BS, use="pairwise") #NAを除外して計算する
## [1] 0.07856771

2-3. 平均値代入法(mean imputation)

欠損値を補正するために、平均値を代入する。Rではmiceパッケージを使用することで実施可能。

imp <- mice(subset(sample, select=c(Wei, BP, BS)), method="mean", m=1, maxit=1)
imp

## Multiply imputed data set
## Call:
## mice(data = subset(sample, select = c(Wei, BP, BS)), m = 1, method = "mean", 
##     maxit = 1)
## Number of multiple imputations:  1
## Missing cells per column:
## Wei  BP  BS 
##   2   2   1 
## Imputation methods:
##    Wei     BP     BS 
## "mean" "mean" "mean" 
## VisitSequence:
## Wei  BP  BS 
##   1   2   3 
## PredictorMatrix:
##     Wei BP BS
## Wei   0  1  1
## BP    1  0  1
## BS    1  1  0
## Random generator seed value:  NA
#欠損値のある列に関して補完した行列の表示
complete(imp)

##      Wei      BP        BS
## 1  35.00 100.000  99.00000
## 2  39.00 109.125 100.00000
## 3  43.00 121.000 101.00000
## 4  40.00 109.125  90.22222
## 5  29.00  98.000 106.00000
## 6  32.75 111.000  80.00000
## 7  29.00 119.000  67.00000
## 8  22.00 101.000  88.00000
## 9  25.00  93.000  92.00000
## 10 32.75 130.000  79.00000

#元のデータを編集
sample.mean <- sample
sample.mean$Wei <- complete(imp)$Wei
sample.mean$BP <- complete(imp)$BP
sample.mean$BS <- complete(imp)$BS

#完成したデータの表示
sample.mean

##    ID Hei   Wei Age      BP        BS
## 1   1 123 35.00  11 100.000  99.00000
## 2   2 134 39.00  10 109.125 100.00000
## 3   3 156 43.00   9 121.000 101.00000
## 4   4 147 40.00   8 109.125  90.22222
## 5   5 120 29.00   9  98.000 106.00000
## 6   6 135 32.75   9 111.000  80.00000
## 7   7 129 29.00  10 119.000  67.00000
## 8   8 119 22.00  11 101.000  88.00000
## 9   9 131 25.00   8  93.000  92.00000
## 10 10 125 32.75   8 130.000  79.00000

年齢と血糖を例にして、リストワイズ法およびペアワイズ法、平均代入法を用いた場合を比較する。

#リストワイズ法
cor(sample.lw$Age, sample.lw$BS)
## [1] -0.2335631

#ペアワイズ法
cor(sample$Age, sample$BS, use="pairwise")
## [1] 0.07856771

#平均値代入法
cor(sample.mean$Age, sample.mean$BS)
## [1] 0.07221444

2-4. 回帰代入法(regression imputation)

#回帰式の推定
fit.lm <- lm(BS ~ Age, data=sample)
pred <- predict(fit.lm, subset(sample, is.na(BS)))

#欠損値の補完
sample.regression <- sample
sample.regression$BS[is.na(sample$BS)] <- pred

#結果の確認
cor(sample.regression$Age, sample.regression$BS)
## [1] 0.0854315

2-5. 多重代入法(multiple imputation)

欠損値を代入したデータを複数作成し、それぞれについて分析した後、その結果を統合する つまり、代入、分析、統合の3ステップからなっている Rでは、miceパッケージおよびmiパッケージを使用できる

#欠損値を補完したデータセットの作成
imp <- mice(subset(sample, select=c(BP, BS)), seed=123, print=FALSE)
#分析
fit <- with(imp, lm(BS ~ Age, data=sample))
#分析結果の統合
pool.fit <- pool(fit)
sum.pf <- summary(pool.fit)
tab <- round(sum.pf, 3)
tab[,c(1:3, 5)]

##                est     se     t Pr(>|t|)
## (Intercept) 81.804 40.627 2.014    0.094
## Age          0.891  4.275 0.209    0.842
#欠損値の補完
slope <- pool.fit$qbar[2]
intercept <- pool.fit$qbar[1]
imputed <- (slope*sample$Age + intercept)
sample.mi <- sample
is.missing <- is.na(sample.mi$BS)
sample.mi$BS[is.missing] <- imputed[is.missing]
mean(sample$BS, na.rm=TRUE)
## [1] 90.22222

mean(sample.mi$BS)
## [1] 90.09348

参考文献

Useful R(2) データ分析プロセス P.46-60
www.kyoritsu-pub.co.jp


20170203
RF