あなたが不動産の価格予測モデルを構築しているとしましょう。手元には「面積」「築年数」「駅からの距離」「階数」「周辺の学校数」「日照時間」「犯罪率」「公園までの距離」「バス路線数」「コンビニの数」など、数十もの候補変数があります。これらの全てをモデルに投入すべきでしょうか、それとも予測に本当に役立つ変数だけを選ぶべきでしょうか。

全変数を投入すると、予測に無関係な変数のノイズまでモデルが学習してしまい（過学習）、新しいデータに対する予測精度が低下します。一方、重要な変数を見落とすと、予測に必要な情報が失われます。変数選択（variable selection, feature selection）とは、予測に有用な変数のみを選び出し、簡潔で高精度なモデルを構築する技術です。

変数選択の理論と実践は、以下のような幅広い分野で重要です。

• 医学・疫学: 数百の候補バイオマーカーから疾患リスクの予測に有効な少数の変数を同定する
• 経済学: GDP成長率の予測に寄与するマクロ経済指標を特定する
• 工学: 製品品質に影響する工程パラメータを特定し、効率的な品質管理を行う
• 環境科学: 生態系の状態を予測する上で重要な環境変数を選択する

本記事では、変数選択の古典的手法であるステップワイズ法を3つのバリエーション（前進選択法、後退消去法、双方向ステップワイズ法）に分けて解説し、その選択基準としてAIC・BICの理論を導出します。さらにステップワイズ法の理論的・実用的な問題点を明らかにし、代替手法との比較を行います。

本記事の内容

• 変数選択問題の定式化 — なぜ全探索は難しいのか
• 前進選択法・後退消去法・双方向ステップワイズ法のアルゴリズム
• 選択基準: F検定・AIC・BICの理論と導出
• ステップワイズ法の6つの問題点
• 代替手法の概要と比較（正則化、ベストサブセット選択、交差検証）
• Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 重回帰分析の理論 — 回帰モデルの基本的な枠組み
• リッジ回帰とLASSO — 正則化による変数選択の代替アプローチ

変数選択問題の定式化

モデル空間の爆発

$p$ 個の候補説明変数 $x_1, x_2, \dots, x_p$ があるとき、各変数を「使う/使わない」の2択で考えると、可能なモデルの数は $2^p$ 通りです（空のモデルを除くと $2^p – 1$ 通り）。

$p = 20$ でもすでに100万通りを超え、全てのモデルを評価するのは現実的ではありません。$p = 50$ になると天文学的な数になり、完全探索は不可能です。

この組合せ爆発の問題を効率的に（ただし近似的に）解決するのがステップワイズ法の役割です。

最適モデルの基準

全探索が困難な場合でも、「どのモデルが良いか」を判断する基準は必要です。よく使われる基準には以下があります。

• 残差平方和（RSS）: データへの適合度を測る。小さいほど良いが、変数を増やすほど必ず減少するため、そのままでは過学習を防げない
• 自由度調整済み決定係数 $R^2_{\text{adj}}$: 変数の数に応じたペナルティを課した決定係数
• AIC（赤池情報量規準）: 予測精度を重視する情報理論的基準
• BIC（ベイズ情報量規準）: モデルの一致性を重視するベイズ的基準
• 交差検証誤差（CV）: データ駆動的な予測性能評価

AICとBICは理論的に明確な根拠を持ち、かつ計算効率が高いため、ステップワイズ法の基準として最もよく使われます。これらの詳細は後のセクションで解説します。

まずは、ステップワイズ法の具体的なアルゴリズムを見ていきましょう。

前進選択法（Forward Selection）

アルゴリズム

前進選択法は、最もシンプルなステップワイズ系の手法です。空のモデル（切片のみ）から始めて、1ステップごとに最も「有益な」変数を1つずつ追加していきます。

アルゴリズム:

1. 初期化: 空のモデル $M_0 = \{(\text{切片のみ})\}$ から開始。選択済み変数集合 $S = \emptyset$、候補変数集合 $C = \{x_1, \dots, x_p\}$
2. 変数評価: 候補集合 $C$ の各変数 $x_j$ について、現在のモデルに $x_j$ を追加したモデルの評価基準（AIC, BIC, F検定統計量など）を計算する
3. 変数追加の判定: 最も評価基準を改善する変数 $x_{j^*}$ を選ぶ。改善が基準を満たす場合（例: AICが減少する、F検定のp値が有意水準以下）、$x_{j^*}$ をモデルに追加し、$S \leftarrow S \cup \{x_{j^*}\}$, $C \leftarrow C \setminus \{x_{j^*}\}$ とする
4. 終了判定: 追加しても改善が見られなくなったら終了。そうでなければステップ2に戻る

F検定による判定

伝統的には、前進選択の各ステップでF検定を用いて変数追加の有意性を判定します。現在のモデルに変数 $x_j$ を追加したときの偏F統計量は:

$$ \begin{equation} F_j = \frac{(\text{RSS}_{\text{current}} – \text{RSS}_{\text{+}j}) / 1}{\text{RSS}_{\text{+}j} / (n – k – 2)} \end{equation} $$

ここで $\text{RSS}_{\text{current}}$ は現在のモデルの残差平方和、$\text{RSS}_{\text{+}j}$ は $x_j$ を追加したモデルの残差平方和、$k$ は現在のモデルの変数の数です。帰無仮説のもとで $F_j \sim F(1, n-k-2)$ に従います。

$F_j$ が閾値 $F_{\text{in}}$（例えば $F_{0.05}(1, n-k-2)$）を超えたら変数を追加し、全ての候補変数で閾値を超えなくなったら停止します。

前進選択法の特徴

前進選択法には次の利点があります。

• 計算効率: 最大 $p$ ステップで完了し、各ステップで最大 $p$ 個の候補を評価するため、計算量は $O(p^2)$（モデル適合の計算を除く）
• 実装の容易さ: アルゴリズムが単純で理解しやすい
• $p > n$ でも適用可能: 変数を1つずつ追加するため、$p$ がサンプル数を超えていても使える（後退消去法はこの場合不可）

一方、前進選択法の限界として、一度追加した変数は削除されないという点があります。ステップ3で変数 $x_j$ を追加した後、他の変数が追加された結果 $x_j$ が不要になっても、$x_j$ はモデルに残り続けます。この問題を解決するのが後退消去法と双方向ステップワイズ法です。

後退消去法（Backward Elimination）

アルゴリズム

後退消去法は、前進選択法とは反対に、全変数を含むフルモデルから始めて、1ステップごとに最も「不要な」変数を1つずつ削除していきます。

アルゴリズム:

1. 初期化: フルモデル $M_{\text{full}} = \{x_1, \dots, x_p\}$ から開始
2. 変数評価: モデル内の各変数 $x_j$ について、$x_j$ を削除したモデルの評価基準を計算する
3. 変数削除の判定: 削除しても評価基準の悪化が最も小さい変数 $x_{j^*}$ を選ぶ。悪化が許容範囲内であれば（例: AICが増加しない、F検定のp値が有意水準を超える）、$x_{j^*}$ をモデルから除外する
4. 終了判定: これ以上削除できなくなったら終了。そうでなければステップ2に戻る

F検定による判定の場合、偏F統計量が閾値 $F_{\text{out}}$ を下回る変数を削除します。

後退消去法の特徴

後退消去法の利点は:

• 変数間の交互作用を考慮: フルモデルから始めるため、変数間の共同効果が最初から考慮される
• 抑制変数の発見: ある変数が他の変数の存在下で初めて重要性を示す場合（抑制変数）を発見できる

限界は:

• $p > n$ では使えない: フルモデルの推定に $p < n$ が必要（最小二乗法の可解性条件）
• 計算コスト: フルモデルから始めるため、$p$ が大きいと初期のモデル推定が重くなる
• 削除した変数は復帰しない: 前進選択法と対称的な問題

前進選択法と後退消去法はそれぞれ一方向的な探索であるため、最適解を見逃す可能性があります。双方向ステップワイズ法はこの問題を緩和します。

双方向ステップワイズ法（Stepwise Selection）

アルゴリズム

双方向ステップワイズ法は、前進選択と後退消去を組み合わせた方法です。各ステップで変数の追加と削除の両方を検討します。

アルゴリズム:

1. 初期化: 空のモデルまたは適当な初期モデルから開始
2. 追加候補の評価: モデルに含まれていない各変数について、追加した場合の評価基準（AIC等）を計算する
3. 削除候補の評価: モデルに含まれている各変数について、削除した場合の評価基準を計算する
4. 最良の操作を実行: 追加・削除の全候補の中で、評価基準を最も改善する操作（追加または削除）を実行する
5. 終了判定: これ以上改善する操作がなければ終了。そうでなければステップ2に戻る

この手法では、前進選択法で一度追加した変数が、後のステップで他の変数が追加された結果不要になった場合、削除されることがあります。これにより、前進選択法の「一方向性」の問題を緩和できます。

追加と削除の閾値

F検定を基準とする場合、追加の閾値 $F_{\text{in}}$ と削除の閾値 $F_{\text{out}}$ を設定する必要があります。通常 $F_{\text{in}} \geq F_{\text{out}}$ とし、追加の基準を削除の基準より厳しくすることで、変数が追加と削除を繰り返すサイクルを防ぎます。

AIC/BICを基準とする場合は、単に「AIC/BICを最も改善する操作を実行し、改善しなくなったら停止」というシンプルなルールになります。

これらのステップワイズ法で用いられるAICとBICの理論的背景を、次のセクションで詳しく見ていきましょう。

情報量規準の理論

AIC（赤池情報量規準）の導出

AIC（Akaike Information Criterion）は、赤池弘次が1973年に提案した情報量規準です。AICはモデルの予測性能を評価する基準であり、カルバック・ライブラー情報量（KLダイバージェンス）の推定に基づいています。

真の分布を $g(\bm{y})$、候補モデルの分布を $f(\bm{y}|\hat{\bm{\theta}})$ としたとき、モデルの良さはKLダイバージェンスで測られます。

$$ \begin{equation} D_{\text{KL}}(g \| f) = \int g(\bm{y})\ln\frac{g(\bm{y})}{f(\bm{y}|\hat{\bm{\theta}})}d\bm{y} \end{equation} $$

$D_{\text{KL}}$ の中で $g$ のみに依存する項（$\int g\ln g$）はモデル比較では定数なので、最小化すべきは $-\mathbb{E}_g[\ln f(\bm{y}|\hat{\bm{\theta}})]$ です。

赤池はこの期待対数尤度を推定するために、対数尤度 $\ell(\hat{\bm{\theta}}) = \ln f(\bm{y}|\hat{\bm{\theta}})$ からバイアス補正項を引くべきことを示しました。最大対数尤度は楽観的なバイアスを持ち、そのバイアスは漸近的にパラメータ数 $k$ に等しくなります。

$$ \begin{equation} \mathbb{E}[\ell(\hat{\bm{\theta}})] – \mathbb{E}_g[\mathbb{E}_{\bm{y}_{\text{new}}}[\ln f(\bm{y}_{\text{new}}|\hat{\bm{\theta}})]] \approx k \end{equation} $$

したがって、予測対数尤度の不偏推定量として次のAICが得られます。

$$ \begin{equation} \text{AIC} = -2\ell(\hat{\bm{\theta}}) + 2k \end{equation} $$

$-2$ のファクターは歴史的慣習であり、AICが小さいほど良いモデルです。

正規線形回帰モデルの場合、対数尤度は $\ell(\hat{\bm{\theta}}) = -\frac{n}{2}\ln(2\pi) – \frac{n}{2}\ln(\hat{\sigma}^2) – \frac{n}{2}$ です。定数項を除くと:

$$ \begin{equation} \text{AIC} = n\ln\left(\frac{\text{RSS}}{n}\right) + 2(p’ + 1) \end{equation} $$

ここで $p’$ は説明変数の数、$+1$ は分散パラメータ $\sigma^2$ の分を含めたものです。

BIC（ベイズ情報量規準）の導出

BIC（Bayesian Information Criterion）は、シュワルツ（Schwarz）が1978年にベイズ的観点から導出した基準です。

ベイズ的モデル選択では、事後確率が最大のモデルを選びます。モデル $M_k$ の事後確率は:

$$ P(M_k|\bm{y}) \propto P(\bm{y}|M_k)P(M_k) $$

$P(\bm{y}|M_k)$ は周辺尤度（エビデンス）であり、パラメータを積分消去した量です。

$$ P(\bm{y}|M_k) = \int f(\bm{y}|\bm{\theta}_k, M_k)\pi(\bm{\theta}_k|M_k)d\bm{\theta}_k $$

ラプラス近似（被積分関数を最大値の周りで2次近似する方法）を用いると、$\ln P(\bm{y}|M_k)$ は次のように近似されます。

$$ \ln P(\bm{y}|M_k) \approx \ell(\hat{\bm{\theta}}_k) + \frac{k}{2}\ln(2\pi) – \frac{1}{2}\ln|\bm{H}(\hat{\bm{\theta}}_k)| + \ln\pi(\hat{\bm{\theta}}_k) $$

ここで $\bm{H}$ はヘッセ行列です。$n \to \infty$ のとき $|\bm{H}| = O(n^k)$ なので、主要項のみを残すと:

$$ \begin{equation} \ln P(\bm{y}|M_k) \approx \ell(\hat{\bm{\theta}}_k) – \frac{k}{2}\ln n + O(1) \end{equation} $$

$-2$ を掛けてBICを定義します。

$$ \begin{equation} \text{BIC} = -2\ell(\hat{\bm{\theta}}) + k\ln n \end{equation} $$

AICとの違いはペナルティ項です。AICのペナルティは $2k$ ですが、BICのペナルティは $k\ln n$ です。$\ln n > 2$（つまり $n \geq 8$）のとき、BICはAICより重いペナルティを課し、より簡素なモデルを好みます。

AICとBICの理論的性質の比較

AICとBICは数学的に異なる目標を追求しています。

実用的な指針として:

• 予測が主目的の場合 → AICが適している
• 真の生成モデルの同定が主目的の場合 → BICが適している
• 判断に迷う場合 → 交差検証を併用して両基準の結果を比較する

情報量規準の理論を踏まえた上で、ステップワイズ法の問題点を見ていきましょう。

ステップワイズ法の問題点

ステップワイズ法は計算効率が高く直感的に理解しやすいため、長年にわたって広く使われてきました。しかし、近年の統計学では、ステップワイズ法には深刻な問題があることが広く認識されています。ここでは6つの主要な問題点を解説します。

問題1: 多重検定の問題

ステップワイズ法は各ステップで統計的検定を繰り返します。$p$ 個の候補変数に対して最大 $p$ ステップ、各ステップで最大 $p$ 回の検定を行うため、全体での検定回数は最大 $O(p^2)$ にもなります。

各ステップでの有意水準が $\alpha = 0.05$ であっても、$p = 20$ の場合に全体として少なくとも1つの偽陽性が生じる確率は $1 – (1 – 0.05)^{20} \approx 0.64$ にもなります。つまり、実際にはほぼ確実にいくつかの無関係な変数が誤って選択されます。

問題2: p値と信頼区間の歪み

ステップワイズ法で選択されたモデルのp値や信頼区間は、変数選択の過程を無視しています。選択後の推論（post-selection inference）は、選択バイアスにより通常のp値や信頼区間が無効になることが知られています。

具体的には、ステップワイズ法で選択された変数のp値は真のp値よりも過度に小さくなります（過度に有意に見える）。信頼区間も過度に狭くなり、カバレッジ確率が名目水準を大きく下回ります。

問題3: 最適解の非保証

ステップワイズ法は貪欲法（greedy algorithm）であり、各ステップで局所最適な選択を行います。しかし、各ステップでの最良の選択が全体としての最良の選択であるとは限りません。

例えば、$x_1$ と $x_2$ はそれぞれ単独では目的変数に対する説明力が弱いが、両方を同時にモデルに含めると高い説明力を示す場合（交互作用や抑制効果）、前進選択法はこの組み合わせを発見できません。

問題4: 不安定性

ステップワイズ法の結果は、データの小さな変動に対して非常に不安定であることが知られています。サンプルのわずかな変更（1つのデータ点の追加や除外）で、全く異なる変数の組み合わせが選択されることがあります。

この不安定性は、変数間に多重共線性がある場合に特に顕著です。相関の高い変数群の中でどの変数が選ばれるかは、データのわずかな変動に依存し、実質的に偶然の選択になってしまいます。

問題5: 自由度の過小評価

ステップワイズ法で選択されたモデルの有効自由度は、モデルに含まれる変数の数よりも大きくなります。変数選択の過程自体が追加の自由度を消費しているためです。

$p$ 個の候補変数から $k$ 個を選んだ場合、実効的な自由度は $k$ ではなく $k + c$（$c > 0$ は選択過程による追加自由度）です。この過小評価により、$R^2_{\text{adj}}$ やF検定の結果が過度に楽観的になります。

問題6: 変数の重要度の順位が歪む

ステップワイズ法は「最初に選ばれた変数が最も重要」という暗黙の解釈を与えがちですが、選択順序は変数の真の重要度を必ずしも反映しません。変数間の相関構造により、真に重要な変数が選ばれず、代わりに相関の高い別の変数が選ばれることがあります。

これらの問題は理論的に十分に文書化されており、多くの統計学者がステップワイズ法の使用に対して強い警告を発しています。では、どのような代替手法があるのでしょうか。

代替手法との比較

LASSO・Elastic Net（正則化手法）

LASSOとElastic Netは、目的関数にペナルティ項を追加することで連続的に変数選択を行います。

$$ \hat{\bm{\beta}} = \arg\min_{\bm{\beta}}\left\{\|\bm{y} – \bm{X}\bm{\beta}\|^2 + \lambda\|\bm{\beta}\|_1\right\} $$

正則化手法の利点:

• ペナルティの強さ $\lambda$ は交差検証で客観的に選べる
• 全ての変数を同時に考慮するため、変数間の関係が自然に反映される
• $p > n$ の高次元問題にも適用可能
• 理論的な保証（oracle不等式など）がある

ステップワイズ法との比較では、正則化手法の方がほぼ全ての面で優れています。

ベストサブセット選択

ベストサブセット選択は、全ての $2^p$ 通りのモデルを（何らかの方法で）評価し、AIC/BICが最小のモデルを選ぶ方法です。

計算量が問題ですが、近年の研究で混合整数最適化（MIO: Mixed Integer Optimization）を用いることで、$p$ が数百の問題でも実用的な時間内にベストサブセット選択が可能になりました（Bertsimas, King, Mazumder 2016）。

ベストサブセット選択はステップワイズ法の「最適解非保証」問題を解決しますが、多重検定やp値の歪みの問題は残ります。

交差検証に基づく選択

ステップワイズ法やAIC/BICの代わりに、$K$-fold 交差検証を直接的な選択基準として使う方法もあります。各候補モデルの予測誤差を交差検証で推定し、予測誤差が最小のモデルを選びます。

交差検証の利点は、AIC/BICのような漸近近似に頼らず、有限サンプルでの予測性能を直接評価できる点です。一方、計算コストが高い（各モデルに対して $K$ 回の学習が必要）という欠点があります。

手法の比較表

それでは、Pythonでこれらの手法を実装し、動作を比較してみましょう。

Pythonによる実装と可視化

データ生成と前進選択法の実装

import numpy as np
import matplotlib.pyplot as plt
from itertools import combinations

np.random.seed(42)

# --- データ生成 ---
n = 100
p_total = 10

# 真のモデル: y = 3*x1 - 2*x2 + 1.5*x3 + noise
# x4-x10 は無関係なノイズ変数
X = np.random.randn(n, p_total)
beta_true = np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0, 0, 0])
y = X @ beta_true + np.random.randn(n)

def compute_aic(X_sub, y, n):
    """サブモデルのAICを計算"""
    X_design = np.column_stack([np.ones(n), X_sub])
    beta_hat = np.linalg.lstsq(X_design, y, rcond=None)[0]
    RSS = np.sum((y - X_design @ beta_hat)**2)
    k = X_design.shape[1]
    aic = n * np.log(RSS / n) + 2 * k
    return aic, RSS, beta_hat

def compute_bic(X_sub, y, n):
    """サブモデルのBICを計算"""
    X_design = np.column_stack([np.ones(n), X_sub])
    beta_hat = np.linalg.lstsq(X_design, y, rcond=None)[0]
    RSS = np.sum((y - X_design @ beta_hat)**2)
    k = X_design.shape[1]
    bic = n * np.log(RSS / n) + k * np.log(n)
    return bic

# --- 前進選択法（AIC基準） ---
selected = []
remaining = list(range(p_total))
forward_history = []

# 空モデルのAIC
X_null = np.ones((n, 1))
beta_null = np.linalg.lstsq(X_null, y, rcond=None)[0]
RSS_null = np.sum((y - X_null @ beta_null)**2)
aic_null = n * np.log(RSS_null / n) + 2 * 1
forward_history.append(({'selected': [], 'aic': aic_null}))

for step in range(p_total):
    best_aic = np.inf
    best_var = None

    for var in remaining:
        if selected:
            X_sub = X[:, selected + [var]]
        else:
            X_sub = X[:, [var]]
        aic, _, _ = compute_aic(X_sub, y, n)
        if aic < best_aic:
            best_aic = aic
            best_var = var

    current_aic = forward_history[-1]['aic']
    if best_var is not None and best_aic < current_aic:
        selected.append(best_var)
        remaining.remove(best_var)
        forward_history.append({'selected': list(selected), 'aic': best_aic})
    else:
        break

print("=== 前進選択法（AIC基準） ===")
for i, h in enumerate(forward_history):
    vars_str = [f"x{v+1}" for v in h['selected']]
    print(f"Step {i}: 変数={vars_str}, AIC={h['aic']:.2f}")
print(f"最終選択: {[f'x{v+1}' for v in selected]}")
print(f"真の非ゼロ変数: x1, x2, x3")

前進選択法の結果から、重要なことが確認できます。まず、AIC基準での前進選択は、最初の3ステップで真に重要な変数 $x_1, x_2, x_3$ を正しく選択していることが期待されます。各ステップでAICが改善し、それ以降は追加してもAICが改善しなくなるため停止します。この例では、ノイズ変数が偶然に有意に見えない限り、正しいモデルが選択されます。

ベストサブセット選択との比較

import numpy as np
import matplotlib.pyplot as plt
from itertools import combinations

np.random.seed(42)
n = 100
p_total = 10
X = np.random.randn(n, p_total)
beta_true = np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0, 0, 0])
y = X @ beta_true + np.random.randn(n)

# --- ベストサブセット選択（全探索） ---
all_results = []
for k in range(1, p_total + 1):
    for subset in combinations(range(p_total), k):
        X_sub = np.column_stack([np.ones(n), X[:, list(subset)]])
        beta_hat = np.linalg.lstsq(X_sub, y, rcond=None)[0]
        RSS = np.sum((y - X_sub @ beta_hat)**2)
        n_params = len(subset) + 1
        aic = n * np.log(RSS / n) + 2 * n_params
        bic = n * np.log(RSS / n) + n_params * np.log(n)
        all_results.append({
            'subset': subset, 'k': len(subset),
            'RSS': RSS, 'AIC': aic, 'BIC': bic
        })

# 各kでの最良モデルを抽出
best_by_k_aic = {}
best_by_k_bic = {}
for r in all_results:
    k = r['k']
    if k not in best_by_k_aic or r['AIC'] < best_by_k_aic[k]['AIC']:
        best_by_k_aic[k] = r
    if k not in best_by_k_bic or r['BIC'] < best_by_k_bic[k]['BIC']:
        best_by_k_bic[k] = r

# 全体最良
best_aic = min(all_results, key=lambda x: x['AIC'])
best_bic = min(all_results, key=lambda x: x['BIC'])

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# (a) AIC/BICの変化
ax = axes[0]
ks = sorted(best_by_k_aic.keys())
aic_vals = [best_by_k_aic[k]['AIC'] for k in ks]
bic_vals = [best_by_k_bic[k]['BIC'] for k in ks]
ax.plot(ks, aic_vals, 'bo-', linewidth=2, markersize=7, label='Best AIC')
ax.plot(ks, bic_vals, 'rs-', linewidth=2, markersize=7, label='Best BIC')
ax.axvline(3, color='gray', linestyle='--', alpha=0.7, label='True # of vars')
ax.set_xlabel('Number of variables', fontsize=11)
ax.set_ylabel('Information criterion', fontsize=11)
ax.set_title('(a) Best AIC/BIC by model size', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) 全サブセットのAIC散布図
ax = axes[1]
for r in all_results:
    ax.scatter(r['k'], r['AIC'], color='lightblue', alpha=0.2, s=10)
for k in ks:
    ax.scatter(k, best_by_k_aic[k]['AIC'], color='blue', s=50, zorder=5)
ax.scatter(best_aic['k'], best_aic['AIC'], color='red', s=100, marker='*',
           zorder=10, label=f"Best: vars={best_aic['subset']}")
ax.set_xlabel('Number of variables', fontsize=11)
ax.set_ylabel('AIC', fontsize=11)
ax.set_title('(b) All subsets AIC', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\n=== ベストサブセット選択の結果 ===")
print(f"AIC最良モデル: 変数{tuple(v+1 for v in best_aic['subset'])}, "
      f"AIC={best_aic['AIC']:.2f}")
print(f"BIC最良モデル: 変数{tuple(v+1 for v in best_bic['subset'])}, "
      f"BIC={best_bic['BIC']:.2f}")

ベストサブセット選択の結果から、以下の知見が得られます。

1. 図(a): AICとBICの両方が変数数3で最小値を取り、真のモデルの変数数と一致しています。変数数が3を超えると、不要な変数のペナルティにより情報量規準が悪化します。BICはAICよりも急激に上昇しており、不要な変数に対するペナルティが重いことが視覚的にわかります

2. 図(b): 同じ変数数でも、どの変数の組み合わせを選ぶかによってAICに大きな差があります。変数数3のモデルの中でも、真の変数$(x_1, x_2, x_3)$ を含むモデルだけがAICの最小値を達成し、他の3変数の組み合わせはAICがはるかに大きくなっています

ステップワイズ法の不安定性の実証

ステップワイズ法の不安定性を、ブートストラップ法で確認します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 100
p_total = 10
X = np.random.randn(n, p_total)
beta_true = np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0, 0, 0])
y = X @ beta_true + np.random.randn(n)

def forward_selection_aic(X, y):
    """前進選択法（AIC基準）"""
    n = X.shape[0]
    p = X.shape[1]
    selected = []
    remaining = list(range(p))

    X_null = np.ones((n, 1))
    beta_null = np.linalg.lstsq(X_null, y, rcond=None)[0]
    current_aic = n * np.log(np.sum((y - X_null @ beta_null)**2) / n) + 2

    for _ in range(p):
        best_aic, best_var = np.inf, None
        for var in remaining:
            X_sub = np.column_stack([np.ones(n), X[:, selected + [var]]])
            beta_hat = np.linalg.lstsq(X_sub, y, rcond=None)[0]
            RSS = np.sum((y - X_sub @ beta_hat)**2)
            aic = n * np.log(RSS / n) + 2 * X_sub.shape[1]
            if aic < best_aic:
                best_aic, best_var = aic, var
        if best_var is not None and best_aic < current_aic:
            selected.append(best_var)
            remaining.remove(best_var)
            current_aic = best_aic
        else:
            break
    return selected

# ブートストラップ実験
n_bootstrap = 500
selection_counts = np.zeros(p_total)

for b in range(n_bootstrap):
    idx = np.random.choice(n, size=n, replace=True)
    X_boot, y_boot = X[idx], y[idx]
    selected = forward_selection_aic(X_boot, y_boot)
    for var in selected:
        selection_counts[var] += 1

selection_freq = selection_counts / n_bootstrap

fig, ax = plt.subplots(figsize=(10, 5))
colors = ['#e41a1c' if beta_true[j] != 0 else '#377eb8' for j in range(p_total)]
bars = ax.bar(range(1, p_total + 1), selection_freq, color=colors, alpha=0.7)
ax.set_xlabel('Variable', fontsize=11)
ax.set_ylabel('Selection frequency', fontsize=11)
ax.set_title('Bootstrap selection frequency (Forward Selection, 500 replicates)',
             fontsize=12)
ax.set_xticks(range(1, p_total + 1))
ax.set_xticklabels([f'$x_{{{j+1}}}$' for j in range(p_total)])
ax.axhline(1.0, color='gray', linestyle='--', alpha=0.5)

# 凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='#e41a1c', alpha=0.7, label='True variable'),
                   Patch(facecolor='#377eb8', alpha=0.7, label='Noise variable')]
ax.legend(handles=legend_elements, fontsize=10)
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print("=== 選択頻度 ===")
for j in range(p_total):
    status = "TRUE" if beta_true[j] != 0 else "noise"
    print(f"x{j+1} ({status}): {selection_freq[j]:.1%}")

ブートストラップ実験の結果は、ステップワイズ法の不安定性を定量的に示しています。

1. 真の変数（赤）: $x_1, x_2, x_3$ はほぼ全てのブートストラップ標本で選択されていますが、100%ではありません。特に効果の小さい $x_3$（$\beta_3 = 1.5$）は、$x_1$（$\beta_1 = 3.0$）よりも選択頻度が低くなる傾向があります

2. ノイズ変数（青）: 理想的にはゼロに近い選択頻度であるべきですが、実際にはいくつかのノイズ変数が5〜15%の頻度で選択されています。500回のブートストラップのうち数十回は誤って選択されているということです。これが多重検定問題の影響です

3. 解釈: サンプルのわずかな変動でモデルが変わりうることがわかります。特にノイズ変数の選択はデータ依存的で不安定であり、ステップワイズ法の結果を過度に信頼すべきでないことを示しています

まとめ

本記事では、変数選択法の理論と実践を包括的に解説しました。

• ステップワイズ法: 前進選択法・後退消去法・双方向法の3バリエーション。計算効率は高いが、局所最適にとどまる貪欲法である
• AIC: KLダイバージェンスの推定に基づく予測志向の基準。ペナルティは $2k$。効率性を持つ
• BIC: ベイズ周辺尤度のラプラス近似に基づく基準。ペナルティは $k\ln n$。一致性を持つ
• 問題点: 多重検定問題、p値と信頼区間の歪み、最適解の非保証、不安定性、自由度の過小評価、重要度順位の歪み
• 推奨代替手法: LASSO/Elastic Net（連続的な変数選択）、ベストサブセット選択（完全探索）、交差検証（予測性能の直接評価）

現代の統計学・機械学習では、ステップワイズ法よりも正則化手法（特にLASSOやElastic Net）が推奨されることが多くなっています。ステップワイズ法を使う場合でも、その結果の不安定性を認識し、ブートストラップやクロスバリデーションによる確認を併用すべきです。

次のステップとして、以下の記事も参考にしてください。

• リッジ回帰とLASSOの理論 — 正則化による変数選択の理論的基盤
• Elastic Net — L1とL2を組み合わせた正則化手法
• 多重共線性とVIF — 変数間の相関の診断方法