kaggleで予測モデルを構築してみた (5) - Rで行うMultipleImputation

ということで，前回で触れたように，データの前処理を実際に行っていきたいと思います．その中でも今回は，欠損値補完についての話をしていきます．

今回のデータでは，NAが含まれているageのデータを補完する必要があります．とはいえ，欠損値を補完するにもいくつか方法があって，どの補完を行うのが妥当かというのを考えなければいけません．そこでまず，欠損値がどういう性質を持っているかについてみていきましょう．

欠損のメカニズム

欠損のパターン

データの欠損には大きく分けて3つのパターンがあります．

Missing Completely At Random(MCAR)：完全にランダムに欠損が生じているもの
Missing At Random(MAR) ：データ欠損が，データに含まれるほかの変数と関連はしているが，その影響を取り除いた自分自身の値とは関連していないもの
Missing Not At Random(MNAR)：上記の2例に当てはまらないもの

欠損値補完の方法

これらのパターンによって，どのような補完を行うのが妥当かが変わってきます．ではどのような補完方法があるかということになりますが，大きく分けると以下の通りです．

リストワイズ法：1つでも欠損値のあるサンプルは分析から除外
ペアワイズ法：分析対象の変数に欠損がない場合はすべてサンプルに含める
単一代入法：なんらかの値で補完する（1サンプルにつき1つの補完値）
1. 平均値・中央値など
2. 回帰による推定値：データ内のほかの変数から回帰分析で予測した値で補完する
3. 完全情報最尤推定（Full-Information Maximum Likelihood: FIML）による推定値：EMアルゴリズムを用いて推定した値で補完する
多重代入法（Multiple Imputation: MI）：欠損値に異なる値を補完した複数のデータセットを作成して，そのすべてを用いて分析を行い，その結果を統合する

どれを使うか

個人的な印象では，だいたいこんな感じです．できることならFIMLかMIを，そうでなくても回帰推定食らいはしておいた方がよくて，そうでなければリストワイズを使うのが手っ取り早いと思います．Rに限らず，SASとかSPSSとかだいたいの統計ソフトは何もしなければリストワイズを使います．

手法	ポイント	MCAR	MAR	MNAR
リストワイズ法	データ量減る．偏りに弱い	△	×	×
ペアワイズ法	計算面倒．偏りに弱い	△	×	×
平均値・中央値など	そもそも補完の意味があまりない	△	×	×
回帰による推定	偏りに弱い	○	△	△
FIML	今のところイケてる	◎	◎	△
MI	今のところイケてる	◎	◎	△

2000年代以降の応用統計学の界隈では，おおむねFIMLかMIによる欠損値補完が，最も実用的に優れているものといわれているようです．もちろん最先端の事例ではもっとあるんでしょうが，パッケージにあるものを簡便に使えるという意味では，FIMLかMIのどちらかを使うのがよいと思います．そのどちらにしても，ちゃんと使うのは割と面倒なんですけどね...

ということで，FIMLとMIについてもう少し詳しくみていきましょう．

Multiple Imputation

RでMIを行うためには，miパッケージを使えばできるそうです．詳細は以下の論文に全部書いてあるので，詳しくは読んでください．

http://www.jstatsoft.org/v45/i02/paper

ということで，実際にデータを補完してみました．コードは以下の通りです．

# パッケージのインストールとライブラリ呼び出し
install.packages('mi')
library('mi')

# データの読み込みと不要変数の削除
train <- read.table('~/Documents/workspace/R/kaggle/data/titanic/train.csv', header=T, sep=',')
train.mi <- train
train.mi$name <- NULL
train.mi$ticket <- NULL
train.mi$cabin <- NULL

# NAデータのプロット
mp.plot(train.mi, clustered=FALSE) 

# データの情報取得と確認
train.info <- mi.info(train.mi)
train.info$imp.formula

# 欠損値補完の前処理と確認
train.pre <- mi.preprocess(train.mi)
attr(train.pre, 'mi.info')

# 欠損値補完とデータの取得
train.imp <- mi(train.pre, R.hat=1.8)
train.dat.all <- mi.completed(train.imp)

MIはデフォルトでは全変数を予測に使ってしまうので，必要のない変数ははじいておいてください．今回の場合だとname, ticket, cabinの3変数は，年齢の予測に役立ちそうもないので削除しておきます（特に名前のような変数は，すべての値が違うカテゴリ変数として認識されてしまうので，MIの計算に非常に時間がかかることになります）．train.info$imp.formulaで確認できるように，基本のモデル式は全変数を使った回帰になります．

また，MIでの予測に適した形に前変換をかけたうえで，欠損値補完を行います．その結果はtrain.dat.allとして得ることができます．しかしMIはそれぞれ異なる値を補完した複数のデータセットを作成するため，train.dat.allをみれば3つのデータセットが含まれているのがわかると思います（いくつのデータセットを作成するかは変更可能）．

これらのうちどれかひとつだけ得られればよいのであれば．以下のようにしてデータフレームとして切り出すことができます．

train.dat <- mi.data.frame(train.imp, m=1)

この中身をみてみれば，

id	survived	pclass	sex	age	sibsp	fare	embarked
1	0	3	male	22.00	1	7.2500	S
2	1	1	female	38.00	1	71.2833	C
3	1	3	female	26.00	0	7.9250	S
4	1	1	female	35.00	1	53.1000	S
5	0	3	male	35.00	0	8.0500	S
6	0	3	male	NA	0	8.4583	Q
7	0	1	male	54.00	0	51.8625	S

だったのが，

id	survived	pclass	sex	age	sibsp	fare	embarked
1	0	3	male	22.00	1	7.2500	S
2	1	1	female	38.00	1	71.2833	C
3	1	3	female	26.00	0	7.9250	S
4	1	1	female	35.00	1	53.1000	S
5	0	3	male	35.00	0	8.0500	S
6	0	3	male	46.966050	0	8.4583	Q
7	0	1	male	54.00	0	51.8625	S