欠損値の種類と補完とRでのワークフロー
1. はじめに
実験データを扱う時に欠損値が出てしまうことはよくあると思います。欠損値を除いて解析をしたりすることも多いと思いますが、どうしても欠損値を含むデータを使用しなければいけないこともあると思います。欠損値を補完する方法はいろいろと開発されているようです。
http://fixxman.hatenablog.com/entry/2015/10/06/045417fixxman.hatenablog.com
途中まで書いていて、まったく同じ記事で質が段違いの記事があるじゃないかと思いましたが、殆ど書いてしまったので引用させてもらいながら続けます。
2. 欠損データのパターン
データが欠損する際にはいくつかパターンが存在します。Sunny side up!によると
- 完全にランダムに欠損 -> Missing Completely At Random ( MCAR )
- 観測データに依存する欠損 ->Missing At Random ( MAR )
体重が大きい人が体重を測ることを拒否- 欠損データに依存する欠損 -> Missing Not At Random ( MNAR )
女性が生理などで欠席することが多く、そのことによる欠損
とあります。実験のミスなどで欠損してしまうのは1.ということになりますかね。
3. 欠損値の保管方法
このセクションはSE Can't Codeから引用しています。
リストワイズ法: 欠損値をもつサンプルの削除
ペアワイズ法: 相関係数や分散等の算出において、2変数のいずれかが欠損値をもつサンプルを削除
平均値代入法: 平均値により欠損値を補完
回帰代入法: 欠損値のないサンプルに回帰分析を行い、欠損値を含む項目の推定式を元に欠損値を補完
確率的回帰代入法: 回帰代入法により確定した値にランダムに誤差を加えて欠損値を補完
完全情報最尤推定法: サンプル毎に欠損パターンに応じた尤度関数を仮定して最尤推定を実施して得られる多変量正規分布を用いて平均値や分散共分散行列を推定
多重代入法: 欠損値に代入したデータセットを複数作成し、各データセットに対して分析を実行し、その結果を統合することにより欠損値を補完
4. パターン毎の補完法
Sunny side up!を参照。ざっくりまとめると
MCAR
- ランダムが偶然に起きているのでどのような対処をしても母集団への推定に影響しない
- ペアワイズ削除やリストワイズ削除でも母数の推定に影響しない
- 全部のデータを使わないリストワイズ削除は単純に「もったいない」
- 平均値代入などの単一の値を欠損値に代入する法
- MCARであっても推定値にバイアスを受ける
- 推定精度を過剰によく推定しまうためタイプⅠエラーを犯す危険性
- 基本的には使わないほうがいい
MAR
- リストワイズ削除でも推定にバイアスが生じる
- 特定の傾向を持つ人がデータからいなくなってしまうため
- 解決策
- 観測データをすべて用いた推定方法を使う
- 多重代入による推定を行う
- この二つの推定精度やバイアスはほぼ同程度だと言われている
1. 完全情報最尤法 ( Full Information Maximum Likelihood method )
- FIMLと呼ばれたり、ただ単に最尤法と呼ばれたりする
2. 多重代入法
MNAR
- 有効な方法はあまりない
5. Rでの実際の作業
じゃあ実際にRでどう作業していくかということを記載します。この記事のスクリプトおよび解説は下記サイトから使用させていただいております。 datascienceplus.com
データの作成
data <- airquality # airqualityにNAを強制挿入する data[4:10,3] <- rep(NA,7) data[1:5,4] <- NA data <- data[-c(5,6)] # summaryを使用すると各variableのNAの数が表示される summary(data)
{mice} と {VIM} を使用した全体像の把握
欠損値を含むデータの場合、miceパッケージが便利です。まずmd.pattern()で欠損値を含むデータの全体像を知りましょう。
# miceをinstall
install.packages("mice")
library(mice)
md.pattern(data)
Temp Solar.R Wind Ozone
104 1 1 1 1 0
34 1 1 1 0 1
4 1 0 1 1 1
3 1 1 0 1 1
3 0 1 1 1 1
1 1 0 1 0 2
1 1 1 0 0 2
1 1 0 0 1 2
1 0 1 0 1 2
1 0 0 0 0 4
5 7 7 37 56
104個はすべて揃っているデータで、34個はOzoneに欠損値があるデータで、4個はWindに欠損値が有るデータという様に結果を返してくれます。 さらにこれをVIMパッケージを使用することで視覚化することができます。VIMはVisualization and Imputation of Missing Valuesの略で、欠損値を含むデータの可視化のために作られたパッケージです。
install.packages("VIM")
library(VIM)
aggr_plot <- aggr(data,
col=c('navyblue','red'),
numbers=TRUE, sortVars=TRUE, labels=names(data),
cex.axis=.7, gap=3,
ylab=c("Histogram of missing data","Pattern"))
marginplot(data[c(1,2)])
データを補完する
mice()を使用します。データセットの候補をいくつ設定するかは-mで設定します。
tempData <- mice(data,
m=10, # refers to the number of imputed datasets. Five is the default value.
maxit=50,
meth='pmm', # refers to the imputation method, pmm: predictive mean matching
seed=500)
summary(tempData)
ちなみにmethods(mice)で補完する方法を列挙してくれます。
methods(mice) [1] mice.impute.2l.norm mice.impute.2l.pan mice.impute.2lonly.mean mice.impute.2lonly.norm [5] mice.impute.2lonly.pmm mice.impute.cart mice.impute.fastpmm mice.impute.lda [9] mice.impute.logreg mice.impute.logreg.boot mice.impute.mean mice.impute.norm [13] mice.impute.norm.boot mice.impute.norm.nob mice.impute.norm.predict mice.impute.passive [17] mice.impute.pmm mice.impute.polr mice.impute.polyreg mice.impute.quadratic [21] mice.impute.rf mice.impute.ri mice.impute.sample mice.mids [25] mice.theme
補完したデータを確認するには
tempData$imp$Ozone
とします。行が欠損値のvariable、列が作成した補完データセットの種類です。mice()の-mを10にしておけば10列作成されます。この中から好きな ( ? ) データを選んでデータをcomplete()でデータを確定させます。
completedData <- complete(tempData,3) # 選択した番号
補完後のデータの確認を行う
保管したデータ変な値になっていないかどうかチェックをしておく必要があります。
Scatter plot
まずはxyplot()でscatter plotで全体像を把握します。赤が補完したデータで青が本来計測されたデータです。だいたい分布が一緒になっており最もらしい補完データであることがわかります。
library(lattice) xyplot(tempData,Ozone ~ Wind+Temp+Solar.R,pch=18,cex=1)
Density plot
作成したデータセット分の補完値のdensity plotが描画される。赤が補完値、青が実測値。
densityplot(tempData)
Strip plot
同上
stripplot(tempData, pch = 20, cex = 1.2)
