[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

コメント

タイトルとURLをコピーしました