AICによるモデル選択が信頼区間の被覆確率に与える影響

AICなどの基準を用いた変数選択を行うことがあると思います．
しかし，この処理は特に信頼区間の被覆確率に影響を与えます．
実験を通して説明します．

実験

設定

以下のようなデータ生成過程に対し実験をします．
$$
y_i = \beta_0 + x_{i1}\beta_1 + x_{i2}\beta_2 + \varepsilon_i, \varepsilon_i \sim \mathcal{N}(0, 1), i = 1, \ldots, N
$$
ここで，回帰係数の真値は$\beta_0=1, \beta_1=1, \beta_2=0$とします．すなわち，$x_{i2}$は不要な説明変数です．
また，説明変数は

\left[\begin{array}{c}
    x_{i1} \\
    x_{i2}
\end{array}\right] \sim \mathcal{N}\left(\left[\begin{array}{c}
    0 \\
    0
\end{array}\right], \left[\begin{array}{cc}
    1.0 & 0.8 \\
    0.8 & 1.0
\end{array}\right]\right)
, i = 1, \ldots, N

と発生させます．

このようなデータ生成過程に対し，以下の3つの手法を比較します．

• 手法1：全ての説明変数を利用した線形回帰
• 手法2：$x_{i1}$のみ利用した線形回帰
• 手法3：全ての説明変数を利用したモデルと$x_{i1}$のみ利用したのみ利用したモデルをAICにより選択

これの3つの手法に対して，$\beta_1$について$95$%信頼区間に真の値が含まれている割合を評価します．

実験の繰り返し回数は$100000$回，$N$は$100$とします．

実験の結果

手法1の実験結果

まず，全ての説明変数を利用した線形回帰について実験します．
実験コードは以下です．

import numpy as np
from scipy.stats import t as t_dist

np.random.seed(1234)

cov = np.array([
    [1.0, 0.8],
    [0.8, 1.0]
])
N = 100
num_iter = 100000
beta = np.array([1.0, 1.0, 0.0])

prob_in_ci = 0
for _ in range(num_iter):
    X = np.concatenate([np.ones((N, 1)), np.random.multivariate_normal(mean=np.zeros(2), cov=cov, size=N)], axis=1)
    y = X @ beta + np.random.randn(N)

    model = sm.OLS(y, X)
    results = model.fit()

    prob_in_ci += float((results.conf_int()[1][0] < beta[1]) & (results.conf_int()[1][1] > beta[1])) / num_iter

print('信頼区間に真の値が含まれる割合： ', prob_in_ci)

信頼区間に真の値が含まれる割合：  0.9499299999983116

手法2の実験結果

次に，$x_{i1}$のみ利用した線形回帰について実験します．
実験コードは以下です．

prob_in_ci = 0
for _ in range(num_iter):
    X = np.concatenate([np.ones((N, 1)), np.random.multivariate_normal(mean=np.zeros(2), cov=cov, size=N)], axis=1)
    y = X @ beta + np.random.randn(N)

    model = sm.OLS(y, X[:, 0:2])
    results = model.fit()

    prob_in_ci += float((results.conf_int()[1][0] < beta[1]) & (results.conf_int()[1][1] > beta[1])) / num_iter

print('信頼区間に真の値が含まれる割合： ', prob_in_ci)

信頼区間に真の値が含まれる割合：  0.9506399999983084

手法3の実験結果

最後に，全ての説明変数を利用したモデルと$x_{i1}$のみ利用したのみ利用したモデルをAICにより選択する手法について実験します．
実験コードは以下です．

prob_in_ci = 0

for _ in range(num_iter):
    X = np.concatenate([np.ones((N, 1)), np.random.multivariate_normal(mean=np.zeros(2), cov=cov, size=N)], axis=1)
    y = X @ beta + np.random.randn(N)

    model1 = sm.OLS(y, X)
    results1 = model1.fit()

    model2 = sm.OLS(y, X[:, 0:2])
    results2 = model2.fit()

    if results1.aic < results2.aic:
        prob_in_ci += float((results1.conf_int()[1][0] < beta[1]) & (results1.conf_int()[1][1] > beta[1])) / num_iter
    else:
        prob_in_ci += float((results2.conf_int()[1][0] < beta[1]) & (results2.conf_int()[1][1] > beta[1])) / num_iter

print('信頼区間に真の値が含まれる割合： ', prob_in_ci)

信頼区間に真の値が含まれる割合：  0.9202699999984466

まとめ

実験結果をまとめると以下です．

実験の結果

	$95$%信頼区間に 真値が含まれる割合
全変数利用	0.950
$x_{i1}$のみ利用	0.951
AICによるモデル選択	0.920

この結果から，全変数を利用したモデルも$x_{i1}$のみ利用したモデルも$95$%信頼区間の被覆確率は$95$%になるが，AICによるモデル選択という処理を挟むと，被覆確率が$95$%ではなくなることが分かります．
直感的には，AICにより当てはまりが良いモデルが選択されることで，信頼区間が本来よりも狭くなっていると理解できます．
この現象は，AICだけでなく，より一般に，同じデータを用いたモデル選択により発生します．
モデル選択が推定量・信頼区間に与える影響とその補正方法はPost Selection InferenceやSelective Inferenceと呼ばれる分野で研究されており，現在においてもホットなトピックです．