【python】sklearnのIterativeImputerで欠損値補完

　注意：IterativeImputerは本記事の執筆時点（2019年11月）で実験的な実装とされており、最新の仕様等はこの記事の内容と異なる可能性があります。常にstable版の公式のドキュメントを確認してください。

　公式のドキュメント
sklearn.impute.IterativeImputer — scikit-learn 0.21.3 documentation

はじめに

　説明変数の欠損値補完には様々な方法があり、いろいろな研究がなされています。

　直感的に思いつくのは、欠損部分は他の変数から回帰すれば良いという発想です。それを発展させて高級にすると、多重代入法というような手法になったりするらしいのですが、詳しくないので説明できません。どうせよく知らなくても機械学習モデルは構築できます。

　scikit-learnではIterativeImputerなるモデルを「開発中」です。ここで開発中と書いたのは、公式ですら

Note This estimator is still experimental for now: the predictions and the API might change without any deprecation cycle.

　という扱いになっているからです。「いつの間にか仕様が変わったり、最悪消えたりしてても文句言うなよ！」といったところでしょうか。それでも一応masterブランチにマージされていて、普通にsklearnをインストールすれば（2019年11月時点では）使えます。

　ということで、便利そうなので紹介します。ただし、そこそこ癖があるので気をつけます。

使い方

　まずリファレンスにも書いてあることですが、こいつを使うためには余計なものもimportする必要があります。

from sklearn.experimental import enable_iterative_imputer  # 本当におまじない
from sklearn.impute import IterativeImputer  # 使うのこっち

　使っているlintによっては「imported but unused」とか言われたりしますが、無視してください。

　あとはハイパーパラメータをどうするかですが、あんまり凝ったことはしないのをおすすめします。と言いつつ、何も設定しないで使うのもおすすめしません。

　重要なパラメータを重要な順番に並べてみます。

estimator=None,

　内部で回帰を行って補完するので、回帰モデルを渡せるという訳です。デフォルトはベイジアンリッジ（sklearn.linear_model.BayesianRidge）になります。なんでデフォルトがこれなのかは後述します。

sample_posterior=False

　そのままだと単一代入法、Trueにすると多重代入法になるのだと思いますが、いかんせん私が知識不足で、またドキュメントの記述がそっけないので細かいアルゴリズムまではよくわかりません。多重代入法の方がかっこよくて性能も良さそうに思えますが、実はその場合は

Estimator must support return_std in its predict method if set to True.

　という制約があり、これは割ときつい制約です（対応しているモデルはほとんどありません）。デフォルトのBayesianRidgeはベイズ法なので予測の標準偏差が出せるということから選定されていると思われます。

max_iter=10
tol=0.001

　特に説明は要らないと思います。なんか上手く行かないときはこの辺をいじる（増やしてみる）。

　ということで、議論の余地があるのは多重代入法にするのか単一代入法にするのかです。多重代入法の場合、現実的なestimatorの選択肢はたぶんベイジアンリッジだけで、これは「ベイズでリッジだからつよつよ」と見せかけて基本的にただの線形モデルです（回帰式が一次式のやつ）。機械学習に突っ込むようなデータだとこれは微妙なことが多そうなので、多重代入法を諦めて単一代入法にすることで自由に好きな非線形回帰手法を使うということにした方が良いケースが多いと思います。

　あとは何で回帰するか？　ですが、ここの回帰モデルのパラメータチューニングのことなんて考えたくもないので、適当な設定でもそこそこ動く奴を使います。と言って真っ先に思いつくのはランダムフォレストですね。

　以上をまとめると、だいたいこんな感じで定義することになると思います。

from sklearn.ensemble import RandomForestRegressor
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

rf = RandomForestRegressor(n_estimators=500, n_jobs=-1)
imp = IterativeImputer(estimator=rf)

　他のパラメータは気になったら適当に弄ってください。大抵は弄らなくても大勢に影響しません。あと、ランダムフォレストはそこそこ遅いので、気になったら木を浅くして本数を減らすことで高速向きのチューニングにするか、他のものも試してみてください。
　

実験

　SimpleImputerのときと同じことを試します。

【python】sklearnのSimpleImputerで欠損値補完をしてみる - 静かなる名辞

　irisを欠損させて補完して予測精度を見ます。

import numpy as np
from sklearn.datasets import load_iris
from sklearn.neighbors import KNeighborsClassifier    
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

# 再現性確保
np.random.seed(0)

# データ生成
dataset = load_iris()
X_train, X_test, y_train, y_test = train_test_split(
    dataset.data, dataset.target, stratify=dataset.target, random_state=0)
mask = ~np.random.randint(10, size=X_train.shape).astype(np.bool)
X_train[mask] = np.nan

# モデル定義
knn = KNeighborsClassifier()

# 処理
def drop():
    print("# drop")
    idx = ~(np.isnan(X_train).any(axis=1))
    X_train_, y_train_ = X_train[idx], y_train[idx]
    knn.fit(X_train_, y_train_)
    y_pred = knn.predict(X_test)
    print(classification_report(
        y_test, y_pred, digits=3, 
        target_names=dataset.target_names))
                    
def impute():
    print("# impute")

    rf = RandomForestRegressor(
        n_estimators=200, max_depth=3, n_jobs=-1)
    imp = IterativeImputer(estimator=rf, max_iter=20)
    pl = Pipeline([("imputer", imp), ("KNN", knn)])
    pl.fit(X_train, y_train)
    y_pred = pl.predict(X_test)
    print(classification_report(
        y_test, y_pred, digits=3, 
        target_names=dataset.target_names))
    
drop()
impute()

# drop
              precision    recall  f1-score   support

      setosa      1.000     1.000     1.000        13
  versicolor      0.929     1.000     0.963        13
   virginica      1.000     0.917     0.957        12

    accuracy                          0.974        38
   macro avg      0.976     0.972     0.973        38
weighted avg      0.976     0.974     0.974        38

# impute
              precision    recall  f1-score   support

      setosa      1.000     1.000     1.000        13
  versicolor      1.000     1.000     1.000        13
   virginica      1.000     1.000     1.000        12

    accuracy                          1.000        38
   macro avg      1.000     1.000     1.000        38
weighted avg      1.000     1.000     1.000        38

　単に欠損値のある行を学習データから落とすより良いわけです。まあ、irisのときは平均値補完でやっても同じ結果になりましたが。

　データセットをwineにしてみます。

# drop
              precision    recall  f1-score   support

     class_0      0.812     0.867     0.839        15
     class_1      0.750     0.667     0.706        18
     class_2      0.538     0.583     0.560        12

    accuracy                          0.711        45
   macro avg      0.700     0.706     0.702        45
weighted avg      0.714     0.711     0.711        45

# impute
              precision    recall  f1-score   support

     class_0      0.800     0.800     0.800        15
     class_1      0.733     0.611     0.667        18
     class_2      0.533     0.667     0.593        12

    accuracy                          0.689        45
   macro avg      0.689     0.693     0.686        45
weighted avg      0.702     0.689     0.691        45

　dropより良くなるかと思いましたが、そうでもないようです。平均値で補完するよりは良いみたいですが……。

考察

　そもそも欠損値補完というアプローチは、

欠損値が多すぎるとそもそもまともに予測できないし、予測の悪さが大勢に影響を及ぼして全体のパフォーマンスを悪化させかねない
欠損値が少ないならdropしても大勢に影響はない

　という矛盾を抱えている訳ですが、そこそこの欠損があるときにdropするよりちょっと良くなるかな？　という可能性を追求するためにあるものだと思います。

まとめ

　使うの面倒くさいしこの先どうなるかは読めないけど、使えるみたいです。ただし、できるからといってパフォーマンスの上で凄い効能があるとか、そういうものではないみたいです。