[python] statsmodelsでMICEを用いた重回帰分析を行う

MICEはMultiple Imputation by Chained Equationの略.
近年, よく使われる欠損値の処理の方法の一つ.
欠損値が一つでもあったら, そのデータを弾くリストワイズ法やペアワイズ法では, バイアスがかかってしまう状況でも処理できる方法.

この記事では, MICEの実装をstatsmodels[1]で行う.
だが, 残念なことに, 非常に残念なことに現状では, statsmodelsの中ではカテゴリカルのデータに対するコードが出来ておらず, ガウス分布での適応のみ.

そんな状況なわけですが, statsmodelsでの使い方紹介…
理論については, 統計数理研究所の野間久史(2017)の”連鎖方程式による多重代入法”が入門者には分かりやすかった.

準備
実装とグラフ化
おまけ
参考文献

準備

動作はjupyter notebook上を想定.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(context="paper" , style ="whitegrid",rc={"figure.facecolor":"white"})
# for MICE and multiple regression model
from statsmodels.imputation import mice
import statsmodels.api as sm

データはstatsmodelsが提供するもの.
データを構造化された形に変形.
今回は, さらにデータに欠損値を各項目に5つずつ挿入.

Data = sm.datasets.spector.load(as_pandas=True)
df = Data.exog
df["res"] = Data.endog

np.random.seed(10)
for k in df.keys():
    rand = np.random.randint(0,32,size = 5)
    df.loc[rand,k] = np.nan
print(df.shape)
df.head().T

＊ここで注意！ patsyのversionを”0.5.1″にアップデートしておかないと[2]の記事のようなエラーが起こります.

実装とグラフ化

式はR形式の書き方で書くことになる. ここでは, 切片が自動的に挿入される.

mice_.fitの第一引数は, burn-inの数(最初の実行は誤差が大きいため打ち切る必要があるから), 二つ目は何個のデータセットを用いるかを表している. 本当は, 100~1000回くらいが必要だが, ここでは10回ですぐに検証が終わるようにしている.

imp = mice.MICEData(df)
fml = "res ~ GPA + TUCE + PSI"
mice_ = mice.MICE(fml, sm.OLS, imp)
results = mice_.fit(10,10)
print(results.summary())

＊サンプルに載っているコードでは, moduleのimportとしてのmiceがあるにも関わらず, 上からclassとしてのmiceを入れているから, 一回コードをチャンクレベルで回すとエラーが起きる. 気をつけて…

おまけ

statsmodels は, 欠損値の補完がどのように行われているのか, いくつかmethodを用意している. fig のサイズの調整は, 線形の重回帰分析をstatsmodelsのグラフ機能を出来るだけ使って行うで書いた方法で対処すればいける.

imp.plot_bivariate("TUCE","GPA")
plt.show()

imp.plot_fit_obs("TUCE")
plt.show()

imp.plot_imputed_hist("GPA")
plt.show()

———-雑感(｀・ω・´)———-
MICEをpythonでやろうとしたけれども現状では, 自分でコードを書かないと厳しそう.
pythonもRも無料で統計解析のライブラリ充実してきているけれども, 典型的な構造化されたデータに関する解析なら有料のsasとかspss使った方が遥かに簡単に出来るし, 欲しい統計量全部出してくれる…
sasを真面目に勉強しようと誓った今回の記事でした…

参考文献

[1]Multiple Imputation with Chained Equations, https://www.statsmodels.org/stable/imputation.html
[2] statsmodels.formula.api tokenize issue in Python 3.7.1 #5343, https://github.com/statsmodels/statsmodels/issues/5343