多重検定問題を回避した変数選択(Bonferroni法)

多重検定問題を回避する方法の一つとして知られるBonferroni法(ボンフェローニ法)についてまとめた。

多重の多重性について

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

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

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

多重検定問題を回避する(Famillywise Error rateを調整する)方法

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

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

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

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

今回はBonferroni法を使って回帰分析における変数選択を実施した。

Bonferroni法について

検定の全体を踏まえてp値を評価する方法。検定の回数をM回とすると、はじめに設定した有意水準 $\alpha$ を $M$ で割った $\alpha '$ を新たな有意水準として採用する。

$$ \alpha' = \frac{\alpha}{M} $$

Bonferroni法の問題点

検定の回数が多いと有意水準が極端に小さくなってしまい「本当は変数間の関係があるのに有意水準が小さすぎて検出できない」という問題がある)

このように第2種の過誤の可能性が高くなるため、「p値 > 補正α」となった場合は、帰無仮説は採択されるのではなく棄却が保留されると考えるのが妥当。

コード

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

#Bonferroni法で補正したの有意水準と補正なしの有意水準(比較のため)を設定する
a = 0.05
a_bon = 0.05 / (df_X.shape[1])

print('有意水準:',a)
print('補正有意水準:',round(a,5))
"""output
有意水準: 0.05
補正有意水準: 0.00385
"""

変数選択

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

#全ての変数と、目的変数「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)
print('それぞれの変数のp値')
for f, p in zip(df_X.columns, p_values):
print('{}:{}'.format(f,round(p,5)))
"""output
それぞれの変数のp値
CRIM:0.01602
ZN:0.00391
INDUS:0.00078
CHAS:0.02801
NOX:0.00439
RM:0.0
AGE:0.00591
DIS:0.08387
RAD:0.05783
TAX:0.00542
PTRATIO:0.00283
B:0.04196
LSTAT:0.0
"""

selected_features = df_X.columns[np.where(p_values < a)[0]]
selected_features_bon = df_X.columns[np.where(p_values < a_bon)[0]]
print('補正なしαで選ばれた特徴量:',selected_features)
print('補正αで選ばれた特徴量:',selected_features_bon)
"""output
補正なしαで選ばれた特徴量: Index(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'TAX', 'PTRATIO', 'B', 'LSTAT'],
      dtype='object')
補正αで選ばれた特徴量: Index(['INDUS', 'RM', 'PTRATIO', 'LSTAT'], dtype='object')
"""

補正なしαでは11変数、補正ありαでは4変数が選ばれた。

予測と精度確認

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

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_bon = df_X.loc[:,selected_features_bon]

#標準化
scaler_X = StandardScaler()
scaler_X_sel = StandardScaler()
scaler_X_bon = StandardScaler()
scaler_y = StandardScaler()
scaled_X = scaler_X.fit_transform(df_X)
scaled_X_sel = scaler_X.fit_transform(df_X_sel)
scaled_X_bon = scaler_X.fit_transform(df_X_bon)
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_bon = cross_val_predict(pls,
                                    scaled_X_bon,
                                    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_bon = scaler_y.inverse_transform(scaled_pred_y_bon)

#精度を計算する
r2 = r2_score(df_y.values, pred_y)
r2_sel = r2_score(df_y.values, pred_y_sel)
r2_bon = r2_score(df_y.values, pred_y_bon)
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_bon = np.sqrt(mean_squared_error(df_y.values, pred_y_bon))
print('r2(変数選択なし):',round(r2,3))
print('r2(補正なしα):',round(r2_sel,3))
print('r2(補正ありα):',round(r2_bon,3))
print('RMSE(変数選択なし):',round(rmse,3))
print('RMSE(補正なしα):',round(rmse_sel,3))
print('RMSE(補正ありα):',round(rmse_bon,3))
"""output
r2(変数選択なし): 0.417
r2(補正なしα): 0.402
r2(補正ありα): 0.678
RMSE(変数選択なし): 6.685
RMSE(補正なしα): 6.768
RMSE(補正ありα): 4.971
"""

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

 

r2RMSE
変数選択なし0.4176.685
補正なしαで変数選択0.4026.768
補正ありαで変数選択0.6784.971

 

学習・訓練データの分割をせずに全データを使って検定を実施という条件の下ではあるが、r2, RMSEともにBonferroni法で補正した有意水準 $\alpha$ を用いて変数選択した方が良好な結果となった

参考

関連書籍