統計とイモリとバーベルと

勉強の備忘録メイン

なぜFGOを始めてしまったのか

Missing Data: ① ありがちな対処法とその影響

個人的に,欠測値の処理はデータ分析にあたって考えたくないことの筆頭です.
とりあえず早く結果を出したいから暫定の処理をして分析を進めるも,「欠測値をなんとかしなければ」というのは常に頭の片隅にこびりついている.それをストレスに感じつつも,目の前のタスクに追われていたり,モデル選択・特徴量設計に取り組む魅力のほうが大きかったりと,つい後回しにしてしまう.そんなこんなで,気づいた頃には分析も佳境に......
治験関係の統計解析では,欠測値の処理も含めてすべての条件を本番データが揃う前に決定してしまうので,ここまでわたわたすることもないのですが,それにしても,欠測処理が退屈で考えたくないことだ,という印象に変わりはありません.

欠測値について考えたくない理由5選
1. 理論はともかく,実データの欠測処理は面白くない(賛否ありそうですが…)
2. 早く何らかの分析結果を出したい
3. うまくやったところで評価されない
4. どのような害がどの程度起こるのか想像しづらい
5. いざ実データに相対すると,何をするのが正しいのかわからない

3は諦めるしかないでしょうね.前処理全般に言えることなのでしょうが…. 今回はシリーズで欠測処理についてまとめる予定です.可能な限り手早く,頭を使わずに欠測値処理を行えるように勉強を進めます.

シリーズを通して,こちらの書籍を参考にしていきます.調査観察データの統計科学で有名な星野崇宏教授も編集・執筆に加わっておられるようです.

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際 第1巻)

今日のお題は「① ありがちな対処法とその影響」. 「とりあえずデータがあるものだけで分析した」「とりあえず平均値を代入しておいた」といった方針がどの程度結果を歪ませるのかを見ていきます.

サンプルデータの生成

今回用いるサンプルデータは,以下のモデルに従って生成しました.


 y = 5 + 2x_1 - 3x_2 + \varepsilon \\\
 x_1, x_2, \varepsilon \sim N(0, 16)

上記と同じ線形モデルでOLSを行って,回帰係数 \betaを推定する状況を想定ます.欠測のパターンには,現実の起こりそうな3つのケースを用意しました.
1. 各サンプルの x1が欠測になるかどうかが,ランダムに決まるケース(independent)
2. 各サンプルの x1が欠測になるかどうかが, x2の値に依存するケース(x2-dependent)
3. 各サンプルの x1が欠測になるかどうかが, yの値に依存するケース(y-dependent)
上記以外に,欠測に影響する変数は存在しないものとします.

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.stats import norm

missing_ratio = 0.3
sample_size = 200
seed = 123

if __name__ == "__main__":
    np.random.seed(seed)
    x1 = np.random.randn(sample_size)
    x2 = np.random.randn(sample_size)
    y = 5 + 3*x1 - 2*x2 + np.random.normal(0, 4, sample_size)

    df = pd.DataFrame({"y": y,
                       "x1": x1,
                       "x2": x2})

    df['x1_case1'] = df.x1.map(
        lambda x: 1 if np.random.rand() >= missing_ratio else 0
    )
    df['x1_case2'] = (df.x2 >= norm.ppf(missing_ratio)).map(
        lambda x: np.int(x)
    )
    df['x1_case3'] = (df.y >= df.y.quantile(missing_ratio)).map(
        lambda x: np.int(x)
    )
    print(df.head(10))

    fig, ((ax11, ax12), (ax21, ax22)) = plt.subplots(nrows=2,
                                                     ncols=2,
                                                     figsize=(9,9),
                                                     dpi=80)
    fig.subplots_adjust(wspace=0.3, hspace=0.3)

    xvar = df.x1
    yvar = df.y
    ax11.plot(xvar, yvar, 'o',
              color='navy',
              markerfacecolor='w',
              alpha=0.7)
    ax11.set_title('complete data')
    ax11.set_xlabel('x1')
    ax11.set_ylabel('y')
    ax11.set_xlim(-3, 3)
    ax11.set_ylim(-10, 20)

    xvar = df.x1.loc[df.x1_case1 == 1]
    yvar = df.y.loc[df.x1_case1 == 1]
    ax12.plot(xvar, yvar, 'o',
              color='navy',
              markerfacecolor='w',
              alpha=0.7)
    ax12.set_title('independent missing')
    ax12.set_xlabel('x1')
    ax12.set_ylabel('y')
    ax12.set_xlim(-3, 3)
    ax12.set_ylim(-10, 20)

    xvar = df.x1.loc[df.x1_case2 == 1]
    yvar = df.x2.loc[df.x1_case2 == 1]
    ax21.plot(xvar, yvar, 'o',
              color='darkgreen',
              markerfacecolor='w',
              alpha=0.7)
    ax21.set_title('x2-dependent missing')
    ax21.set_xlabel('x1')
    ax21.set_ylabel('x2', color="green")
    ax21.set_xlim(-3, 3)
    ax21.set_ylim(-3, 3)
    ax21.hlines([norm.ppf(missing_ratio)], -3, 3, "red", linestyles='dashed')

    xvar = df.x1.loc[df.x1_case3 == 1]
    yvar = df.y.loc[df.x1_case3 == 1]
    ax22.plot(xvar, yvar, 'o',
              color='navy',
              markerfacecolor='w',
              alpha=0.7)
    ax22.set_title('y-dependent missing')
    ax22.set_xlabel('x1')
    ax22.set_ylabel('y')
    ax22.set_xlim(-3, 3)
    ax22.set_ylim(-10, 20)
    ax22.hlines([df.y.quantile(missing_ratio)], -10, 20, "red", linestyles='dashed')

    plt.savefig('../fig/sample_scatter.png', format='png')
    df.to_pickle('../data/sample_data.pkl')

散布図を描くと,次のようになります.x2-dependentおよびy-dependentは,それぞれ依存する値が一定以下の場合,x1が欠測となります. f:id:kaizu_3324:20180113221719p:plain

まずは「欠測のないレコードだけを使う」という発想から

「欠測のないレコード」のことを「完全ケース(complete case)」と呼びます. 用意した各サンプルデータに対して,重回帰分析を実行します.分析にはstatsmodelsを使用しました.

import pandas as pd
import statsmodels.api as sm

if __name__ == "__main__":
    df = pd.read_pickle('../data/sample_data.pkl')

    f = open('../result/regress_listwise.txt', 'w')

    print("----------complete data----------", file=f)
    X = df.loc[:,['x1','x2']].values
    y = df.y.values.reshape(-1, 1)
    model = sm.OLS(y, sm.add_constant(X))
    result = model.fit()
    print(result.summary(), file=f)

    print("independent-----------------------------------------------", file=f)
    d1 = df.loc[df.x1_case1 == 1]
    X = d1.loc[:,['x1','x2']].values
    y = d1.y.values.reshape(-1, 1)
    model = sm.OLS(y, sm.add_constant(X))
    result = model.fit()
    print(result.summary(), file=f)

    print("x2-dependent----------------------------------------------", file=f)
    d1 = df.loc[df.x1_case2 == 1]
    X = d1.loc[:,['x1','x2']].values
    y = d1.y.values.reshape(-1, 1)
    model = sm.OLS(y, sm.add_constant(X))
    result = model.fit()
    print(result.summary(), file=f)

    print("y-dependent-----------------------------------------------", file=f)
    d1 = df.loc[df.x1_case3 == 1]
    X = d1.loc[:, ['x1', 'x2']].values
    y = d1.y.values.reshape(-1, 1)
    model = sm.OLS(y, sm.add_constant(X))
    result = model.fit()
    print(result.summary(), file=f)
Result

それでは結果を見ていきましょう.まずはComplete dataから.coefを見ると,おおよそ真のモデルと同じ値をとっていることがわかります.

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.437
Model:                            OLS   Adj. R-squared:                  0.432
Method:                 Least Squares   F-statistic:                     76.55
Date:                Sat, 13 Jan 2018   Prob (F-statistic):           2.53e-25
Time:                        22:06:50   Log-Likelihood:                -564.20
No. Observations:                 200   AIC:                             1134.
Df Residuals:                     197   BIC:                             1144.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2866      0.292     18.122      0.000       4.711       5.862
x1             2.9879      0.275     10.861      0.000       2.445       3.530
x2            -1.8755      0.313     -5.998      0.000      -2.492      -1.259
==============================================================================
Omnibus:                        2.495   Durbin-Watson:                   1.844
Prob(Omnibus):                  0.287   Jarque-Bera (JB):                2.560
Skew:                          -0.259   Prob(JB):                        0.278
Kurtosis:                       2.800   Cond. No.                         1.17
==============================================================================

次にindependent missing.こちらは標準誤差がやや大きくなったものの,coefを見ると,点推定値は真のモデルに近い値を示しています.欠測値を削除することの影響は,サンプル数が減少したことと同じことが起きたと考えてよいでしょう.

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.476
Model:                            OLS   Adj. R-squared:                  0.468
Method:                 Least Squares   F-statistic:                     63.48
Date:                Sat, 13 Jan 2018   Prob (F-statistic):           2.39e-20
Time:                        22:06:50   Log-Likelihood:                -406.23
No. Observations:                 143   AIC:                             818.5
Df Residuals:                     140   BIC:                             827.4
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.6845      0.352     16.149      0.000       4.989       6.380
x1             3.1826      0.323      9.856      0.000       2.544       3.821
x2            -2.0335      0.373     -5.457      0.000      -2.770      -1.297
==============================================================================
Omnibus:                        2.617   Durbin-Watson:                   1.816
Prob(Omnibus):                  0.270   Jarque-Bera (JB):                2.649
Skew:                          -0.317   Prob(JB):                        0.266
Kurtosis:                       2.796   Cond. No.                         1.19
==============================================================================

次にx2-dependentのケースです.こちらもindependentのケース同様,推定値は正しい値が求められているように見えます.R2乗値もcomplete dataと遜色ありません.

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.433
Model:                            OLS   Adj. R-squared:                  0.424
Method:                 Least Squares   F-statistic:                     50.36
Date:                Sat, 13 Jan 2018   Prob (F-statistic):           5.59e-17
Time:                        22:06:50   Log-Likelihood:                -375.32
No. Observations:                 135   AIC:                             756.6
Df Residuals:                     132   BIC:                             765.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2172      0.391     13.333      0.000       4.443       5.991
x1             2.9913      0.314      9.512      0.000       2.369       3.613
x2            -1.7149      0.513     -3.343      0.001      -2.730      -0.700
==============================================================================
Omnibus:                        2.270   Durbin-Watson:                   1.969
Prob(Omnibus):                  0.321   Jarque-Bera (JB):                2.235
Skew:                          -0.308   Prob(JB):                        0.327
Kurtosis:                       2.870   Cond. No.                         1.85
==============================================================================

最後に,y-dependentのケースを見ていきましょう.先程までと大きく異なり,回帰係数の推定値が真のモデルから大きく外れています.R2乗値も低く,今回のデータを説明するモデルとしては機能していないと考えられます.

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.244
Model:                            OLS   Adj. R-squared:                  0.233
Method:                 Least Squares   F-statistic:                     22.11
Date:                Sat, 13 Jan 2018   Prob (F-statistic):           4.75e-09
Time:                        22:06:50   Log-Likelihood:                -361.03
No. Observations:                 140   AIC:                             728.1
Df Residuals:                     137   BIC:                             736.9
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          7.3014      0.310     23.558      0.000       6.689       7.914
x1             1.9458      0.311      6.251      0.000       1.330       2.561
x2            -1.0487      0.309     -3.389      0.001      -1.660      -0.437
==============================================================================
Omnibus:                        2.202   Durbin-Watson:                   2.146
Prob(Omnibus):                  0.332   Jarque-Bera (JB):                2.132
Skew:                           0.299   Prob(JB):                        0.344
Kurtosis:                       2.906   Cond. No.                         1.69
==============================================================================

後々の回で理論的にも詳しく見ていく予定ですが,実は,シンプルな線形モデルを想定する場合,欠測で気をつけなければならないのは欠測するか否かが目的変数に依存している場合のみ,完全ケース分析の推定値が最尤法による一致性を失うことが知られています(今回はあくまでOLSによるデモンストレーションですが…).
このあたり,経済学方面から統計学を学んだ方などは特に腹落ちしやすい部分ではないかと思います.「欠測するかどうか」という変数rをモデルに組み込んだ場合に,rが内生変数となるケースでは推定量が一致性を失う,ということですね.

次に,「代表値で欠測を埋めたらいいんじゃないの?」という発想について

欠測した部分を,その変数の代表値(平均値,メディアン等)で埋めるという発想は自然なことだと思いますし, 感覚的に,これをやりたくなる気持ちは理解できます.しかも,市販の統計ソフトには(極めて残念ながら)この機能が標準搭載されていることもあるようです.しかし,このあと見ていくように,推定値の分散を過小評価するという重大な欠点があります.

independent データでの実例

まずは,平均値代入を行った場合のデータと,真の分布との間のヒストグラムを比較してみましょう.データは,先ほどのx1 (independent missing)を用いました.
グレーが真のデータのヒストグラム,赤が平均値代入を行ったデータのヒストグラムです. 当然ではありますが,平均値の部分が突出しているのがわかります.

f:id:kaizu_3324:20180113232720p:plain

次に,これらのデータから簿平均を推定し,t検定に基づく信頼区間を算出します. 完全データによる分析結果だけでなく,欠測のない元のデータよりも信頼区間が狭くなっており,分散を過小評価していることがわかります.

# x1 (independent missing) 完全データ分析
Mean: -0.0165
95%CI: (-2.3590149256462909, 2.3260990437590743)

# x1 (independent missing) 平均値代入
Mean: -0.0165
95%CI: (-1.6839755909234084, 1.6510597090361918)

# x1 欠測なし
Mean: 0.0038
95%CI: (-2.1910807760042403, 2.1986542875844539)

ソースコード

import numpy as np
import pandas as pd
from scipy.stats import t, sem
from matplotlib import pyplot as plt

if __name__ == "__main__":
    df = pd.read_pickle('../data/sample_data.pkl')

    x1_mean = df.x1.loc[df.x1_case1 == 1].mean()

    df["x1_imputed"] = df.x1*df.x1_case1 + x1_mean*(1-df.x1_case1)
    df.x1.hist(bins=25, color="lightgray")
    df.x1_imputed.hist(bins=25, color="salmon", alpha=0.5)
    plt.legend()
    plt.savefig("../fig/histogram_mean_imputation.png", format="png")

    arr = df.x1.loc[df.x1_case1 == 1].values
    m = np.mean(arr)
    e = np.var(arr, ddof=1)
    print("Mean: {0:.4f}".format(m))
    print("95%CI:", t.interval(alpha=0.95, df=len(arr), loc=m, scale=e))

    arr = df.x1_imputed.values
    m = np.mean(arr)
    e = np.var(arr, ddof=1)
    print("Mean: {0:.4f}".format(m))
    print("95%CI:", t.interval(alpha=0.95, df=len(arr), loc=m, scale=e))

    arr = df.x1.values
    m = np.mean(arr)
    e = np.var(arr, ddof=1)
    print("Mean: {0:.4f}".format(m))
    print("95%CI:", t.interval(alpha=0.95, df=len(arr), loc=m, scale=e))

今回のまとめ

  • 今後,欠測値処理について勉強を進める予定
  • 線形モデルを当てはめる場合に注意すべきは,欠測状況が目的変数に依存するか否か
  • 平均値代入は分散を過小評価するので,すべきではない
GitHub

github.com