multiple imputation について

Question

Pythonでのmultiple imputationについての質問です。
性別（男女）、年齢（0-105才）、身長（40－200㎝）、体重（3-200㎏)の4変数の10000サンプルあるデータフレーム(df1)に関して、性別以外にところどころNan（10%以下）があります。
これらをMultiple imputationを使って補完しようとしました。
プログラムは以下です。

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
imp = IterativeImputer(max_iter=10, sample_posterior=True)

X = df1.iloc[:, 0:4]
y = df1.iloc[:, 0:4]

imp.fit(X)
a = np.round(imp.transform(y))
feature = df1.columns[0:4]
df2 = pd.DataFrame(data=a ,columns= feature)

df2

このようにすると、Imputationの設定は
IterativeImputer(add_indicator=False, estimator=None,
imputation_order='ascending', initial_strategy='mean',
max_iter=10, max_value=None, min_value=None,
missing_values=nan, n_nearest_features=None, random_state=None,
sample_posterior=True, tol=0.001, verbose=0)
と出ます。
特にエラーなくdf2としてNanがimputeされたデータフレームが返されますが、問題は0才なのに身長160㎝、体重60㎏などと現実的でない数値がimputeされてしまっています。
確かに身長と体重の変数の平均はそれぞれ、そのくらいなのですが、
実際の想定範囲内の身長50-60㎝、体重6㎏などくらいにimputeされていてほしいのです。
（実際、データセットには0歳児の欠損してないデータも60サンプルくらい入っている）

Multiple imputationの原理のためなのか、私のプログラミングミスのせいなのか判断できなくて困っています。

Multiple imputationに詳しい方で、この問題点に関してご意見聞かせて頂けるとありがたいです。
さらにこのデータセットの場合補完法で適していると思われる手法があれば、ご意見くださるとありがたいです。

Comments

[hayataka2049]

ダミーとか似ているものでも構わないので、同じような状況を再現できるデータセットを提供していただけませんか？　検証してみたいのですが、データで苦労するので……

[SyunSyun]

早速のお返事ありがとうございます。
生データに少し手を加えた1000例のExcelファイルをGitHub上に保存しました。
URLは以下です。
https://github.com/yamanaka1234/DummyData/blob/master/年齢、体重、身長、性別.xlsx
GitHubを使用するのは初めてなので、これでhayataka2049様がDataを入手できるか不安ですが、入手不可能であったら、また教えて頂けますか？

[hayataka2049]

ファイルを確認できました。今日中くらいに試してみます。

[SyunSyun]

有難うございます。
もし、何かわかることがあったらまた教えてください。
個人的にはプログラミングの問題というより、Mmultiple imputationの原理のためと思うのですが。

[hayataka2049]

今回答書いてますが、まあその通りです。すぐに投げます。

Answer 1

まず、4つしか変数がないのに一行で2つの欠損があるようなデータは悪条件ですから、ある意味では多少の妥協は必要になります。

---

次に、max_iterを増やしたり、max_value, min_valueを指定したりすることで多少の改善が見込めます。そんなに大げさな改善にはならないと思いますが。

---

最後に、これが一番深刻な問題ですが、IterativeImputerが内部で用いている回帰モデルはベイジアンリッジで、線形モデルになります。性別・年齢・身長・体重の関係が線形ということは考えづらいので、予測に限界があります。

こういうとき、非線形回帰系のモデルを手軽に突っ込めると助かるのですが、IterativeImputerで使える回帰モデルには

If sample_posterior is True, the estimator must support return_std in its predict method.
sklearn.impute.IterativeImputer — scikit-learn 0.21.3 documentation

というきつい制約があり、条件を満たす非線形回帰手法のモデルは軽く探した範囲だとGaussianProcessRegressorくらいしかありません。やってみましたがチューニングがよくわからないし、すぐ警告を吐くのでちょっと難があります。ちゃんと使えばたぶんできます。

一例ですが、私が使ってなんとかまともにうごいたかなという気がした定義です。パラメータの妥当性は保証しません。

# 追加で必要になるimportのみ記載。以下のコードでも同じ
from sklearn.preprocessing import StandardScaler
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.pipeline import Pipeline

ss = StandardScaler()
gp = GaussianProcessRegressor(alpha=1e-5, normalize_y=True)
est = Pipeline([("scaler", ss), ("reg", gp)])
imp = IterativeImputer(est, max_iter=10, sample_posterior=True)

---

multiple imputationを諦めればSVRだろうがランダムフォレストだろうが使えるので、それはそれで楽ではあります。

from sklearn.ensemble import RandomForestRegressor

est = RandomForestRegressor(n_estimators=100, n_jobs=-1)
imp = IterativeImputer(est, max_iter=10, sample_posterior=False)

理論的な良し悪しはともかく、こっちの方が良い結果が得られるなら採用するというのもありだと思います。

---

他の方法は、非線形回帰ができるように然るべき変換をするTransformerを作ってやることです。

どう変換するかというと、選択肢の一つは多項式回帰にして無理やり解くことです。やるとたしかにできますが、しょせん多項式なのでちょっと当てはまりが悪いかもしれません。試しましたが、年齢で100歳とか出ちゃったりして、愉快な感じでした。

from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import BayesianRidge
from sklearn.pipeline import Pipeline

ss = StandardScaler()
pf = PolynomialFeatures(degree=3, include_bias=False)
br = BayesianRidge()
est = Pipeline([("scaler", ss), ("poly_f", pf), ("reg", br)])
imp = IterativeImputer(est, max_iter=10, sample_posterior=True)

もう少し正攻法で考えると、とりあえず身長体重は大雑把に言えば線形とみなしてよさそうで、性別は0か1ですからどうでもよく、加齢に対する成長の効果を適切に反映するような変換だけすればなんとかなる、という発想があり得るかもしれません。ということで、適当に作ったのがこれです。

from sklearn.preprocessing import FunctionTransformer
from sklearn.compose import ColumnTransformer
from sklearn.linear_model import BayesianRidge

def trans(x):
    """ここ見ながら適当に作りました↓
    https://www.graphsketch.com/
    """
    return 20-np.exp(-0.3*x+3)

ft = FunctionTransformer(trans)
ct = ColumnTransformer([("age", ft, [1])], remainder="passthrough")
br = BayesianRidge()
est = Pipeline([("scaler", ct), ("reg", br)])
imp = IterativeImputer(est, max_iter=10, sample_posterior=False)

だいたい20歳以上で成長が止まるという状況を想定した関数を作ってやってますが・・・だいたいそれっぽくはなったけど、もう少し工夫が必要な気もしました。

Comments

[SyunSyun]

hayataka2049様。大変丁寧な解説を有難うございました。
うまくいかない最大の理由は、リグレッサーが線形回帰になっていたことなんですね。
教えて頂いた事をもとに自分なりにさらに、予測性能があがるものに改良したいと思います。
色んな手段を教えて頂いて、感動しました。
これからも地道に頑張ります。