多重検定問題を回避した変数選択（ホルム法｜Holm法）

前回は、多重検定問題を回避する方法の一つとして知られるBonferroni法（ボンフェローニ法）についてまとめたが、ついでに今回はHolm法（ホルム法）についてまとめた。

多重の多重性について

検定の多重性とは、仮設検定を何度も繰り返すことで、発生する問題のことをさす。

例えば、回帰や分類問題の解析において、無相関検定やt検定などを用いてp値が有意水準を下回った場合にその変数を説明変数と使用する場合を考える。
ある１つの変数を対象として検定をする場合には問題ないが、候補となる多くの変数に対して検定を実施する場合では、本来有意ではない（寄与や相関のない）変数のp値がたまたま有意水準を下回りることで有意であるみなされてしまう。

この様に、複数回繰り返された検定全体において帰無仮説が棄却される可能性は、familywise error rateと呼ばれる。 例えば、有意水準 $\alpha = 0.05$ の検定を20回繰り返すと、familywise error rate（1回でも帰無仮説が棄却される可能性）は 0.642となる。

多重検定を回避する(Familly wise Erro rate）を調整する方法。

Familywise error rateを調整する方法は大きく二つに分けられる。

1. F統計量やt統計量等の統計量に基づいた方法
2. 統計量から算出されたp値のみを操作する方法

1 の方法としては、F統計量を用いた(Fisher's LSD)法、ｔ統計量を用いたTukey's HSD法、Dunnet法などが知られている。

2 の方法としては、ボンフェローニ（Bonferroni）法やホルム（Holm）法が提唱されている。これらの方法は、統計量に依存せず、どのような検定に対しても利用できる汎用性が高い方法である。

今回はHolm法を使って回帰分析における変数選択を検証する。

Holm法について

Holm法は非常に保守的なBonferroni法の有意水準の補正方法を、改善したもの。
以下の手順で、p値の大きさに従って有意水準 $\alpha$ を設定する。

1. N個の帰無仮説を、p値の小さい順に並べる。
2. 最もp値が小さい第1順位の帰無仮説の有意水準を $ \frac{\alpha}{N} $ に設定する。
   p値 < $ \frac{\alpha}{N} $ の場合：第1順位の帰無仮説を棄却（対立仮説を採択し）、3に進む。
   p値 >$ \frac{\alpha}{N} $ の場合：第1順位以下のすべての帰無仮説の判定を保留する。
3. （第1順位の帰無仮説が棄却された場合）第2順位の帰無仮説の有意水準を $ \frac{\alpha}{N+1} $ に設定する。
   p値 < $ \frac{\alpha}{N+1} $ の場合：第2順位の帰無仮説を棄却（対立仮説を採択し）、4に進む。
   p値 > $ \frac{\alpha}{N+1} $ の場合：第2順位以下のすべての帰無仮説の判定を保留する。
4. 上記を繰り返す。

つまり、p値が小さい順位検定を行い、棄却できた場合には、有意水準をどんどん小さく（$\alpha$ を$N$, $N+1$, ...で割る）していき、帰無仮説が棄却されにくくしていく。

コード

Bostonの家賃データセットを使った家賃予測において、無相関検定で説明変数を選択する。
元の有意水準は0.05に設定し、変数選択したのちPLS回帰（部分的最小二乗回帰）で予測する。

準備

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_boston
from scipy.stats import pearsonr
from scipy.stats import spearmanr

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])

#サンプル数を絞る
df_X = df_X.iloc[::10,:]
df_y = df_y.iloc[::10,:]
print('変数の数:', df_X.shape[1])
#>>> 変数の数: 13

変数選択

目的変数のRENTとその他の変数で無相関検定を実施し、得られたp値がHolm法で補正した有意水準、および補正なしの有意水準を下回った場合にその変数を使用した。

#全ての変数と、目的変数「RENT」で無相関検定（ピアソン）を実施する。
p_values = []
for feature in df_X.columns:
rent = df_y['RENT'].values
feature = df_X[feature].values

#ピアソンの相関係数とp値を計算する。
result = pearsonr(rent, feature)
p_values.append(result[1])
p_values = np.array(p_values)

#Holm方による変数選択=======
#変数名とp値の辞書を作る
feature_dict = dict(zip(df_X.columns, p_values))

#p値の小さい順位に並び替える
feature_dict = dict(sorted(feature_dict.items(), key=lambda X:X[1]))
print('それぞれの変数のp値（昇順）')
for f, p in zip(feature_dict.keys(), feature_dict.values()):
print('{}:{}'.format(f,round(p,5)))
"""output
それぞれの変数のp値（昇順）
LSTAT:0.0
RM:0.0
INDUS:0.00078
PTRATIO:0.00283
ZN:0.00391
NOX:0.00439
TAX:0.00542
AGE:0.00591
CRIM:0.01602
CHAS:0.02801
B:0.04196
RAD:0.05783
DIS:0.08387
"""

a = 0.05
N = 1
selected_features_holm = []
for i, f in enumerate(feature_dict.keys()):
    a_holm = a/N

#棄却できれば変数をリストに追加して、有意水準を下げる（Nに１をたす）
if feature_dict[f] < a_holm:
    selected_features_holm.append(f)
    N += 1
else:
    break

#Holm方による変数選択ここまで=======
#補正をせずに変数選択（比較用）
selected_features = df_X.columns[np.where(p_values < a)[0]]
print('補正なしのαで選ばれた特徴量:',selected_features)
print('Holm法で選ばれた特徴量:',selected_features_holm)
"""output
補正なしのαで選ばれた特徴量: Index(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'TAX', 'PTRATIO',
       'B', 'LSTAT'],
      dtype='object')
Holm法で選ばれた特徴量: ['LSTAT', 'RM', 'INDUS', 'PTRATIO', 'ZN', 'NOX', 'TAX', 'AGE']
"""

補正なしαでは11変数、Holm法使用時は8変数が選ばれた。

予測と精度検証

定石ではないが、学習/訓練データを分割せずに検定を実施したため、クロスバリデーションで精度を検証した。
plsの成分数は２で固定し、変数選択しない場合と、変数選択した場合（補正ありとなし）とでそれぞれ精度を比較した。

from sklearn.model_selection import cross_val_predict
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error

#選ばれた変数セットのデータを作成
df_X_sel = df_X.loc[:,selected_features]
df_X_holm = df_X.loc[:,selected_features_holm]

#標準化
scaler_X = StandardScaler()
scaler_X_sel = StandardScaler()
scaler_X_holm = StandardScaler()
scaler_y = StandardScaler()
scaled_X = scaler_X.fit_transform(df_X)
scaled_X_sel = scaler_X.fit_transform(df_X_sel)
scaled_X_holm = scaler_X.fit_transform(df_X_holm)
scaled_y = scaler_y.fit_transform(df_y)

#予測モデル作成
pls = PLSRegression(n_components=2)

#全変数使用で予測
scaled_pred_y = cross_val_predict(estimator=pls,
                                X=scaled_X,
                                y=scaled_y,
                                cv=4)

#通常の検定（補正なしα）で選択した変数
scaled_pred_y_sel = cross_val_predict(pls,
                                    scaled_X_sel,
                                    scaled_y,
                                    cv=4)

#補正α使用して選択した変数
scaled_pred_y_holm = cross_val_predict(pls,
                                    scaled_X_holm,
                                    scaled_y,
                                    cv=4)

#標準化された値を元に戻す。
pred_y = scaler_y.inverse_transform(scaled_pred_y)
pred_y_sel = scaler_y.inverse_transform(scaled_pred_y_sel)
pred_y_holm = scaler_y.inverse_transform(scaled_pred_y_holm)

#精度を計算する
r2 = r2_score(df_y.values, pred_y)
r2_sel = r2_score(df_y.values, pred_y_sel)
r2_holm = r2_score(df_y.values, pred_y_holm)
rmse = np.sqrt(mean_squared_error(df_y.values, pred_y))
rmse_sel = np.sqrt(mean_squared_error(df_y.values, pred_y_sel))
rmse_holm = np.sqrt(mean_squared_error(df_y.values, pred_y_holm))
print('r2(変数選択なし):',round(r2,3))
print('r2(補正なしα):',round(r2_sel,3))
print('r2(Holm法):',round(r2_holm,3))
print('RMSE(変数選択なし):',round(rmse,3))
print('RMSE(補正なしα):',round(rmse_sel,3))
print('RMSE(Holm法):',round(rmse_holm,3))
"""output
r2(変数選択なし): 0.417
r2(補正なしα): 0.402
r2(Holm法): 0.618
RMSE(変数選択なし): 6.685
RMSE(補正なしα): 6.768
RMSE(Holm法): 5.412
"""

精度評価指標として決定係数（r²）と二乗平均平方根誤差（RMSE）を算出した。
結果を表にまとめると以下のようになる。

学習・訓練データの分割をせずに全データを使って検定を実施という条件の下ではあるが、r2, RMSEともに変数選択した方が良好な結果となった。
ちなみに、前回記事で検証したBonferroni検定と比べると、今回の検討に関してはBonferroni法の方が良好な結果となった。

参考

• メモ）分析者のためのデータ解釈学入門第二部
• Clinical Journal Club 1.多重比較

関連書籍