はじめに
回帰分析って知ってますか？
こんな質問されたらゾッとしますよね。知らないというわけでもないですが、きっとかなり深いレベルで知っていないとボコボコにされるでしょう。
R言語の lm() や Python の sklearn.LinearRegression() を呼び出して、推定した回帰係数を眺めたことがある、紙とペンでOLSによる係数を求めたことがある、OLS推定量がBLUEであることを確かめた、因果効果を回帰分析で推定した。。。どこまでやれば回帰分析を「知っている」と言って良いのでしょうか。
そんなことは何でも良いのです。とにかく回帰分析はシンプルながら、最も奥深いテーマの一つであり、多層的に理解する必要があります。
この記事では、回帰分析の理解をLevel 0（道具として使うだけ）からLevel 5（高次元データにおけるDMLなど）まで 、6段階に分けて解説します。各レベルで「何が見えるようになるか」「どこに崖があるか」を丁寧に追っていきます。
思いの外分量のある記事になってしまったので、対象ごとにお勧めの章を整理しておきます。
回帰分析への理解は、ダニングクルーガー効果の絶望の谷を味わい続けるのにもってこいです。さて、本記事を通じて一度「絶望の谷」に落ちて、そこから一緒に「啓発の坂」を登っていきましょう。
https://hrmos.co/trend/talent-management/5397/ より引用。
Level 0: 道具としての回帰
最初の段階は、回帰をブラックボックスのツール として使う段階です。
この程度の実装はAIで書けるレベルを通り越してもはやボイラープレート で、書けること自体に価値はありません。
回帰分析をこんな感じで捉えている人は多いでしょう：
年収 = β 0 + β 1 ⋅ 学歴 + β 2 ⋅ 経験年数 + ε
\text{年収} = \beta_0 + \beta_1 \cdot \text{学歴} + \beta_2 \cdot \text{経験年数} + \varepsilon
年収 = β 0 + β 1 ⋅ 学歴 + β 2 ⋅ 経験年数 + ε 回帰係数の符号と大きさを見て「学歴は効いてる」「経験年数も正に効いてる」と語ります。R²が0.7なら「全体的にまあまあ説明できてる」と判断する感じです。
この段階の限界は明らかで、以下の問いに答えられません。
係数が現実世界において何を意味するのか
いつどんな時にその係数を信じていいのか
「効いてる」とはどういう意味か
ここから先のレベルで、これらの問いに段階的に答えられるようになっていきます。
Level 1: 古典的な統計的解釈
このレベルは簡単な統計学の教科書を一通り学び終えた水準です。回帰分析の理論は習得しているものの、深くは踏み込んでいない、統計検定でいうところの2級水準のイメージです。具体的には以下の事柄を理解している状態です。
最小二乗法の幾何学
線形回帰モデルを行列で書きます。
y = X β + ε
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}
y = X β + ε ここで y ∈ R n \mathbf{y} \in \mathbb{R}^n y ∈ R n X ∈ R n × p \mathbf{X} \in \mathbb{R}^{n \times p} X ∈ R n × p β ∈ R p \boldsymbol{\beta} \in \mathbb{R}^p β ∈ R p p p p n n n y y y
最小二乗推定量（OLS）は残差二乗和を最小化することで得られます（今後断りなしに推定したパラメータにはハット（^）をつけます）。
β ^ = arg min β ∥ y − X β ∥ 2 = ( X ⊤ X ) − 1 X ⊤ y
\hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y}
β ^ = arg β min ∥ y − X β ∥ 2 = ( X ⊤ X ) − 1 X ⊤ y OLSの計算詳細
1. 目的関数の定義と展開
最小化したい目的関数（残差平方和）を S ( β ) S(\boldsymbol{\beta}) S ( β )
S ( β ) = ∥ y − X β ∥ 2 = ( y − X β ) ⊤ ( y − X β )
\begin{aligned}
S(\boldsymbol{\beta}) &= \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 \\
&= (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})
\end{aligned}
S ( β ) = ∥ y − X β ∥ 2 = ( y − X β ) ⊤ ( y − X β ) これを展開します。転置の性質 ( A B ) ⊤ = B ⊤ A ⊤ (\mathbf{A}\mathbf{B})^\top = \mathbf{B}^\top\mathbf{A}^\top ( AB ) ⊤ = B ⊤ A ⊤
S ( β ) = ( y ⊤ − ( X β ) ⊤ ) ( y − X β ) = ( y ⊤ − β ⊤ X ⊤ ) ( y − X β ) = y ⊤ y − y ⊤ X β − β ⊤ X ⊤ y + β ⊤ X ⊤ X β
\begin{aligned}
S(\boldsymbol{\beta}) &= (\mathbf{y}^\top - (\mathbf{X}\boldsymbol{\beta})^\top) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\
&= (\mathbf{y}^\top - \boldsymbol{\beta}^\top \mathbf{X}^\top) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \\
&= \mathbf{y}^\top \mathbf{y} - \mathbf{y}^\top \mathbf{X}\boldsymbol{\beta} - \boldsymbol{\beta}^\top \mathbf{X}^\top \mathbf{y} + \boldsymbol{\beta}^\top \mathbf{X}^\top \mathbf{X}\boldsymbol{\beta}
\end{aligned}
S ( β ) = ( y ⊤ − ( X β ) ⊤ ) ( y − X β ) = ( y ⊤ − β ⊤ X ⊤ ) ( y − X β ) = y ⊤ y − y ⊤ X β − β ⊤ X ⊤ y + β ⊤ X ⊤ X β ここで、第2項と第3項はどちらもスカラー（1 × 1 1 \times 1 1 × 1 y ⊤ X β = ( y ⊤ X β ) ⊤ = β ⊤ X ⊤ y \mathbf{y}^\top \mathbf{X}\boldsymbol{\beta} = (\mathbf{y}^\top \mathbf{X}\boldsymbol{\beta})^\top = \boldsymbol{\beta}^\top \mathbf{X}^\top \mathbf{y} y ⊤ X β = ( y ⊤ X β ) ⊤ = β ⊤ X ⊤ y
S ( β ) = y ⊤ y − 2 β ⊤ X ⊤ y + β ⊤ X ⊤ X β
S(\boldsymbol{\beta}) = \mathbf{y}^\top \mathbf{y} - 2\boldsymbol{\beta}^\top \mathbf{X}^\top \mathbf{y} + \boldsymbol{\beta}^\top \mathbf{X}^\top \mathbf{X}\boldsymbol{\beta}
S ( β ) = y ⊤ y − 2 β ⊤ X ⊤ y + β ⊤ X ⊤ X β
2. 行列による微分
この関数 S ( β ) S(\boldsymbol{\beta}) S ( β ) β \boldsymbol{\beta} β 0 0 0
行列微分の公式
∂ ( a ⊤ x ) ∂ x = a \frac{\partial (\mathbf{a}^\top \mathbf{x})}{\partial \mathbf{x}} = \mathbf{a} ∂ x ∂ ( a ⊤ x ) = a
∂ ( x ⊤ A x ) ∂ x = ( A + A ⊤ ) x \frac{\partial (\mathbf{x}^\top \mathbf{A} \mathbf{x})}{\partial \mathbf{x}} = (\mathbf{A} + \mathbf{A}^\top)\mathbf{x} ∂ x ∂ ( x ⊤ Ax ) = ( A + A ⊤ ) x A \mathbf{A} A 2 A x 2\mathbf{Ax} 2 Ax
これらを適用すると：
∂ S ( β ) ∂ β = − 2 X ⊤ y + 2 X ⊤ X β
\frac{\partial S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = -2\mathbf{X}^\top \mathbf{y} + 2\mathbf{X}^\top \mathbf{X} \boldsymbol{\beta}
∂ β ∂ S ( β ) = − 2 X ⊤ y + 2 X ⊤ X β ※ X ⊤ X \mathbf{X}^\top \mathbf{X} X ⊤ X 2 A x 2\mathbf{Ax} 2 Ax
3. 正規方程式の解
勾配を 0 0 0 β \boldsymbol{\beta} β β ^ \hat{\boldsymbol{\beta}} β ^
− 2 X ⊤ y + 2 X ⊤ X β ^ = 0 X ⊤ X β ^ = X ⊤ y （正規方程式）
\begin{aligned}
-2\mathbf{X}^\top \mathbf{y} + 2\mathbf{X}^\top \mathbf{X} \hat{\boldsymbol{\beta}} &= \mathbf{0} \\
\mathbf{X}^\top \mathbf{X} \hat{\boldsymbol{\beta}} &= \mathbf{X}^\top \mathbf{y} \quad \text{（正規方程式）}
\end{aligned}
− 2 X ⊤ y + 2 X ⊤ X β ^ X ⊤ X β ^ = 0 = X ⊤ y （正規方程式） 最後に、X ⊤ X \mathbf{X}^\top \mathbf{X} X ⊤ X
β ^ = ( X ⊤ X ) − 1 X ⊤ y
\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y}
β ^ = ( X ⊤ X ) − 1 X ⊤ y このOLS計算では予測誤差の二乗和を最小にする係数を求めていますが、幾何学的に理解することも重要です。
回帰分析の幾何学的理解。
このアニメーションを説明していきます。
空間設定
このアニメーションの空間は n n n p p p
観測ベクトル y y y n n n n n n R n \mathbb{R}^n R n
データ行列 X X X n × p n \times p n × p x 1 , x 2 , … , x p x_1, x_2, \dots, x_p x 1 , x 2 , … , x p R n \mathbb{R}^n R n
まずは黒矢印に着目してください。これは基底ベクトルです（x 1 , x 2 x_1, x_2 x 1 , x 2 p p p p = 2 p=2 p = 2 X \mathbf{X} X n n n M ( X ) M(X) M ( X )
念の為の補足ですが、これはあるデータ y i y_i y i n n n
データ生成プロセス
次にデータ生成プロセスを考えます。
真の平均 η \eta η
モデル内成分 η 0 \eta_0 η 0 x 1 , x 2 x_1, x_2 x 1 , x 2
バイアス η 1 \eta_1 η 1 y y y x 3 , x 4 x_3, x_4 x 3 , x 4
η = η 0 + η 1 \eta = \eta_0 + \eta_1 η = η 0 + η 1 バイアス と呼ばれるものです。
観測値 y y y η \eta η ϵ \epsilon ϵ
最小二乗法による推定（直交射影）
私たちがやりたいことは、手元にある赤色の矢印 y y y
推定値 y ^ \hat{y} y ^ y y y M ( X ) M(X) M ( X ) y y y バリアンス と呼びます。
なぜ垂直なのか？:y y y
幾何学的視点から見たバイアス・バリアンストレードオフ
説明変数が少なすぎると、モデルの表現力不足に陥ります。真の平均値と、モデルで表現できる平面との間の距離が開いてバイアスが大きくなります。一方、変数を増やしすぎてモデルが複雑になると、ノイズに反応しやすくなってしまいます。その分、平面自体が真の値に近付いてくれるはずなのでバイアスは小さくなります。
回帰分析の道具
この段階の理解をしている人は以下の事柄も理解していることでしょう。
検定と信頼区間
https://zenn.dev/zenkigen_tech/articles/8ffcd1b39e1e5f
最小二乗法そのものは、誤差項の分布に特定の仮定を置かなくても計算可能です。しかし、推定された係数が「偶然得られたものではないか」を統計的に検定するためには、誤差項に関する強い仮定が必要になります。
ここで、誤差項ベクトル ε \boldsymbol{\varepsilon} ε 0 \mathbf{0} 0 σ 2 \sigma^2 σ 2
ε ∼ N ( 0 , σ 2 I )
\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I})
ε ∼ N ( 0 , σ 2 I ) この「正規性の仮定」を置くことで、推定量 β ^ = ( X ⊤ X ) − 1 X ⊤ y \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} β ^ = ( X ⊤ X ) − 1 X ⊤ y
β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 )
\hat{\boldsymbol{\beta}} \sim \mathcal{N}(\boldsymbol{\beta}, \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1})
β ^ ∼ N ( β , σ 2 ( X ⊤ X ) − 1 ) この性質を起点に、回帰分析における代表的な2つの検定が構成されます。
1. 個別係数のt検定
「ある特定の説明変数が、目的変数に対して本当に影響を与えているか」を確認するための検定です。
通常、帰無仮説 H 0 : β j = 0 H_0: \beta_j = 0 H 0 : β j = 0 x j x_j x j σ 2 \sigma^2 σ 2 σ ^ 2 \hat{\sigma}^2 σ ^ 2 S E ( β ^ j ) \mathrm{SE}(\hat{\beta}_j) SE ( β ^ j )
t = β ^ j − 0 S E ( β ^ j ) ∼ t ( n − k − 1 )
t = \frac{\hat{\beta}_j - 0}{\mathrm{SE}(\hat{\beta}_j)} \sim t(n - k - 1)
t = SE ( β ^ j ) β ^ j − 0 ∼ t ( n − k − 1 ) この t t t n − k − 1 n - k - 1 n − k − 1 n n n k k k 0.05 0.05 0.05
2. モデル全体の有意性を評価する F検定
「構築した回帰モデル全体が、単なる定数モデル（平均値予測）よりも意味のある予測をしているか」を確認するための検定です。
帰無仮説として、定数項以外のすべての回帰係数がゼロである（H 0 : β 1 = β 2 = ⋯ = β k = 0 H_0: \beta_1 = \beta_2 = \dots = \beta_k = 0 H 0 : β 1 = β 2 = ⋯ = β k = 0 R 2 R^2 R 2 F F F
F = R 2 / k ( 1 − R 2 ) / ( n − k − 1 ) ∼ F ( k , n − k − 1 )
F = \frac{R^2 / k}{(1 - R^2) / (n - k - 1)} \sim F(k, n - k - 1)
F = ( 1 − R 2 ) / ( n − k − 1 ) R 2 / k ∼ F ( k , n − k − 1 ) この統計量は、自由度 ( k , n − k − 1 ) (k, n - k - 1) ( k , n − k − 1 ) β 1 = β 2 \beta_1 = \beta_2 β 1 = β 2
自由度調整済み決定係数
決定係数 R 2 R^2 R 2
R 2 = 1 − ∑ i ( y i − y ^ i ) 2 ∑ i ( y i − y ˉ ) 2 = 1 − RSS TSS
R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2} = 1 - \frac{\text{RSS}}{\text{TSS}}
R 2 = 1 − ∑ i ( y i − y ˉ ) 2 ∑ i ( y i − y ^ i ) 2 = 1 − TSS RSS ただし R 2 R^2 R 2 R 2 R^2 R 2 自由度調整済み決定係数 です：
R adj 2 = 1 − RSS / ( n − p ) TSS / ( n − 1 )
R^2_{\text{adj}} = 1 - \frac{\text{RSS}/(n-p)}{\text{TSS}/(n-1)}
R adj 2 = 1 − TSS / ( n − 1 ) RSS / ( n − p ) n − k − 1 n - k - 1 n − k − 1 k k k R adj 2 {R}_{\text{adj}}^2 R adj 2
多重共線性
説明変数同士が強く相関していると、X ⊤ X \mathbf{X}^\top\mathbf{X} X ⊤ X ( X ⊤ X ) − 1 (\mathbf{X}^\top\mathbf{X})^{-1} ( X ⊤ X ) − 1
有名な診断指標は VIF （Variance Inflation Factor） です。
VIF j = 1 1 − R j 2
\text{VIF}_j = \frac{1}{1 - R_j^2}
VIF j = 1 − R j 2 1 ここで R j 2 R_j^2 R j 2 j j j
残差診断
最小二乗法（OLS）による推定が適切であるためには、誤差項に関して「線形性・等分散性・正規性・独立性」といった前提条件が満たされている必要があります。これらを視覚的に確認する作業が回帰診断です。
主に以下の4つのグラフを確認します。
Residuals vs Fitted（残差 vs 予測値）：線形性と等分散性の確認。平均 0 の周りにランダムに分散していれば OK。データに順番やIDがあって、誤差に周期性があってもよろしくないです。
Normal Q-Q：誤差項の正規性の確認。点が直線上に並んでいれば OK。
Scale-Location：等分散性の確認。標準化残差の絶対値の平方根をプロットし、ばらつきに偏りがないかを見ます。
Residuals vs Leverage：外れ値と影響力の確認。クックの距離（Cook's distance）が大きく、モデルに過度な影響を与えているデータ点がないかをチェックします。
回帰診断図。
Level 1 のまとめ
このレベルに達した分析者は、回帰分析を「ブラックボックスの関数呼び出し」ではなく、最小二乗法という最適化問題の解として、また n n n
OLS推定量を行列で導出し、その幾何学的意味（列空間への射影）を説明できるか？
t検定・F検定が「何を仮定して」「何を主張している」のかを区別できるか？
R 2 R^2 R 2 R 2 R^2 R 2
こうした道具立てを使いこなし、実務で統計的な分析が一通り回せるようになれば、Level 1 です。多くの統計学の教科書ではこの辺りまで書かれており、ある程度簡単な実務で統計的な処理をする分には問題がない水準と思います。
ただし、ここまでの理解は、データに対する仮定が満たされている場合に過ぎず、その仮定自体を問う段階にはまだ訪れておらず、先は長いです。
Level 2: ガウス=マルコフ定理と仮定の体系
ガウス=マルコフ定理の理論的背景と、その仮定が満たされなかった場合の対処法（処世術）を身につけていると Level 2 です。
ガウス=マルコフ定理
OLSが「良い推定量」である根拠は何か？その答えがガウス=マルコフ定理 です。
BLUEは平たくいうと、不偏（平均的に間違わない）であるシンプルな計算式（線形）の中で、一番ブレが少なく信用できる推定量ということです。手元のデータで線形回帰する場合、いつものOLSを使うのが最も精度が良いということを言っています。
仮定（A1〜A4）：
線形性 ：y = X β + ε \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} y = X β + ε y \mathbf{y} y X \mathbf{X} X β \boldsymbol{\beta} β β \boldsymbol{\beta} β x 2 x^2 x 2 log ( x ) \log(x) log ( x )
外生性（厳密外生性） ：E [ ε ∣ X ] = 0 \mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0} E [ ε ∣ X ] = 0 X \mathbf{X} X ε \boldsymbol{\varepsilon} ε ε \boldsymbol{\varepsilon} ε X \mathbf{X} X 欠落変数バイアス です。例えば、賃金を予測するモデルで「個人の本来の能力」をデータ化できず誤差項に押し付けた場合、能力と学歴（説明変数）には相関があるため、この外生性の仮定が崩壊してしまいます。
等分散性・無相関 ：Var ( ε ∣ X ) = σ 2 I \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \mathbf{I} Var ( ε ∣ X ) = σ 2 I σ 2 I \sigma^2 \mathbf{I} σ 2 I I \mathbf{I} I
等分散性 ：対角成分がすべて σ 2 \sigma^2 σ 2
無相関 ：非対角成分がすべて 0 0 0
X \mathbf{X} X X ⊤ X \mathbf{X}^\top\mathbf{X} X ⊤ X β ^ \hat{\boldsymbol{\beta}} β ^
BLUEであることの証明（導出）
なぜOLS推定量 β ^ = ( X ⊤ X ) − 1 X ⊤ y \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} β ^ = ( X ⊤ X ) − 1 X ⊤ y
任意の線形推定量を β ~ = C y \tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{y} β ~ = Cy C \mathbf{C} C C \mathbf{C} C C = ( X ⊤ X ) − 1 X ⊤ + D \mathbf{C} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top + \mathbf{D} C = ( X ⊤ X ) − 1 X ⊤ + D D \mathbf{D} D
1. 不偏性の条件 β ~ \tilde{\boldsymbol{\beta}} β ~ E [ β ~ ∣ X ] = β \mathbb{E}[\tilde{\boldsymbol{\beta}} \mid \mathbf{X}] = \boldsymbol{\beta} E [ β ~ ∣ X ] = β
E [ β ~ ∣ X ] = E [ C ( X β + ε ) ∣ X ] = C X β + C E [ ε ∣ X ] = C X β ( ∵ A2より )
\begin{aligned}
\mathbb{E}[\tilde{\boldsymbol{\beta}} \mid \mathbf{X}] &= \mathbb{E}[\mathbf{C}(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}) \mid \mathbf{X}] \\
&= \mathbf{C}\mathbf{X}\boldsymbol{\beta} + \mathbf{C}\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] \\
&= \mathbf{C}\mathbf{X}\boldsymbol{\beta} \quad (\because \text{A2より})
\end{aligned}
E [ β ~ ∣ X ] = E [ C ( X β + ε ) ∣ X ] = CX β + C E [ ε ∣ X ] = CX β ( ∵ A2 より ) これが任意の β \boldsymbol{\beta} β C X = I \mathbf{C}\mathbf{X} = \mathbf{I} CX = I C \mathbf{C} C
C X = ( ( X ⊤ X ) − 1 X ⊤ + D ) X = I + D X
\mathbf{C}\mathbf{X} = \left( (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top + \mathbf{D} \right)\mathbf{X} = \mathbf{I} + \mathbf{D}\mathbf{X}
CX = ( ( X ⊤ X ) − 1 X ⊤ + D ) X = I + DX したがって、β ~ \tilde{\boldsymbol{\beta}} β ~ D X = 0 \mathbf{D}\mathbf{X} = \mathbf{0} DX = 0
2. 分散の比較 β ~ \tilde{\boldsymbol{\beta}} β ~
Var ( β ~ ∣ X ) = Var ( C y ∣ X ) = C Var ( ε ∣ X ) C ⊤ = σ 2 C C ⊤ ( ∵ A3より )
\begin{aligned}
\text{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) &= \text{Var}(\mathbf{C}\mathbf{y} \mid \mathbf{X}) \\
&= \mathbf{C}\text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X})\mathbf{C}^\top \\
&= \sigma^2 \mathbf{C}\mathbf{C}^\top \quad (\because \text{A3より})
\end{aligned}
Var ( β ~ ∣ X ) = Var ( Cy ∣ X ) = C Var ( ε ∣ X ) C ⊤ = σ 2 C C ⊤ ( ∵ A3 より ) ここで C C ⊤ \mathbf{C}\mathbf{C}^\top C C ⊤
C C ⊤ = ( ( X ⊤ X ) − 1 X ⊤ + D ) ( X ( X ⊤ X ) − 1 + D ⊤ ) = ( X ⊤ X ) − 1 + ( X ⊤ X ) − 1 X ⊤ D ⊤ + D X ( X ⊤ X ) − 1 + D D ⊤
\begin{aligned}
\mathbf{C}\mathbf{C}^\top &= \left( (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top + \mathbf{D} \right) \left( \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1} + \mathbf{D}^\top \right) \\
&= (\mathbf{X}^\top\mathbf{X})^{-1} + (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{D}^\top + \mathbf{D}\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1} + \mathbf{D}\mathbf{D}^\top
\end{aligned}
C C ⊤ = ( ( X ⊤ X ) − 1 X ⊤ + D ) ( X ( X ⊤ X ) − 1 + D ⊤ ) = ( X ⊤ X ) − 1 + ( X ⊤ X ) − 1 X ⊤ D ⊤ + DX ( X ⊤ X ) − 1 + D D ⊤ 不偏性の条件 D X = 0 \mathbf{D}\mathbf{X} = \mathbf{0} DX = 0 X ⊤ D ⊤ = 0 \mathbf{X}^\top\mathbf{D}^\top = \mathbf{0} X ⊤ D ⊤ = 0
C C ⊤ = ( X ⊤ X ) − 1 + D D ⊤
\mathbf{C}\mathbf{C}^\top = (\mathbf{X}^\top\mathbf{X})^{-1} + \mathbf{D}\mathbf{D}^\top
C C ⊤ = ( X ⊤ X ) − 1 + D D ⊤ したがって、β ~ \tilde{\boldsymbol{\beta}} β ~
Var ( β ~ ∣ X ) = σ 2 ( X ⊤ X ) − 1 + σ 2 D D ⊤ = Var ( β ^ ∣ X ) + σ 2 D D ⊤
\text{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) = \sigma^2 (\mathbf{X}^\top\mathbf{X})^{-1} + \sigma^2 \mathbf{D}\mathbf{D}^\top = \text{Var}(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) + \sigma^2 \mathbf{D}\mathbf{D}^\top
Var ( β ~ ∣ X ) = σ 2 ( X ⊤ X ) − 1 + σ 2 D D ⊤ = Var ( β ^ ∣ X ) + σ 2 D D ⊤ 任意の行列 D \mathbf{D} D D D ⊤ \mathbf{D}\mathbf{D}^\top D D ⊤ Var ( β ~ ∣ X ) ⪰ Var ( β ^ ∣ X ) \text{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) \succeq \text{Var}(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) Var ( β ~ ∣ X ) ⪰ Var ( β ^ ∣ X ) β ^ \hat{\boldsymbol{\beta}} β ^
仮定に対する違反とその影響
Level 2の分析者には、各仮定が満たされなかった場合に、「どの性質が失われるのか」を明確に区別できるようになっていることが求められます。
仮定
満たされない場合に失われる性質・生じる問題
線形性（A1）
モデル自体の誤特定となり、推定された係数の解釈が不適切になる
外生性（A2）
推定量の不偏性および一致性が失われる （内生性バイアス）
等分散性（A3）
不偏性は保たれるが、効率性が失われ、標準誤差の推定が不適切になる
列フルランク（A4）
逆行列 ( X ⊤ X ) − 1 (\mathbf{X}^\top\mathbf{X})^{-1} ( X ⊤ X ) − 1
実務上、A2（外生性）の仮定に対する違反が多く、深刻化しやすいです 。これは点推定自体が真のパラメータから乖離してしまうためです（欠落変数バイアスや同時決定などの問題）。
1. 線形性の仮定が崩れた場合（モデルの誤特定）
データの真の姿が非線形であるにもかかわらず、単純な直線を当てはめてしまうと、係数の解釈が歪み、予測精度も低下します。この場合、OLSの枠組み（パラメータについて線形）を維持しつつ、変数の表現を工夫するか、あるいは別のアルゴリズムに頼るかの選択になります。
変数変換と多項式の追加
実務で最も頻繁に行われるアプローチです。散布図や残差プロットを確認し、関係性が直線でない場合は説明変数を変換してモデルに投入します。
対数変換（log ( x ) \log(x) log ( x ) : 経済データ（賃金、売上、人口など）で頻出します。「x x x y y y
多項式の追加（x 2 , x 3 x^2, x^3 x 2 , x 3 : 「年齢と賃金（中年でピーク）」「気温と電力消費（夏と冬にピーク）」のように、U字型や逆U字型の関係が想定される場合、説明変数に x x x x 2 x^2 x 2 β \boldsymbol{\beta} β
「自分の作ったモデルが適切に線形性を捉えられているか（重要な非線形項を見落としていないか）」を統計的にチェックする処世術です。
まず通常のモデルで回帰し、予測値 y ^ \hat{\mathbf{y}} y ^
次に、その予測値の2乗や3乗（y ^ 2 , y ^ 3 \hat{\mathbf{y}}^2, \hat{\mathbf{y}}^3 y ^ 2 , y ^ 3
もし追加した項の係数が統計的に有意であれば、「モデルにはまだ表現しきれていない非線形な構造が残っている（誤特定がある）」と判断し、変数変換などを再検討します。
機械学習手法への切り替え
分析の目的が「係数から因果関係を解釈すること」ではなく、「未知のデータを高精度で予測すること」である場合の手段です。
2. 外生性の仮定が崩れた場合（内生性問題）
説明変数と誤差項が相関してしまう内生性は、推定量の不偏性・一致性を破綻させる最も深刻な問題です。これを解決するには、相関の根源（交絡や逆因果）を断ち切るアプローチが求められます。
操作変数法（IV）と2段階最小二乗法（2SLS）
内生性を持つ説明変数（X X X
これには、以下の2条件を満たす 操作変数（Z Z Z を見つけてくる必要があります。
関連性 : 操作変数 Z Z Z X X X
除外制約 : 操作変数 Z Z Z ε \boldsymbol{\varepsilon} ε X X X Y Y Y
【2SLSの手順】
第1段階 : X X X Z Z Z X X X X ^ \hat{X} X ^ X ^ \hat{X} X ^ Z Z Z
第2段階 : 本来の目的変数 Y Y Y X X X X ^ \hat{X} X ^ β \boldsymbol{\beta} β
パネルデータ分析（固定効果モデル）
同一の個人や企業を、複数期間にわたって観測したデータ（パネルデータ）がある場合に使える処世術です。観測できないが時間的に変化しない要因 （個体固有効果）によって引き起こされます。
3. 等分散性・無相関の仮定が崩れた場合への処世術
この仮定が崩れても、OLS推定量の係数 β ^ \hat{\boldsymbol{\beta}} β ^ BLUE（最小分散）ではなくなる という点です。これを効率性の喪失と呼びます。誤差項に不均一分散がある場合、誤差が小さいサンプルの情報には重みを置き、誤差が大きいサンプルの情報は軽めに扱った方が、よりブレの少ない（分散の小さい）推定が可能になります。OLSはすべてのデータを一律に扱うため、この情報の最適化ができず、結果として推定値が「たまたまの外れ値」に振り回されやすくなってしまうのです。また、この効率性の喪失だけでなく、「標準誤差の計算式」が狂うため、t検定の結果なども信頼できなくなります。
不均一分散下での効率性のロス
実際に、等分散性の仮定が崩れ、誤差項の分散共分散行列が Var ( ε ∣ X ) = Ω \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \boldsymbol{\Omega} Var ( ε ∣ X ) = Ω Ω ≠ σ 2 I \boldsymbol{\Omega} \neq \sigma^2 \mathbf{I} Ω = σ 2 I
このとき、通常のOLS推定量 β ^ OLS = ( X ⊤ X ) − 1 X ⊤ y \hat{\boldsymbol{\beta}}_{\text{OLS}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} β ^ OLS = ( X ⊤ X ) − 1 X ⊤ y
Var ( β ^ OLS ∣ X ) = ( X ⊤ X ) − 1 X ⊤ Ω X ( X ⊤ X ) − 1
\text{Var}(\hat{\boldsymbol{\beta}}_{\text{OLS}} \mid \mathbf{X}) = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol{\Omega}\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}
Var ( β ^ OLS ∣ X ) = ( X ⊤ X ) − 1 X ⊤ Ω X ( X ⊤ X ) − 1 一方、この不均一分散の構造 Ω \boldsymbol{\Omega} Ω 一般化最小二乗法（GLS） です。GLS推定量 β ^ GLS \hat{\boldsymbol{\beta}}_{\text{GLS}} β ^ GLS
Var ( β ^ GLS ∣ X ) = ( X ⊤ Ω − 1 X ) − 1
\text{Var}(\hat{\boldsymbol{\beta}}_{\text{GLS}} \mid \mathbf{X}) = (\mathbf{X}^\top\boldsymbol{\Omega}^{-1}\mathbf{X})^{-1}
Var ( β ^ GLS ∣ X ) = ( X ⊤ Ω − 1 X ) − 1 エイトケンの定理（ガウス=マルコフ定理の一般化）により、任意の線形推定量の中でGLSが最小分散を持つため、両者の分散共分散行列の間には必ず以下の大小関係が成立します。
Var ( β ^ OLS ∣ X ) − Var ( β ^ GLS ∣ X ) ⪰ 0 （半正定値行列）
\text{Var}(\hat{\boldsymbol{\beta}}_{\text{OLS}} \mid \mathbf{X}) - \text{Var}(\hat{\boldsymbol{\beta}}_{\text{GLS}} \mid \mathbf{X}) \succeq 0 \quad \text{（半正定値行列）}
Var ( β ^ OLS ∣ X ) − Var ( β ^ GLS ∣ X ) ⪰ 0 （半正定値行列）
どれくらい効率が悪くなるのか？
この分散の差（＝無駄にしてしまった効率性）の大きさは、説明変数 X \mathbf{X} X
例として「企業の規模（X X X ω i 2 \omega_i^2 ω i 2 Ω − 1 \boldsymbol{\Omega}^{-1} Ω − 1
ロバスト標準誤差（不均一分散に頑健な標準誤差）
実務分析においてよく使う方法です。OLSのデフォルトの分散共分散行列の計算式 σ 2 ( X ⊤ X ) − 1 \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1} σ 2 ( X ⊤ X ) − 1 σ 2 \sigma^2 σ 2 ε ^ i 2 \hat{\varepsilon}_i^2 ε ^ i 2 Pythonのstatsmodels など）では、引数に cov_type='HC3' などを指定するだけで計算してくれます。
クラスター化ロバスト標準誤差
データの中に「同じクラスの生徒」「同じ都道府県の企業」「同一人物の複数年のデータ（パネルデータ）」といったグループ（クラスター）構造がある場合、同一クラスター内のサンプル同士は似たような未観測のショックを共有する ため、誤差項に相関（空間的・時間的な自己相関）が生じます。
これを無視して通常の標準誤差を計算すると、実質的な情報量を過大評価してしまいます。たとえば「同じクラスの生徒40人」のデータは、完全に独立した40人分の情報を持っているわけではありません。これを独立だと仮定して分析すると、標準誤差が不当に小さく（t値が不当に大きく）計算され、本当は有意でないのに有意と判定してしまう危険性 が非常に高くなります。これをモウルトン・バイアス （Moulton bias）と呼びます。
クラスター構造とサンドイッチ推定量
全データ数を N N N G G G g g g N g N_g N g g = 1 , … , G g = 1, \dots, G g = 1 , … , G
y g = X g β + ε g
\mathbf{y}_g = \mathbf{X}_g\boldsymbol{\beta} + \boldsymbol{\varepsilon}_g
y g = X g β + ε g 通常のOLSはすべての誤差項が無相関だと仮定しますが、クラスター構造がある場合、異なるクラスター間の誤差は無相関（独立）ですが、同じクラスター内の誤差には相関がある と考えます。つまり、（データをクラスターで固めて並び替えれば）誤差項の全体の分散共分散行列 Ω \boldsymbol{\Omega} Ω
Ω = ( Ω 1 0 … 0 0 Ω 2 … 0 ⋮ ⋮ ⋱ ⋮ 0 0 … Ω G )
\boldsymbol{\Omega} = \begin{pmatrix} \boldsymbol{\Omega}_1 & \mathbf{0} & \dots & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\Omega}_2 & \dots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \dots & \boldsymbol{\Omega}_G \end{pmatrix}
Ω = Ω 1 0 ⋮ 0 0 Ω 2 ⋮ 0 … … ⋱ … 0 0 ⋮ Ω G この構造を許容した上で、正しい分散を推定するのがクラスター化ロバスト分散推定量（CRVE: Cluster-Robust Variance Estimator） です。通常のロバスト標準誤差（Huber-White推定量）の考え方を拡張し、個別の残差の2乗の代わりに、クラスターごとの残差ベクトル ε ^ g \hat{\boldsymbol{\varepsilon}}_g ε ^ g を用いてサンドイッチの「具」の部分を構成します。
Var ^ ( β ^ cluster ) = c ( X ⊤ X ) − 1 ( ∑ g = 1 G X g ⊤ ε ^ g ε ^ g ⊤ X g ) ( X ⊤ X ) − 1
\widehat{\text{Var}}(\hat{\boldsymbol{\beta}}_{\text{cluster}}) = c (\mathbf{X}^\top\mathbf{X})^{-1} \left( \sum_{g=1}^G \mathbf{X}_g^\top \hat{\boldsymbol{\varepsilon}}_g \hat{\boldsymbol{\varepsilon}}_g^\top \mathbf{X}_g \right) (\mathbf{X}^\top\mathbf{X})^{-1}
Var ( β ^ cluster ) = c ( X ⊤ X ) − 1 ( g = 1 ∑ G X g ⊤ ε ^ g ε ^ g ⊤ X g ) ( X ⊤ X ) − 1 この式が意味するのは、クラスター内の任意の相関（不均一分散を含む）を許容しつつ、漸近的に正しい標準誤差を計算できるということです。ただし、この大数の法則が機能するためにはクラスターの数 G G G
【実務での処世術・実装方法】
実務の分析では非常によく使います。数式は複雑ですが、実装はオプション一つで解決してしまいます。
Python（statsmodels）の場合 fit() メソッドの引数に cov_type='cluster' を指定し、どの変数がクラスターを示すのかを渡します。
statsmodelsでのクラスター化ロバスト標準誤差
Rの場合 estimatr パッケージの lm_robust() を使うのが最も簡単です。
lm_robustによるクラスター化ロバスト標準誤差
どちらの言語でも、出力される係数の推定値（coef）自体は通常のOLSと全く同じ値になります。しかし、標準誤差（std err）がクラスター構造を加味して適切に大きめに補正され、それに伴いt値とp値がより慎重な値に変化 していることが確認できるはずです。
一般化最小二乗法（GLS / FGLS）
ロバスト標準誤差が「OLSの計算結果をそのまま残し、標準誤差の計算式だけを後から補正する」という対症療法であるのに対し、GLS（一般化最小二乗法）は「データを根本から変換することで、OLSが再びBLUE（最小分散）となる性質を取り戻す」というアプローチです。
なぜデータを変換するのか？
不均一分散や自己相関が存在する場合、誤差項の分散共分散行列は Var ( ε ∣ X ) = Ω \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \boldsymbol{\Omega} Var ( ε ∣ X ) = Ω Ω \boldsymbol{\Omega} Ω
ここで、行列の性質から Ω = P P ⊤ \boldsymbol{\Omega} = \mathbf{P}\mathbf{P}^\top Ω = P P ⊤ P \mathbf{P} P P \mathbf{P} P Ω − 1 / 2 \boldsymbol{\Omega}^{-1/2} Ω − 1/2 Ω − 1 / 2 Ω ( Ω − 1 / 2 ) ⊤ = I \boldsymbol{\Omega}^{-1/2}\boldsymbol{\Omega}(\boldsymbol{\Omega}^{-1/2})^\top = \mathbf{I} Ω − 1/2 Ω ( Ω − 1/2 ) ⊤ = I
元の回帰モデル y = X β + ε \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} y = X β + ε Ω − 1 / 2 \boldsymbol{\Omega}^{-1/2} Ω − 1/2
Ω − 1 / 2 y = Ω − 1 / 2 X β + Ω − 1 / 2 ε
\boldsymbol{\Omega}^{-1/2}\mathbf{y} = \boldsymbol{\Omega}^{-1/2}\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\Omega}^{-1/2}\boldsymbol{\varepsilon}
Ω − 1/2 y = Ω − 1/2 X β + Ω − 1/2 ε これを新しい変数 y ~ = X ~ β + ε ~ \tilde{\mathbf{y}} = \tilde{\mathbf{X}}\boldsymbol{\beta} + \tilde{\boldsymbol{\varepsilon}} y ~ = X ~ β + ε ~ ε ~ \tilde{\boldsymbol{\varepsilon}} ε ~
Var ( ε ~ ∣ X ) = Var ( Ω − 1 / 2 ε ∣ X ) = Ω − 1 / 2 Var ( ε ∣ X ) ( Ω − 1 / 2 ) ⊤ = Ω − 1 / 2 Ω ( Ω − 1 / 2 ) ⊤ = I
\begin{aligned}
\text{Var}(\tilde{\boldsymbol{\varepsilon}} \mid \mathbf{X}) &= \text{Var}(\boldsymbol{\Omega}^{-1/2}\boldsymbol{\varepsilon} \mid \mathbf{X}) \\
&= \boldsymbol{\Omega}^{-1/2} \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) (\boldsymbol{\Omega}^{-1/2})^\top \\
&= \boldsymbol{\Omega}^{-1/2} \boldsymbol{\Omega} (\boldsymbol{\Omega}^{-1/2})^\top \\
&= \mathbf{I}
\end{aligned}
Var ( ε ~ ∣ X ) = Var ( Ω − 1/2 ε ∣ X ) = Ω − 1/2 Var ( ε ∣ X ) ( Ω − 1/2 ) ⊤ = Ω − 1/2 Ω ( Ω − 1/2 ) ⊤ = I 分散共分散行列が単位行列 I \mathbf{I} I 変換後のデータにおいては、完全に等分散かつ無相関の仮定（A3）が満たされている ことを意味します。
この推定量 β ^ GLS \hat{\boldsymbol{\beta}}_{\text{GLS}} β ^ GLS
β ^ GLS = ( X ~ ⊤ X ~ ) − 1 X ~ ⊤ y ~ = ( ( Ω − 1 / 2 X ) ⊤ ( Ω − 1 / 2 X ) ) − 1 ( Ω − 1 / 2 X ) ⊤ ( Ω − 1 / 2 y ) = ( X ⊤ Ω − 1 X ) − 1 X ⊤ Ω − 1 y
\begin{aligned}
\hat{\boldsymbol{\beta}}_{\text{GLS}} &= (\tilde{\mathbf{X}}^\top\tilde{\mathbf{X}})^{-1}\tilde{\mathbf{X}}^\top\tilde{\mathbf{y}} \\
&= ((\boldsymbol{\Omega}^{-1/2}\mathbf{X})^\top(\boldsymbol{\Omega}^{-1/2}\mathbf{X}))^{-1}(\boldsymbol{\Omega}^{-1/2}\mathbf{X})^\top(\boldsymbol{\Omega}^{-1/2}\mathbf{y}) \\
&= (\mathbf{X}^\top\boldsymbol{\Omega}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol{\Omega}^{-1}\mathbf{y}
\end{aligned}
β ^ GLS = ( X ~ ⊤ X ~ ) − 1 X ~ ⊤ y ~ = (( Ω − 1/2 X ) ⊤ ( Ω − 1/2 X ) ) − 1 ( Ω − 1/2 X ) ⊤ ( Ω − 1/2 y ) = ( X ⊤ Ω − 1 X ) − 1 X ⊤ Ω − 1 y この式の直観的な理解としては、Ω − 1 \boldsymbol{\Omega}^{-1} Ω − 1
実行可能一般化最小二乗法（FGLS）による実務的アプローチ
理論上は完璧なGLSですが、実務においては真の分散共分散行列 Ω \boldsymbol{\Omega} Ω N N N Ω \boldsymbol{\Omega} Ω N ( N + 1 ) / 2 N(N+1)/2 N ( N + 1 ) /2
そこで実務では、Ω \boldsymbol{\Omega} Ω Ω \boldsymbol{\Omega} Ω Ω ^ \hat{\boldsymbol{\Omega}} Ω ^ 実行可能一般化最小二乗法（FGLS: Feasible GLS） が用いられます。
【FGLSの一般的なステップ】
とりあえずOLSを実行 : 通常のOLSで y \mathbf{y} y X \mathbf{X} X ε ^ \hat{\boldsymbol{\varepsilon}} ε ^
分散の構造をモデリング : 残差の2乗 ε ^ i 2 \hat{\varepsilon}_i^2 ε ^ i 2 x 1 i x_{1i} x 1 i log ( ε ^ i 2 ) = α 0 + α 1 x 1 i + v i \log(\hat{\varepsilon}_i^2) = \alpha_0 + \alpha_1 x_{1i} + v_i log ( ε ^ i 2 ) = α 0 + α 1 x 1 i + v i
重みの作成 : 補助回帰の予測値から、各サンプルの分散の推定値 ω ^ i 2 \hat{\omega}_i^2 ω ^ i 2 Ω ^ \hat{\boldsymbol{\Omega}} Ω ^
FGLSの実行 : 推定した Ω ^ \hat{\boldsymbol{\Omega}} Ω ^ β ^ FGLS = ( X ⊤ Ω ^ − 1 X ) − 1 X ⊤ Ω ^ − 1 y \hat{\boldsymbol{\beta}}_{\text{FGLS}} = (\mathbf{X}^\top\hat{\boldsymbol{\Omega}}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\hat{\boldsymbol{\Omega}}^{-1}\mathbf{y} β ^ FGLS = ( X ⊤ Ω ^ − 1 X ) − 1 X ⊤ Ω ^ − 1 y
大標本においては、分散構造のモデルが正しい場合、Ω ^ \hat{\boldsymbol{\Omega}} Ω ^ Ω \boldsymbol{\Omega} Ω
4. 列フルランクの仮定が崩れた場合（多重共線性）への処世術
完全な多重共線性は、データに由来する避けられない問題というよりも、分析者によるミスです。
冗長なダミー変数
完全な多重共線性の代表例として、ダミー変数の種類の数だけダミー変数を用意してしまうものがあります。
これの対処は簡単で、カテゴリが N N N N − 1 N-1 N − 1
※補足：「完全ではないが、強い」多重共線性への処世術
実務でより頭を悩ませるのは、相関係数が0.9を超えるような「強いが完全ではない多重共線性」です。計算自体は実行できてしまうものの、推定量 β ^ \hat{\boldsymbol{\beta}} β ^
VIF（分散拡大係数）の確認 : 各説明変数が他の説明変数からどれくらい説明されてしまうか（多重共線性の深刻度）を測る指標です。一般にVIFが10を超える変数は、多重共線性の原因として警戒されます。
変数の統合・削除 : 「体重」と「BMI」のように、意味的に重複している変数はどちらか一方を削除するか、主成分分析（PCA）などで1つの指標に統合します。
正則化手法（Ridge回帰・Lasso回帰） : 予測を目的とする場合、ペナルティ項を設けることで係数が極端な値をとることを防ぐリッジ回帰などが極めて有効です。あえて推定量にわずかな「バイアス（偏り）」を入れる代償として、分散を劇的に縮小させ、安定したモデルを構築する処世術です。
Level 2のまとめ
このレベルに達した分析者は、OLSがBLUEであることの恩恵と脆さを同時に理解しています。ガウス＝マルコフの仮定は、回帰分析が成立するための土台ですが、どれが崩れたときに何が失われるのかを区別できることが重要です。
仮定A1〜A4のうち、どれが崩れると「不偏性」が失われ、どれが崩れると「効率性」だけが失われるのかを即答できるか？
内生性、不均一分散、自己相関、多重共線性に対して、それぞれ操作変数法・ロバスト標準誤差・クラスター化標準誤差・GLS/FGLSといった処世術を、状況に応じて選択できるか？
「ロバスト標準誤差をつけておけば安心」ではなく、なぜそれが必要なのかをサンドイッチ分散の式から説明できるか？
仮定の意味とその違反への対処を体系的に語れるようになれば、Level 2です。しかしここまでは「不偏推定量の中で最良を目指す」範疇であり、あえてバイアスを受け入れるという発想はまだ登場していません。
Level 3: 予測と汎化への拡張
Level 2までの議論は、不偏推定量（E [ β ^ ] = β \mathbb{E}[\hat{\beta}] = \beta E [ β ^ ] = β
Level 3では、回帰の本質を「条件付き期待値の推定」と再定義し、バイアスとバリアンスのトレードオフを能動的に制御する手法群を扱います。
回帰の上位概念
統計的学習の枠組みにおいて、回帰問題は「入力 X X X Y Y Y f ( X ) f(X) f ( X )
f ( x ) = arg min g E [ ( Y − g ( X ) ) 2 ∣ X = x ]
f(x) = \arg\min_{g} \mathbb{E}[(Y - g(X))^2 \mid X = x]
f ( x ) = arg g min E [( Y − g ( X ) ) 2 ∣ X = x ] この解は、理論的には条件付き期待値 f ( x ) = E [ Y ∣ X = x ] f(x) = \mathbb{E}[Y \mid X = x] f ( x ) = E [ Y ∣ X = x ]
E [ ( Y − g ( X ) ) 2 ∣ X = x ] = E [ ( Y − E [ Y ∣ X ] + E [ Y ∣ X ] − g ( X ) ) 2 ∣ X = x ] = E [ ( Y − E [ Y ∣ X ] ) 2 ∣ X = x ] + ( E [ Y ∣ X ] − g ( X ) ) 2
\begin{aligned}
\mathbb{E}[(Y - g(X))^2 \mid X = x] &= \mathbb{E}[(Y - \mathbb{E}[Y|X] + \mathbb{E}[Y|X] - g(X))^2 \mid X = x] \\
&= \mathbb{E}[(Y - \mathbb{E}[Y|X])^2 \mid X = x] + (\mathbb{E}[Y|X] - g(X))^2
\end{aligned}
E [( Y − g ( X ) ) 2 ∣ X = x ] = E [( Y − E [ Y ∣ X ] + E [ Y ∣ X ] − g ( X ) ) 2 ∣ X = x ] = E [( Y − E [ Y ∣ X ] ) 2 ∣ X = x ] + ( E [ Y ∣ X ] − g ( X ) ) 2 第1項はデータの本質的なノイズであり、g ( X ) g(X) g ( X ) g ( x ) = E [ Y ∣ X = x ] g(x) = \mathbb{E}[Y \mid X = x] g ( x ) = E [ Y ∣ X = x ] f ( X ) f(X) f ( X ) X β \mathbf{X}\boldsymbol{\beta} X β
正則化：Ridge、Lasso
高次元データ（p ≈ n p \approx n p ≈ n X ⊤ X \mathbf{X}^\top\mathbf{X} X ⊤ X β ^ = ( X ⊤ X ) − 1 X ⊤ y \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} β ^ = ( X ⊤ X ) − 1 X ⊤ y
正則化（Regularization）は、目的関数に残差平方和（RSS）以外のペナルティ項を加えることで、この分散を抑制します。
Ridge回帰（L2正則化）
係数の二乗和をペナルティとして加えます。
β ^ Ridge = arg min β { ∥ y − X β ∥ 2 + λ ∥ β ∥ 2 2 }
\hat{\boldsymbol{\beta}}_{\text{Ridge}} = \arg\min_{\boldsymbol{\beta}} \left\{ \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \lambda \|\boldsymbol{\beta}\|_2^2 \right\}
β ^ Ridge = arg β min { ∥ y − X β ∥ 2 + λ ∥ β ∥ 2 2 } 正規方程式を解くと、以下の解が導かれます。
β ^ Ridge = ( X ⊤ X + λ I ) − 1 X ⊤ y
\hat{\boldsymbol{\beta}}_{\text{Ridge}} = (\mathbf{X}^\top\mathbf{X} + \lambda \mathbf{I})^{-1}\mathbf{X}^\top \mathbf{y}
β ^ Ridge = ( X ⊤ X + λ I ) − 1 X ⊤ y ここで λ I \lambda \mathbf{I} λ I X ⊤ X \mathbf{X}^\top\mathbf{X} X ⊤ X
Lasso回帰（L1正則化）
係数の絶対値の和をペナルティとして加えます。
β ^ Lasso = arg min β { ∥ y − X β ∥ 2 + λ ∥ β ∥ 1 }
\hat{\boldsymbol{\beta}}_{\text{Lasso}} = \arg\min_{\boldsymbol{\beta}} \left\{ \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \lambda \|\boldsymbol{\beta}\|_1 \right\}
β ^ Lasso = arg β min { ∥ y − X β ∥ 2 + λ ∥ β ∥ 1 } L1ノルムの幾何学的性質により、多くの係数が正確に 0 になる「スパース性」が生じます。これは予測精度の向上と同時に、自動的な変数選択としての機能を果たします。なお、Ridgeと異なり解析的に解けないので数値的に解くことになります。
正則化の解釈：ベイズ推定との関連
正則化はベイズ統計学における最大事後確率（MAP）推定 として解釈できます。
Ridge : 係数 β \boldsymbol{\beta} β
Lasso : 係数 β \boldsymbol{\beta} β
線形回帰 y i = x i ⊤ β + ε i y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i y i = x i ⊤ β + ε i ε i ∼ i.i.d. N ( 0 , σ 2 ) \varepsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2) ε i ∼ i.i.d. N ( 0 , σ 2 )
p ( β ∣ y , X ) ∝ p ( y ∣ X , β ) p ( β )
p(\boldsymbol{\beta} \mid \boldsymbol{y}, X) \propto p(\boldsymbol{y} \mid X, \boldsymbol{\beta}) \, p(\boldsymbol{\beta})
p ( β ∣ y , X ) ∝ p ( y ∣ X , β ) p ( β ) MAP 推定は事後分布を最大化する点なので、対数をとって、
β ^ MAP = arg max β [ log p ( y ∣ X , β ) + log p ( β ) ]
\hat{\boldsymbol{\beta}}_{\text{MAP}} = \arg\max_{\boldsymbol{\beta}} \left[ \log p(\boldsymbol{y} \mid X, \boldsymbol{\beta}) + \log p(\boldsymbol{\beta}) \right]
β ^ MAP = arg β max [ log p ( y ∣ X , β ) + log p ( β ) ] 尤度項はガウス誤差の仮定から、
log p ( y ∣ X , β ) = − 1 2 σ 2 ∥ y − X β ∥ 2 2 + const
\log p(\boldsymbol{y} \mid X, \boldsymbol{\beta}) = -\frac{1}{2\sigma^2} \| \boldsymbol{y} - X\boldsymbol{\beta} \|_2^2 + \text{const}
log p ( y ∣ X , β ) = − 2 σ 2 1 ∥ y − X β ∥ 2 2 + const したがって MAP 推定は次の最小化問題と等価になります。
β ^ MAP = arg min β [ 1 2 σ 2 ∥ y − X β ∥ 2 2 − log p ( β ) ]
\hat{\boldsymbol{\beta}}_{\text{MAP}} = \arg\min_{\boldsymbol{\beta}} \left[ \frac{1}{2\sigma^2} \| \boldsymbol{y} - X\boldsymbol{\beta} \|_2^2 - \log p(\boldsymbol{\beta}) \right]
β ^ MAP = arg β min [ 2 σ 2 1 ∥ y − X β ∥ 2 2 − log p ( β ) ] Ridge の場合 ：β j ∼ i.i.d. N ( 0 , τ 2 ) \beta_j \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \tau^2) β j ∼ i.i.d. N ( 0 , τ 2 )
− log p ( β ) = 1 2 τ 2 ∥ β ∥ 2 2 + const
-\log p(\boldsymbol{\beta}) = \frac{1}{2\tau^2} \| \boldsymbol{\beta} \|_2^2 + \text{const}
− log p ( β ) = 2 τ 2 1 ∥ β ∥ 2 2 + const となり、λ = σ 2 / τ 2 \lambda = \sigma^2 / \tau^2 λ = σ 2 / τ 2
∥ y − X β ∥ 2 2 + λ ∥ β ∥ 2 2
\| \boldsymbol{y} - X\boldsymbol{\beta} \|_2^2 + \lambda \| \boldsymbol{\beta} \|_2^2
∥ y − X β ∥ 2 2 + λ ∥ β ∥ 2 2 と一致します。
Lasso の場合 ：β j ∼ i.i.d. Laplace ( 0 , b ) \beta_j \overset{\text{i.i.d.}}{\sim} \text{Laplace}(0, b) β j ∼ i.i.d. Laplace ( 0 , b ) p ( β j ) = 1 2 b exp ( − ∣ β j ∣ / b ) p(\beta_j) = \frac{1}{2b} \exp(-|\beta_j| / b) p ( β j ) = 2 b 1 exp ( − ∣ β j ∣/ b )
− log p ( β ) = 1 b ∥ β ∥ 1 + const
-\log p(\boldsymbol{\beta}) = \frac{1}{b} \| \boldsymbol{\beta} \|_1 + \text{const}
− log p ( β ) = b 1 ∥ β ∥ 1 + const となり、λ = 2 σ 2 / b \lambda = 2\sigma^2 / b λ = 2 σ 2 / b
∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1
\| \boldsymbol{y} - X\boldsymbol{\beta} \|_2^2 + \lambda \| \boldsymbol{\beta} \|_1
∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 と一致します。
事前分布の裾の重さ と原点での尖り が、推定量の性質を決めています。ガウス事前は原点で滑らかなので係数を 0 方向に「縮小」するだけですが、ラプラス事前は原点で微分不可能な尖りを持つため、解が厳密に β j = 0 \beta_j = 0 β j = 0 スパース性 が生じます。
バイアス-バリアンス分解
正則化がなぜ有効かを理解するための概念が、バイアス-バリアンス分解です。特定の点 x 0 x_0 x 0
Err ( x 0 ) = E [ ( Y − f ^ ( x 0 ) ) 2 ] = σ 2 ⏟ 不可避な誤差 + ( E [ f ^ ( x 0 ) ] − f ( x 0 ) ) 2 ⏟ バイアス 2 + E [ ( f ^ ( x 0 ) − E [ f ^ ( x 0 ) ] ) 2 ] ⏟ バリアンス
\text{Err}(x_0) = \mathbb{E}[(Y - \hat{f}(x_0))^2] = \underbrace{\sigma^2}_{\text{不可避な誤差}} + \underbrace{\left( \mathbb{E}[\hat{f}(x_0)] - f(x_0) \right)^2}_{\text{バイアス}^2} + \underbrace{\mathbb{E}[(\hat{f}(x_0) - \mathbb{E}[\hat{f}(x_0)])^2]}_{\text{バリアンス}}
Err ( x 0 ) = E [( Y − f ^ ( x 0 ) ) 2 ] = 不可避な誤差 σ 2 + バイアス 2 ( E [ f ^ ( x 0 )] − f ( x 0 ) ) 2 + バリアンス E [( f ^ ( x 0 ) − E [ f ^ ( x 0 )] ) 2 ] OLSは不偏推定量であるため、バイアスは 0 ですが、モデルが複雑であったりデータが少なかったりするとバリアンスが非常に大きくなります。λ \lambda λ 意図的にバイアスを導入し、それ以上にバリアンスを削減する 手法です。その結果、全体のMSE（Bias 2 + Variance \text{Bias}^2 + \text{Variance} Bias 2 + Variance
非線形への一般化
線形モデルの枠組みを保ちつつ、データの非線形な関係を捉えるための拡張手法です。
多項式回帰とスプライン回帰
説明変数を ϕ ( x ) = ( x , x 2 , … , x k ) \phi(x) = (x, x^2, \ldots, x^k) ϕ ( x ) = ( x , x 2 , … , x k ) スプライン回帰 は領域を複数の区間に分割し、接続点（ノット）で連続かつ滑らか（一般に2次微分まで連続）になるよう低次多項式を当てはめます。
多項式回帰と自然キュービックスプラインの比較。高次の多項式回帰（赤）は学習データの範囲外で予測が極端に暴れる（ルンゲ現象 ）のに対し、自然スプライン（緑）はノット間で滑らかに結合しつつ、境界外では線形に保たれるためバリアンスが抑えられています。
特にB-スプライン や自然スプライン は、計算の安定性と解釈性のバランスに優れており、非線形なトレンドを捉えるための標準的な手法です。
一般化加法モデル（GAM）
複数の説明変数がある場合、各変数ごとの非線形効果を足し合わせた形でモデル化します。
y = β 0 + ∑ j = 1 p f j ( x j ) + ε
y = \beta_0 + \sum_{j=1}^p f_j(x_j) + \varepsilon
y = β 0 + j = 1 ∑ p f j ( x j ) + ε 各 f j f_j f j
一般化加法モデル（GAM）による成分プロットの概念図。各変数の非線形な影響（U字型、対数型など）を個別の関数 f j ( x j ) f_j(x_j) f j ( x j )
カーネル回帰
特定の関数形を仮定せず、予測したい点 x x x ( x i , y i ) (x_i, y_i) ( x i , y i ) Nadaraya-Watson推定量 は次のように定義されます。
f ^ ( x ) = ∑ i = 1 n K h ( x − x i ) y i ∑ i = 1 n K h ( x − x i )
\hat{f}(x) = \frac{\sum_{i=1}^n K_h(x - x_i) y_i}{\sum_{i=1}^n K_h(x - x_i)}
f ^ ( x ) = ∑ i = 1 n K h ( x − x i ) ∑ i = 1 n K h ( x − x i ) y i ここで K h K_h K h x x x E [ Y ∣ X ] \mathbb{E}[Y|X] E [ Y ∣ X ]
Level 3のまとめ
このレベルに達した分析者は、不偏性に捉われず、予測の高性能化を目指せます。回帰の目的が「未知のデータをうまく予測する」ことに力点を置き、あえてバイアスを導入してバリアンスを削るというトレードオフが視野に入ってきます。
回帰の上位概念が「条件付き期待値 E [ Y ∣ X ] \mathbb{E}[Y \mid X] E [ Y ∣ X ]
Ridge・Lassoを単なる「過学習対策」ではなく、ベイズ的にはガウス事前・ラプラス事前を置いたMAP推定として説明できるか？
バイアス-バリアンス分解の各項を式で書き下し、λ \lambda λ
スプライン・GAM・カーネル回帰など、線形性の枠を保ったまま非線形性を取り込む手法を、目的に応じて使い分けられるか？
（自由度調整済み）決定係数 R 2 R^2 R 2
さて、ここまでの議論は「相関」の話であり、X X X Y Y Y
Level 4: 因果推論へ
Level 3までで、私たちは「手元のデータに対して、予測するための妥当な線を引く方法」を習得しました。次の壁は、「回帰係数が有意であること、予測性能が高いこと」と「変数 X X X Y Y Y
Level4では、予測モデルとしての回帰を、因果推論のためのツール として見ていきます。
観測データの限界と内生性
私たちが手にするデータの多くは、ランダム化比較試験（RCT）によって得られたものではなく、単なる観測データです。観測データにおいて、Level2で学んだガウス＝マルコフの仮定 A2（外生性：E [ ε ∣ X ] = 0 \mathbb{E}[\varepsilon \mid \mathbf{X}] = 0 E [ ε ∣ X ] = 0
この仮定が崩れた状態を 内生性（Endogeneity） があると言います。内生性があるとき、OLS推定量は不偏性も一致性も失い、因果推論が適切にできなくなります。
内生性を引き起こす3パターン
欠落変数バイアス（OVB） : X X X Y Y Y
同時決定（逆の因果） : X X X Y Y Y Y Y Y X X X
測定誤差 : 説明変数 X X X
欠落変数バイアス（OVB）
まずは説明変数が欠落した場合（必要な変数を入力できていない場合）、推定した回帰係数が信頼できなくなることを確認します。
真のモデルが以下の2変数モデルであるとします。
y = β 1 x 1 + β 2 x 2 + ε
y = \beta_1 x_1 + \beta_2 x_2 + \varepsilon
y = β 1 x 1 + β 2 x 2 + ε ここで、分析者が x 2 x_2 x 2 x 1 x_1 x 1
y = α 1 x 1 + u
y = \alpha_1 x_1 + u
y = α 1 x 1 + u このとき、推定される係数 α ^ 1 \hat{\alpha}_1 α ^ 1
α ^ 1 = Cov ( x 1 , y ) Var ( x 1 ) \hat{\alpha}_1 = \frac{\text{Cov}(x_1, y)}{\text{Var}(x_1)} α ^ 1 = Var ( x 1 ) Cov ( x 1 , y ) ここで、手元のデータ y y y y = β 1 x 1 + β 2 x 2 + ε y = \beta_1 x_1 + \beta_2 x_2 + \varepsilon y = β 1 x 1 + β 2 x 2 + ε α ^ 1 \hat{\alpha}_1 α ^ 1
α ^ 1 = Cov ( x 1 , β 1 x 1 + β 2 x 2 + ε ) Var ( x 1 ) = β 1 Cov ( x 1 , x 1 ) + β 2 Cov ( x 1 , x 2 ) + Cov ( x 1 , ε ) Var ( x 1 ) = β 1 Var ( x 1 ) Var ( x 1 ) + β 2 Cov ( x 1 , x 2 ) Var ( x 1 ) + Cov ( x 1 , ε ) Var ( x 1 ) = β 1 + β 2 Cov ( x 1 , x 2 ) Var ( x 1 ) + Cov ( x 1 , ε ) Var ( x 1 )
\begin{aligned}
\hat{\alpha}_1 &= \frac{\text{Cov}(x_1, \beta_1 x_1 + \beta_2 x_2 + \varepsilon)}{\text{Var}(x_1)} \\
&= \frac{\beta_1 \text{Cov}(x_1, x_1) + \beta_2 \text{Cov}(x_1, x_2) + \text{Cov}(x_1, \varepsilon)}{\text{Var}(x_1)} \\
&= \beta_1 \frac{\text{Var}(x_1)}{\text{Var}(x_1)} + \beta_2 \frac{\text{Cov}(x_1, x_2)}{\text{Var}(x_1)} + \frac{\text{Cov}(x_1, \varepsilon)}{\text{Var}(x_1)} \\
&= \beta_1 + \beta_2 \frac{\text{Cov}(x_1, x_2)}{\text{Var}(x_1)} + \frac{\text{Cov}(x_1, \varepsilon)}{\text{Var}(x_1)}
\end{aligned}
α ^ 1 = Var ( x 1 ) Cov ( x 1 , β 1 x 1 + β 2 x 2 + ε ) = Var ( x 1 ) β 1 Cov ( x 1 , x 1 ) + β 2 Cov ( x 1 , x 2 ) + Cov ( x 1 , ε ) = β 1 Var ( x 1 ) Var ( x 1 ) + β 2 Var ( x 1 ) Cov ( x 1 , x 2 ) + Var ( x 1 ) Cov ( x 1 , ε ) = β 1 + β 2 Var ( x 1 ) Cov ( x 1 , x 2 ) + Var ( x 1 ) Cov ( x 1 , ε ) 最後に、期待値はをとると以下のようになります（ここで「誤差項 ε \varepsilon ε x 1 x_1 x 1 E [ Cov ( x 1 , ε ) ] = 0 \mathbb{E}[\text{Cov}(x_1, \varepsilon)] = 0 E [ Cov ( x 1 , ε )] = 0
E [ α ^ 1 ] = β 1 + β 2 ⋅ Cov ( x 1 , x 2 ) Var ( x 1 ) ⏟ x 2 を x 1 に回帰した際の係数
\mathbb{E}[\hat{\alpha}_1] = \beta_1 + \beta_2 \cdot \underbrace{\frac{\text{Cov}(x_1, x_2)}{\text{Var}(x_1)}}_{x_2 \text{を} x_1 \text{に回帰した際の係数}}
E [ α ^ 1 ] = β 1 + β 2 ⋅ x 2 を x 1 に回帰した際の係数 Var ( x 1 ) Cov ( x 1 , x 2 ) 右辺の第2項が欠落変数バイアス です。
β 2 ≠ 0 \beta_2 \neq 0 β 2 = 0 x 2 x_2 x 2 y y y
Cov ( x 1 , x 2 ) ≠ 0 \text{Cov}(x_1, x_2) \neq 0 Cov ( x 1 , x 2 ) = 0 x 2 x_2 x 2 x 1 x_1 x 1
この2条件が揃ったとき、α ^ 1 \hat{\alpha}_1 α ^ 1 β 1 \beta_1 β 1
Frisch-Waugh-Lovell (FWL) 定理
次に、回帰分析において変数を追加する（コントロールする）と、因果推論に近づくのかを見ていきます。FWL定理 です。
FWL定理 y = X 1 β 1 + X 2 β 2 + ε y = \mathbf{X}_1 \boldsymbol{\beta}_1 + \mathbf{X}_2 \boldsymbol{\beta}_2 + \varepsilon y = X 1 β 1 + X 2 β 2 + ε β ^ 1 \hat{\boldsymbol{\beta}}_1 β ^ 1
y y y X 2 \mathbf{X}_2 X 2 y ~ \tilde{y} y ~
X 1 \mathbf{X}_1 X 1 X 2 \mathbf{X}_2 X 2 X ~ 1 \tilde{\mathbf{X}}_1 X ~ 1
y ~ \tilde{y} y ~ X ~ 1 \tilde{\mathbf{X}}_1 X ~ 1
この定理の主張は、「コントロールする」とは、悪影響を及ぼす変数の影響を事前に削除する ことです。
このプロセスを、情報（分散）の重なりを示すベン図で視覚的に理解してみましょう。
Step 1 & 2の意味（フィルター）y y y x 1 x_1 x 1 x 2 x_2 x 2 x 2 x_2 x 2 x 2 x_2 x 2 x 2 x_2 x 2 x 2 x_2 x 2 y ~ \tilde{y} y ~ x ~ 1 \tilde{x}_1 x ~ 1
Step 3の意味（純粋な相関）x 1 x_1 x 1 x 2 x_2 x 2
FWL定理を知っていると、重回帰分析をすれば因果関係がわかるという安易な考えを持たなくなるでしょう。私たちが本当に知りたい因果を抽出するためには、「モデルに入れていない未観測の変数」の影響までも、残差から完璧に削ぎ落とさなければならないのですが、データが存在しない以上それは不可能です。
だからこそ因果推論では、単なる回帰分析を諦め、「自然実験」や「操作変数」といった特殊なアプローチ（次節の識別戦略）を用いて、この「純粋な部分」を頑張って抽出しようと試みるのです。
識別戦略（Identification Strategy）
ここで紹介する戦略はすべて、擬似的なランダム化 を作り出すための工夫です。
1. 操作変数法（IV / 2SLS）
X X X ε \varepsilon ε Z Z Z Z Z Z X X X
2段階最小二乗法 (2SLS) : 第1段階で X ^ = Proj ( X ∣ Z ) \hat{X} = \text{Proj}(X \mid Z) X ^ = Proj ( X ∣ Z ) y y y X ^ \hat{X} X ^
2. 固定効果モデル（Fixed Effects）
パネルデータを用いて、個体固有の「変わらない特性（性格、企業の社風など）」を、個体ごとの平均からの差分をとることで取り除きます。
( y i t − y ˉ i ) = β ( x i t − x ˉ i ) + ( ε i t − ε ˉ i )
(y_{it} - \bar{y}_i) = \beta (x_{it} - \bar{x}_i) + (\varepsilon_{it} - \bar{\varepsilon}_i)
( y i t − y ˉ i ) = β ( x i t − x ˉ i ) + ( ε i t − ε ˉ i ) これにより、観測不可能な時不変の交絡因子によるバイアスを取り除きます。
3. 差分の差分法（DID）
「政策が行われたグループ」と「行われなかったグループ」の、前後の変化の「差」を比較します。
τ ^ D I D = ( 介入群の後 − 介入群の前 ) − ( 対照群の後 − 対照群の前 )
\hat{\tau}_{DID} = (\text{介入群の後} - \text{介入群の前}) - (\text{対照群の後} - \text{対照群の前})
τ ^ D I D = ( 介入群の後 − 介入群の前 ) − ( 対照群の後 − 対照群の前 ) 共通のトレンドさえ仮定（平行トレンド仮定）できれば、説得力のある方法です。
4. 回帰不連続デザイン（RDD）
「テストが60点以上なら合格」といった閾値（カットオフ）を利用します。59点の人と61点の人は「ほぼ同じ」であるという仮定に基づき、カットオフ付近での y y y
Level 4のまとめ
このレベルに達した分析者は、出力された回帰係数を安易に「効果」と解釈しません。
この係数は内生性によって過大評価されていないか？
どの変数をコントロールすべきで、どの変数を入れてはいけないのか（合流点バイアスなど）？
そもそもこの効果は識別可能なのか？
こうした因果的な問いを持ち、バイアスを適切に取り除けられれば（取り除けないことを理解できていれば）、Level 4です。
Level 5: 超高次元
Level 3 では、Lasso や Ridge を「予測精度を上げるための正則化」として見ました。p p p
p > n p > n p > n
線形回帰モデル
y = X β + ε , X ∈ R n × p
\mathbf y = \mathbf X \boldsymbol{\beta} + \boldsymbol{\varepsilon},
\qquad
\mathbf X \in \mathbb{R}^{n \times p}
y = X β + ε , X ∈ R n × p を考えます。ここで p > n p > n p > n X \mathbf X X n n n
r a n k ( X ) ≤ n < p
\mathrm{rank}(\mathbf X) \le n < p
rank ( X ) ≤ n < p です。したがって X ⊤ X \mathbf X^\top \mathbf{X} X ⊤ X
β ^ O L S = ( X ⊤ X ) − 1 X ⊤ y
\hat{\boldsymbol{\beta}}_{\mathrm{OLS}}
=
(\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top \mathbf y
β ^ OLS = ( X ⊤ X ) − 1 X ⊤ y は使えません。
使えないというのは、最小二乗問題の解が存在しないのではなく、解が一意に定まらない のです。 実際、β ^ 0 \hat{\boldsymbol{\beta}}_0 β ^ 0 v ∈ N u l l ( X ) \mathbf v \in \mathrm{Null}(\mathbf X) v ∈ Null ( X )
X ( β ^ 0 + v ) = X β ^ 0
\mathbf X(\hat{\boldsymbol{\beta}}_0 + \mathbf v)
=
\mathbf X\hat{\boldsymbol{\beta}}_0
X ( β ^ 0 + v ) = X β ^ 0 なので、β ^ 0 + v \hat{\boldsymbol{\beta}}_0 + \mathbf v β ^ 0 + v 同じ y ^ \hat{\mathbf y} y ^ のです。
このとき必要になるのが、「その中でどの解をもっともらしいと見なすか」という追加の最適化問題 です。そこで重要になってくるのがスパース性（sparsity） です。
以後、必ずしも p > n p > n p > n p p p
スパース性と Lasso
高次元では、説明変数の総数 p p p s s s
s = ∣ s u p p ( β ∗ ) ∣ ≪ n
s = |\mathrm{supp}(\boldsymbol{\beta}^\ast)| \ll n
s = ∣ supp ( β ∗ ) ∣ ≪ n であると仮定することが多いです。これを スパース性 と呼びます。β ∗ \boldsymbol{\beta}^\ast β ∗ s u p p ( β ∗ ) \mathrm{supp}(\boldsymbol{\beta}^\ast) supp ( β ∗ ) β ∗ \boldsymbol{\beta}^\ast β ∗
この仮定のもとで代表的なのが Lasso です。
β ^ L a s s o = arg min β { 1 2 n ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 }
\hat{\boldsymbol{\beta}}_{\mathrm{Lasso}}
=
\arg\min_{\boldsymbol{\beta}}
\left\{
\frac{1}{2n}\|\mathbf y - \mathbf X\boldsymbol{\beta}\|_2^2
+
\lambda \|\boldsymbol{\beta}\|_1
\right\}
β ^ Lasso = arg β min { 2 n 1 ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 } L 1 L^1 L 1 ( ∥ β ∥ 1 = ∑ j ∣ β j ∣ (\|\boldsymbol{\beta}\|_1 = \sum_j |\beta_j| ( ∥ β ∥ 1 = ∑ j ∣ β j ∣
Level 3 ではこれを「予測のための正則化」と見ましたが、Level 5 ではさらに一歩進めて、どの条件のもとで、Lasso が真の β ∗ \boldsymbol{\beta}^\ast β ∗
1. 変数推定と変数選択は別物
Lasso について議論するとき、最初に区別すべきは次の 2 つです。
(A) 推定された係数 β ^ \hat\beta β ^ β ∗ \beta^\ast β ∗
(B) 推定されたサポート S ^ = supp ( β ^ ) \hat S = \operatorname{supp}(\hat\beta) S ^ = supp ( β ^ ) S = supp ( β ∗ ) S = \operatorname{supp}(\beta^\ast) S = supp ( β ∗ )
(A) は予測・推定誤差の問題、(B) は変数選択の一貫性 (selection consistency) の問題で、要求される条件のレベルが全く違います。
1-1. 推定誤差が小さくなるための条件: Restricted Eigenvalue 条件
Lasso の解 β ^ L a s s o \hat\beta_{\mathrm{Lasso}} β ^ Lasso
β ^ L a s s o = arg min β ∈ R p { 1 2 n ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 }
\hat\beta_{\mathrm{Lasso}}
= \arg\min_{\beta \in \mathbb R^p}
\left\{ \frac{1}{2n}\|\mathbf y - \mathbf X \beta\|_2^2 + \lambda \|\beta\|_1 \right\}
β ^ Lasso = arg β ∈ R p min { 2 n 1 ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 } を解いて得られます。p > n p > n p > n Σ ^ = X ⊤ X / n \hat\Sigma = \mathbf X^\top \mathbf X / n Σ ^ = X ⊤ X / n
Lasso の最適性条件 (KKT 条件)
Lasso の解 β ^ L a s s o \hat\beta_{\mathrm{Lasso}} β ^ Lasso KKT 条件 (Karush–Kuhn–Tucker 条件) を考えます。
【KKT 条件とは】
凸最適化問題 min x f ( x ) \min_x f(x) min x f ( x ) x ^ \hat x x ^ ∇ f ( x ^ ) = 0 \nabla f(\hat x) = 0 ∇ f ( x ^ ) = 0 最適性の必要十分条件 です。
Lasso の目的関数
F ( β ) = 1 2 n ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1
F(\beta) = \frac{1}{2n}\|\mathbf y - \mathbf X \beta\|_2^2 + \lambda \|\beta\|_1
F ( β ) = 2 n 1 ∥ y − X β ∥ 2 2 + λ ∥ β ∥ 1 は、ℓ 1 \ell_1 ℓ 1 ∥ β ∥ 1 = ∑ j ∣ β j ∣ \|\beta\|_1 = \sum_j |\beta_j| ∥ β ∥ 1 = ∑ j ∣ β j ∣ β j = 0 \beta_j = 0 β j = 0 劣微分 (subdifferential) を使った KKT 条件を考えます。j j j β ^ = β ^ L a s s o \hat\beta = \hat\beta_{\mathrm{Lasso}} β ^ = β ^ Lasso j = 1 , … , p j = 1, \ldots, p j = 1 , … , p
1 n x j ⊤ ( y − X β ^ ) ⏟ 特徴量 j と残差の内積 / n = λ s j , s j ∈ { { sign ( β ^ j ) } ( β ^ j ≠ 0 のとき ) [ − 1 , 1 ] ( β ^ j = 0 のとき )
\underbrace{\frac{1}{n}\,\mathbf x_j^\top (\mathbf y - \mathbf X \hat\beta)}_{\text{特徴量 } j \text{ と残差の内積}/n}
=
\lambda\, s_j,
\qquad
s_j \in
\begin{cases}
\{\operatorname{sign}(\hat\beta_j)\} & (\hat\beta_j \ne 0 \text{ のとき}) \\
[-1, 1] & (\hat\beta_j = 0 \text{ のとき})
\end{cases}
特徴量 j と残差の内積 / n n 1 x j ⊤ ( y − X β ^ ) = λ s j , s j ∈ { { sign ( β ^ j )} [ − 1 , 1 ] ( β ^ j = 0 のとき ) ( β ^ j = 0 のとき ) が成り立ちます。ここで x j ∈ R n \mathbf x_j \in \mathbb R^n x j ∈ R n X \mathbf X X j j j j j j y − X β ^ \mathbf y - \mathbf X \hat\beta y − X β ^
この式の意味を考えます。∣ s j ∣ ≤ 1 |s_j| \le 1 ∣ s j ∣ ≤ 1 j j j
∣ 1 n x j ⊤ ( y − X β ^ ) ∣ ≤ λ
\left|\frac{1}{n}\,\mathbf x_j^\top (\mathbf y - \mathbf X \hat\beta)\right| \le \lambda
n 1 x j ⊤ ( y − X β ^ ) ≤ λ が成立します。1 n x j ⊤ ( y − X β ^ ) \frac{1}{n}\mathbf x_j^\top (\mathbf y - \mathbf X \hat\beta) n 1 x j ⊤ ( y − X β ^ ) x j \mathbf x_j x j 残差と特徴量の標本相関に対応する量 です。つまり KKT 条件は「Lasso 解 β ^ \hat\beta β ^ どの特徴量とも相関の大きさが λ \lambda λ になっている」ということを述べています。これは「Lasso が残差をできるだけ各特徴量と無相関にするように β ^ \hat\beta β ^
コーン制約の導出
次に、誤差ベクトル δ = β ^ − β ∗ \boldsymbol\delta = \hat\beta - \beta^\ast δ = β ^ − β ∗
β ^ \hat\beta β ^ F ( β ^ ) ≤ F ( β ∗ ) F(\hat\beta) \le F(\beta^\ast) F ( β ^ ) ≤ F ( β ∗ )
1 2 n ∥ y − X β ^ ∥ 2 2 + λ ∥ β ^ ∥ 1 ≤ 1 2 n ∥ y − X β ∗ ∥ 2 2 + λ ∥ β ∗ ∥ 1
\frac{1}{2n}\|\mathbf y - \mathbf X \hat\beta\|_2^2 + \lambda \|\hat\beta\|_1
\le
\frac{1}{2n}\|\mathbf y - \mathbf X \beta^\ast\|_2^2 + \lambda \|\beta^\ast\|_1
2 n 1 ∥ y − X β ^ ∥ 2 2 + λ ∥ β ^ ∥ 1 ≤ 2 n 1 ∥ y − X β ∗ ∥ 2 2 + λ ∥ β ∗ ∥ 1 が成り立ちます。y = X β ∗ + ε \mathbf y = \mathbf X \beta^\ast + \boldsymbol\varepsilon y = X β ∗ + ε
1 2 n ∥ X δ ∥ 2 2 ≤ 1 n ε ⊤ X δ + λ ( ∥ β ∗ ∥ 1 − ∥ β ^ ∥ 1 )
\frac{1}{2n}\|\mathbf X \boldsymbol\delta\|_2^2
\le
\frac{1}{n}\,\boldsymbol\varepsilon^\top \mathbf X \boldsymbol\delta
+ \lambda \bigl( \|\beta^\ast\|_1 - \|\hat\beta\|_1 \bigr)
2 n 1 ∥ X δ ∥ 2 2 ≤ n 1 ε ⊤ X δ + λ ( ∥ β ∗ ∥ 1 − ∥ β ^ ∥ 1 ) が得られます。チューニングパラメータ λ \lambda λ 1 n ∥ X ⊤ ε ∥ ∞ \frac{1}{n}\|\mathbf X^\top \boldsymbol\varepsilon\|_\infty n 1 ∥ X ⊤ ε ∥ ∞ λ ≥ 2 n ∥ X ⊤ ε ∥ ∞ \lambda \ge \frac{2}{n}\|\mathbf X^\top \boldsymbol\varepsilon\|_\infty λ ≥ n 2 ∥ X ⊤ ε ∥ ∞ ℓ 1 \ell_1 ℓ 1 ∥ β ∥ 1 = ∥ β S ∥ 1 + ∥ β S c ∥ 1 \|\beta\|_1 = \|\beta_S\|_1 + \|\beta_{S^c}\|_1 ∥ β ∥ 1 = ∥ β S ∥ 1 + ∥ β S c ∥ 1
1 2 n ∥ X δ ∥ 2 2 ≤ 3 λ 2 ∥ δ S ∥ 1 − λ 2 ∥ δ S c ∥ 1
\frac{1}{2n} \|\mathbf X \boldsymbol\delta\|_2^2
\le
\frac{3\lambda}{2} \|\boldsymbol\delta_S\|_1 - \frac{\lambda}{2} \|\boldsymbol\delta_{S^c}\|_1
2 n 1 ∥ X δ ∥ 2 2 ≤ 2 3 λ ∥ δ S ∥ 1 − 2 λ ∥ δ S c ∥ 1 という basic inequality が導けます。
左辺は ∥ X δ ∥ 2 2 ≥ 0 \|\mathbf X \boldsymbol\delta\|_2^2 \ge 0 ∥ X δ ∥ 2 2 ≥ 0
λ 2 ∥ δ S c ∥ 1 ≤ 3 λ 2 ∥ δ S ∥ 1 ⟺ ∥ δ S c ∥ 1 ≤ 3 ∥ δ S ∥ 1
\frac{\lambda}{2} \|\boldsymbol\delta_{S^c}\|_1 \le \frac{3\lambda}{2} \|\boldsymbol\delta_S\|_1
\quad \Longleftrightarrow \quad
\|\boldsymbol\delta_{S^c}\|_1 \le 3 \|\boldsymbol\delta_S\|_1
2 λ ∥ δ S c ∥ 1 ≤ 2 3 λ ∥ δ S ∥ 1 ⟺ ∥ δ S c ∥ 1 ≤ 3∥ δ S ∥ 1 を得ます。これが、本記事で扱う コーン制約 C ( S , 3 ) \mathcal C(S, 3) C ( S , 3 )
コーン制約の意味
この不等式は、誤差ベクトル δ = β ^ − β ∗ \boldsymbol\delta = \hat\beta - \beta^\ast δ = β ^ − β ∗
∥ δ S c ∥ 1 ⏟ 非サポートでの誤差の総量 ≤ 3 ⋅ ∥ δ S ∥ 1 ⏟ サポートでの誤差の総量
\underbrace{\|\boldsymbol\delta_{S^c}\|_1}_{\text{非サポートでの誤差の総量}}
\le 3 \cdot
\underbrace{\|\boldsymbol\delta_S\|_1}_{\text{サポートでの誤差の総量}}
非サポートでの誤差の総量 ∥ δ S c ∥ 1 ≤ 3 ⋅ サポートでの誤差の総量 ∥ δ S ∥ 1 を主張しており、「Lasso の誤差は、真サポート方向の大きさの 高々 3 倍までしか 、非サポート方向に広がらない」ことを意味します。これは ℓ 1 \ell_1 ℓ 1
この結果は Lasso の誤差解析にとって重要です。p > n p > n p > n X \mathbf X X R p \mathbb R^p R p ∥ X δ ∥ 2 2 ≳ ∥ δ ∥ 2 2 \|\mathbf X \boldsymbol\delta\|_2^2 \gtrsim \|\boldsymbol\delta\|_2^2 ∥ X δ ∥ 2 2 ≳ ∥ δ ∥ 2 2 誤差ベクトルがコーン C ( S , 3 ) \mathcal C(S, 3) C ( S , 3 ) 以上、コーンの内側で ∥ X δ ∥ 2 2 ≳ ∥ δ ∥ 2 2 \|\mathbf X \boldsymbol\delta\|_2^2 \gtrsim \|\boldsymbol\delta\|_2^2 ∥ X δ ∥ 2 2 ≳ ∥ δ ∥ 2 2
Restricted Eigenvalue (RE) 条件
真のサポートを S S S ∣ S ∣ = s |S| = s ∣ S ∣ = s κ > 0 \kappa > 0 κ > 0
C ( S , 3 ) = { δ ∈ R p : ∥ δ S c ∥ 1 ≤ 3 ∥ δ S ∥ 1 }
\mathcal C(S, 3) = \left\{ \boldsymbol\delta \in \mathbb R^p : \|\boldsymbol\delta_{S^c}\|_1 \le 3 \|\boldsymbol\delta_S\|_1 \right\}
C ( S , 3 ) = { δ ∈ R p : ∥ δ S c ∥ 1 ≤ 3∥ δ S ∥ 1 } 上のすべての δ ≠ 0 \boldsymbol\delta \ne 0 δ = 0
∥ X δ ∥ 2 n ≥ κ ∥ δ S ∥ 2
\frac{\|\mathbf X \boldsymbol\delta\|_2}{\sqrt n} \ge \kappa \|\boldsymbol\delta_S\|_2
n ∥ X δ ∥ 2 ≥ κ ∥ δ S ∥ 2 が成り立ちます。コーンの広さを決める定数として、ここでは 3 を採用しています(よく使われる値で、5 を使う流儀もあります)。
1-2. RE 条件のもとでの誤差レート
RE 条件と、適切なチューニングパラメータ
λ ≍ σ log p n
\lambda \asymp \sigma \sqrt{\frac{\log p}{n}}
λ ≍ σ n log p のもとで、Lasso は典型的に
∥ β ^ L a s s o − β ∗ ∥ 1 = O p ( s log p n ) , ∥ β ^ L a s s o − β ∗ ∥ 2 = O p ( s log p n )
\|\hat\beta_{\mathrm{Lasso}} - \beta^\ast\|_1 = O_p\!\left( s\sqrt{\frac{\log p}{n}} \right),
\qquad
\|\hat\beta_{\mathrm{Lasso}} - \beta^\ast\|_2 = O_p\!\left( \sqrt{\frac{s\log p}{n}} \right)
∥ β ^ Lasso − β ∗ ∥ 1 = O p ( s n log p ) , ∥ β ^ Lasso − β ∗ ∥ 2 = O p ( n s log p ) というレートで真の係数に近づくことが示されます。これはサンプル数の対数オーダーで次元 p p p ということです。
ここで強調すべきは、RE 条件は「Lasso の推定誤差を抑える条件」であり、変数選択の完全一致を保証する条件ではない ということです。
1-3. 真の変数集合を当てるための条件: Irrepresentable Condition
S ^ = s u p p ( β ^ L a s s o ) \hat S = \mathrm{supp}(\hat\beta_{\mathrm{Lasso}}) S ^ = supp ( β ^ Lasso ) S S S Irrepresentable Condition (IC) です。グラム行列を Σ = X ⊤ X / n \boldsymbol\Sigma = \mathbf X^\top \mathbf X / n Σ = X ⊤ X / n S S S S c S^c S c
∥ Σ S c S Σ S S − 1 sign ( β S ∗ ) ∥ ∞ ≤ 1 − η ( η > 0 )
\left\| \boldsymbol\Sigma_{S^c S}\, \boldsymbol\Sigma_{SS}^{-1}\, \operatorname{sign}(\beta^\ast_S) \right\|_\infty \le 1 - \eta
\qquad (\eta > 0)
Σ S c S Σ S S − 1 sign ( β S ∗ ) ∞ ≤ 1 − η ( η > 0 ) を要求します。直感的な理解としては、「本当は不要な変数が、本当に必要な変数を強く模倣できないこと」を要求している感じです。Σ S c S Σ S S − 1 \boldsymbol\Sigma_{S^c S}\boldsymbol\Sigma_{SS}^{-1} Σ S c S Σ S S − 1
さらに、信号強度も十分大きくなければいけません。具体的には、
min j ∈ S ∣ β j ∗ ∣ ≳ λ
\min_{j \in S} |\beta_j^\ast| \gtrsim \lambda
j ∈ S min ∣ β j ∗ ∣ ≳ λ のような beta-min 条件 が必要です。これは「真に非ゼロな係数は、ノイズに紛れない程度に大きい」という条件です。
1-4. 強さの順序
二つの条件には強弱の関係 があります。
IC ⟹ RE ( 一般に逆は成立しない )
\text{IC} \;\Longrightarrow\; \text{RE} \quad (\text{一般に逆は成立しない})
IC ⟹ RE ( 一般に逆は成立しない ) つまり、相関が中程度まで強くなると、予測と推定はうまくいくが、変数選択は失敗する というレジームに入ります。これが「Lasso は予測には強いが、選ばれた変数の解釈には注意」という実務の鉄則の理論的背景です。
1-5. 数値実験で確認する
理屈だけだとピンとこないので、実際に確認してみましょう。equicorrelation 構造（全特徴量が一様に相関 ρ \rho ρ ρ \rho ρ
複数の ρ \rho ρ
相関 ρ \rho ρ ρ \rho ρ ∥ β ^ − β ∗ ∥ 2 \|\hat\beta - \beta^\ast\|_2 ∥ β ^ − β ∗ ∥ 2 S ^ = S \hat S = S S ^ = S
この実験結果には 3 つのレジームがあり、それぞれ理論と対応しています。
(A) 低相関領域 (ρ ≲ 0.1 \rho \lesssim 0.1 ρ ≲ 0.1 λ \lambda λ
(B) 中相関領域 (0.1 ≲ ρ ≲ 0.6 0.1 \lesssim \rho \lesssim 0.6 0.1 ≲ ρ ≲ 0.6 最も興味深いレジーム です。サポート完全回復率は急落し、ρ ≈ 0.4 \rho \approx 0.4 ρ ≈ 0.4 ところが推定誤差はむしろ減少 しており、ρ ≈ 0.4 \rho \approx 0.4 ρ ≈ 0.4 選択は失敗しているが推定は成功している 」というレジームで、Lasso と多重共線性に関する最も実務的に重要な現象です。
(C) 強相関領域 (ρ ≳ 0.6 \rho \gtrsim 0.6 ρ ≳ 0.6 ρ → 1 \rho \to 1 ρ → 1 1 − ρ 1 - \rho 1 − ρ ρ → 1 \rho \to 1 ρ → 1 κ \kappa κ
1-6. 実務的教訓
この実験から学べることをまとめます。
「Lasso は変数選択をしてくれる」は誤解 : 変数間に相関があると IC が成り立たず、サポート回復は信頼できなくなります。
予測・推定だけが目的なら Lasso はそこそこ頑健 : 中程度の相関までは推定誤差は大きく崩れません。
極端な多重共線性は両方を壊す : ρ → 1 \rho \to 1 ρ → 1
この区別を理解していないと、「Lasso で選ばれた変数は重要な変数」などと誤った主張をする恐れがあります。
2. Post-Selection Inference の罠
Lasso に関して、実務でついやってしまいがちな罠があります。
Lasso で変数選択して、選ばれた変数だけで OLS を回し、普通の標準誤差で有意性検定する
これはそのままではNGです 。
2-1. なぜダメなのか
理由は単純で、同じデータを使って「選ぶ」ことと「検定する」ことを両方やっている からです。選択イベント自体が y \mathbf y y
直感的には、「たまたま大きく見えた係数だけを拾って、そのあとで“やっぱり大きいですね”と検定している」状態です。これはいわゆる winner's curse であり、p 値は過小評価され、信頼区間の被覆率もいい加減になります。
数式で書くと、通常の信頼区間 C I j naive CI_j^{\text{naive}} C I j naive
P ( β j ∗ ∈ C I j naive ∣ j ∈ S ^ ) ≠ 1 − α
P\!\left( \beta_j^\ast \in CI_j^{\text{naive}} \mid j \in \hat S \right) \neq 1 - \alpha
P ( β j ∗ ∈ C I j naive ∣ j ∈ S ^ ) = 1 − α です。つまり、「選ばれた後」の条件付きでは、通常の 95% 信頼区間は 95% を保たない のです。
2-2. どう対処するか
大きく次の 3 つの方法があります。
方法
アイデア
データ分割 (sample splitting) 片方のデータで変数選択、もう片方のデータで推論する
選択的推論 (selective inference) 「その変数が選ばれた」という事象に条件付けて推論する
デバイアスド Lasso Lasso のバイアスを解析的に補正して、漸近正規性を回復する
上2つはイメージしやすいと思うので、次節では3 つ目のデバイアスド Lasso を詳説します。
3. デバイアスド Lasso
Lasso は L 1 L^1 L 1 縮小による一次のバイアスを打ち消す補正 を加えるのがデバイアスド Lasso のアイデアです。
3-1. 補正
デバイアスド Lasso の推定量は
β ^ d b = β ^ L a s s o + Θ ^ X ⊤ ( y − X β ^ L a s s o ) n
\hat\beta_{\mathrm{db}}
= \hat\beta_{\mathrm{Lasso}}
+ \hat{\boldsymbol\Theta}\, \frac{\mathbf X^\top (\mathbf y - \mathbf X \hat\beta_{\mathrm{Lasso}})}{n}
β ^ db = β ^ Lasso + Θ ^ n X ⊤ ( y − X β ^ Lasso ) と書けます。ここで Θ ^ \hat{\boldsymbol\Theta} Θ ^ X ⊤ X / n \mathbf X^\top \mathbf X / n X ⊤ X / n nodewise Lasso で構成します。
理解を深めるために展開すると:
β ^ d b − β ∗ = ( β ^ L a s s o − β ∗ ) + Θ ^ X ⊤ X n ( β ∗ − β ^ L a s s o ) + Θ ^ X ⊤ ϵ n
\hat\beta_{\mathrm{db}} - \beta^\ast
= (\hat\beta_{\mathrm{Lasso}} - \beta^\ast)
+ \hat{\boldsymbol\Theta} \frac{\mathbf X^\top \mathbf X}{n}(\beta^\ast - \hat\beta_{\mathrm{Lasso}})
+ \hat{\boldsymbol\Theta} \frac{\mathbf X^\top \boldsymbol\epsilon}{n}
β ^ db − β ∗ = ( β ^ Lasso − β ∗ ) + Θ ^ n X ⊤ X ( β ∗ − β ^ Lasso ) + Θ ^ n X ⊤ ϵ 第 1 項と第 2 項を合わせると ( I − Θ ^ Σ ^ ) ( β ^ L a s s o − β ∗ ) \bigl(I - \hat{\boldsymbol\Theta} \hat\Sigma\bigr)(\hat\beta_{\mathrm{Lasso}} - \beta^\ast) ( I − Θ ^ Σ ^ ) ( β ^ Lasso − β ∗ ) Θ ^ ≈ Σ ^ − 1 \hat{\boldsymbol\Theta} \approx \hat\Sigma^{-1} Θ ^ ≈ Σ ^ − 1 Θ ^ X ⊤ ϵ / n \hat{\boldsymbol\Theta} \mathbf X^\top \boldsymbol\epsilon / n Θ ^ X ⊤ ϵ / n
3-2. 漸近正規性
適切な正則条件のもとで、各係数 j j j
n β ^ d b , j − β j ∗ τ ^ j ⟶ d N ( 0 , 1 )
\sqrt n \, \frac{\hat\beta_{\mathrm{db}, j} - \beta_j^\ast}{\hat\tau_j} \overset{d}{\longrightarrow} \mathcal N(0, 1)
n τ ^ j β ^ db , j − β j ∗ ⟶ d N ( 0 , 1 ) が成り立ちます。ここで τ ^ j \hat\tau_j τ ^ j 高次元でも個別係数について信頼区間や仮説検定ができる ようになります。
もちろん、これは無条件に成立する魔法ではありません。設計行列の条件、真の係数の疎性、さらには Θ \boldsymbol\Theta Θ
3-3. 実装と被覆率の確認
scikit-learn は直接サポートしていませんが、nodewise Lasso を使って自前で実装できます。
Naive (Lasso → OLS → 通常の CI) と Debiased Lasso で、それぞれ 95% 信頼区間の被覆率と平均幅を比較した結果が次の図です。
左: 95% 信頼区間の被覆率。点線が公称値 95%、これに近いほど良い結果です。Naive (赤) はどの ρ \rho ρ 、ρ \rho ρ ρ ≤ 0.6 \rho \le 0.6 ρ ≤ 0.6 ρ = 0.8 \rho = 0.8 ρ = 0.8
メッセージを 3 つに整理すると:
Naive は失敗 : 95% を一貫して下回り、相関とともに悪化していることが確認できます。これが正に Post-Selection Inference の罠です。
Debiased Lasso は理論に近い : 中相関までは95% に近い被覆を保ちます。
極端な相関には限界がある : 極端な高相関 (ρ = 0.8 \rho = 0.8 ρ = 0.8 Θ \Theta Θ
4. Double / Debiased Machine Learning (DML)
最後に、デバイアスドの考え方を機械学習による因果効果推定 にまで拡張した DML を紹介します。これは Chernozhukov et al. (2018) で体系化された枠組みです。
4-1. 部分線形モデル
シンプルな部分線形モデル を考えます。
Y = θ 0 D + g 0 ( X ) + U , D = m 0 ( X ) + V
Y = \theta_0 D + g_0(X) + U, \qquad D = m_0(X) + V
Y = θ 0 D + g 0 ( X ) + U , D = m 0 ( X ) + V
Y Y Y
D D D
X X X
θ 0 \theta_0 θ 0
g 0 , m 0 g_0, m_0 g 0 , m 0
U , V U, V U , V E [ U ∣ X , D ] = 0 E[U|X, D] = 0 E [ U ∣ X , D ] = 0 E [ V ∣ X ] = 0 E[V|X] = 0 E [ V ∣ X ] = 0
ここで g 0 ( X ) g_0(X) g 0 ( X ) D D D X X X Y Y Y
m 0 ( X ) = E [ D ∣ X ] , ℓ 0 ( X ) = E [ Y ∣ X ]
m_0(X) = E[D \mid X],
\qquad
\ell_0(X) = E[Y \mid X]
m 0 ( X ) = E [ D ∣ X ] , ℓ 0 ( X ) = E [ Y ∣ X ] このとき、
ℓ 0 ( X ) = θ 0 m 0 ( X ) + g 0 ( X )
\ell_0(X) = \theta_0 m_0(X) + g_0(X)
ℓ 0 ( X ) = θ 0 m 0 ( X ) + g 0 ( X ) であり、一般には g 0 ( X ) ≠ ℓ 0 ( X ) g_0(X) \neq \ell_0(X) g 0 ( X ) = ℓ 0 ( X )
θ 0 \theta_0 θ 0 g 0 , m 0 g_0, m_0 g 0 , m 0
したがって、
Y − ℓ 0 ( X ) = θ 0 { D − m 0 ( X ) } + U
Y - \ell_0(X)
=
\theta_0 \{D - m_0(X)\} + U
Y − ℓ 0 ( X ) = θ 0 { D − m 0 ( X )} + U となります。Y Y Y D D D X X X θ 0 \theta_0 θ 0
4-2. なぜ素朴に ML を入れるだけではダメなのか
g 0 ( X ) g_0(X) g 0 ( X ) m 0 ( X ) m_0(X) m 0 ( X )
θ ^ = ∑ i D i ( Y i − g ^ ( X i ) ) ∑ i D i 2
\hat\theta = \frac{\sum_i D_i (Y_i - \hat g(X_i))}{\sum_i D_i^2}
θ ^ = ∑ i D i 2 ∑ i D i ( Y i − g ^ ( X i )) のように θ 0 \theta_0 θ 0
g ^ \hat g g ^
同じデータで g ^ \hat g g ^
が θ ^ \hat\theta θ ^ 予測器としては優秀でも、推論対象 θ 0 \theta_0 θ 0 n \sqrt n n
4-3. DML の核心: 直交化と cross-fitting
DML はこの問題を 2 つのアイデアで解決します。
(1) Neyman 直交化
局外パラメータ η = ( g , m ) \eta = (g, m) η = ( g , m )
ψ ( W ; θ , η ) = ( D − m ( X ) ) { Y − g ( X ) − θ ( D − m ( X ) ) }
\psi(W; \theta, \eta) = (D - m(X))\, \bigl\{ Y - g(X) - \theta(D - m(X)) \bigr\}
ψ ( W ; θ , η ) = ( D − m ( X )) { Y − g ( X ) − θ ( D − m ( X )) } を使います。重要な性質は、η = ( g 0 , m 0 ) \eta = (g_0, m_0) η = ( g 0 , m 0 )
∂ ∂ η E [ ψ ( W ; θ 0 , η ) ] ∣ η = η 0 = 0
\left. \frac{\partial}{\partial \eta} E[\psi(W; \theta_0, \eta)] \right|_{\eta = \eta_0} = 0
∂ η ∂ E [ ψ ( W ; θ 0 , η )] η = η 0 = 0 となること(Neyman 直交性)です。つまり g , m g, m g , m θ \theta θ 一次では消える ように設計されています。
(2) Cross-fitting
データを K K K g ^ , m ^ \hat g, \hat m g ^ , m ^ 別 fold で学習 します。これにより、「同じデータで学習した複雑な予測器を、同じデータにそのまま当てる」ことによる過学習バイアスを取り除きます。
このとき、部分線形モデルの DML 推定量は実質的に
θ ^ D M L = 1 n ∑ i = 1 n ( Y i − g ^ ( X i ) ) ( D i − m ^ ( X i ) ) 1 n ∑ i = 1 n ( D i − m ^ ( X i ) ) 2
\hat\theta_{\mathrm{DML}}
= \frac{\frac{1}{n}\sum_{i=1}^n (Y_i - \hat g(X_i))(D_i - \hat m(X_i))}{\frac{1}{n}\sum_{i=1}^n (D_i - \hat m(X_i))^2}
θ ^ DML = n 1 ∑ i = 1 n ( D i − m ^ ( X i ) ) 2 n 1 ∑ i = 1 n ( Y i − g ^ ( X i )) ( D i − m ^ ( X i )) となります。これは見方を変えると、
X X X Y Y Y D D D
ことに相当します。この意味で DML は、Frisch–Waugh–Lovell (FWL) 定理の機械学習版 と見なせます。
4-4. 理論的保証
十分条件の一つとして、局外パラメータの推定誤差が
∥ g ^ − g 0 ∥ ⋅ ∥ m ^ − m 0 ∥ = o p ( n − 1 / 2 )
\|\hat g - g_0\|\cdot \|\hat m - m_0\| = o_p(n^{-1/2})
∥ g ^ − g 0 ∥ ⋅ ∥ m ^ − m 0 ∥ = o p ( n − 1/2 ) を満たせば、θ ^ D M L \hat\theta_{\mathrm{DML}} θ ^ DML n \sqrt n n n − 1 / 4 n^{-1/4} n − 1/4 o p ( 1 ) o_p(1) o p ( 1 ) n \sqrt n n
4-5. 実装で確認する
実装は驚くほどシンプルで、本質は数十行で書けます。
データ生成過程として g 0 , m 0 g_0, m_0 g 0 , m 0 θ 0 = 2.0 \theta_0 = 2.0 θ 0 = 2.0
Naive plug-in (赤) は θ 0 = 2.0 \theta_0 = 2.0 θ 0 = 2.0 θ ^ ≈ 0.43 \hat\theta \approx 0.43 θ ^ ≈ 0.43 − 1.57 -1.57 − 1.57 θ ^ ≈ 1.90 \hat\theta \approx 1.90 θ ^ ≈ 1.90 − 0.10 -0.10 − 0.10
差は劇的です。Naive plug-in は真値から大きく外れた値に集中しており、いくらサンプル数を増やしてもこのバイアスは消えません。一方 DML は真値の周りに集中しており、n \sqrt n n
4-6. DML で何が嬉しいのか
DML は、
柔軟な機械学習 で局外パラメータ g 0 , m 0 g_0, m_0 g 0 , m 0
興味のあるパラメータ θ 0 \theta_0 θ 0 厳密な統計推論 （漸近正規性、信頼区間）ができる
というところがとても嬉しいです。これは因果推論を機械学習時代にアップデートする上で重要な枠組みです。R/Python ともに DoubleML
https://docs.doubleml.org/
Level 5のまとめ
このレベルに達した分析者は、p p p
推定誤差を抑える条件（Restricted Eigenvalue）と、真の変数集合を当てる条件（Irrepresentable Condition + beta-min）の強弱関係を理解しているか？
「Lassoで変数を選んでからOLSで検定する」というナイーブな手続きがなぜダメなのか、Post-Selection Inferenceの罠を説明できるか？
デバイアスドLassoが、L 1 L^1 L 1
Double/Debiased Machine Learningが、Neyman直交化とcross-fittingという2つの仕掛けで、「機械学習の柔軟性」と「n \sqrt{n} n
「Lassoが選んだ変数だから重要だ」「選択後のp値だから信頼できる」といった主張に対して、それは推定誤差の話なのか変数選択の話なのか、相関構造のもとでICは成立しているのかと立ち止まれるようになれば、Level 5です。ここまで来ると、回帰分析はもはや単一の手法ではなく、予測・推論・因果という別々の目的に応じて使い分けるべき道具立ての集合体として見えてくるはずです。
Level 6: 哲学的問いへ
最後のレベルでは、これまでの「回帰」という概念そのものを問い直します。が、私がこのレベルには到底到達していないので書けません笑
こんなのがあるのかもしれませんし、ないのかもしれません。Level 5より上位の理解、または別観点の理解がある方は是非執筆していただけますと嬉しいです。
結び
本記事では、回帰分析の全体像を1つの地図として描くことを優先しました。そのため、各トピックの厳密な証明や実装上の細部は参考文献等に譲っています。もし途中で分からない箇所があれば、そこが次に深掘りすべきポイントです。
私はこの記事の執筆を通じて一層回帰分析への理解が深まりました。皆さんが回帰分析の奥深さを改めて知り、勉強のモチベーションが上がればこの記事としては大成功です。
参考文献
佐和 隆光, 回帰分析 , 朝倉書店, 2020.
安井翔太, 効果検証入門 , 技術評論社, 2020.
一般社団法人 日本統計学会, 日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック , 学術図書出版社, 2020.
久保川達也, 現代数理統計学の基礎 , 共立出版株式会社, 2017.
鈴木 讓, 統計的機械学習の数理100問 with R , 共立出版株式会社, 2020.
Peter J. Bickel, et. al., Simultaneous analysis of Lasso and Dantzig selector , 2008.
Zhao, P., & Yu, B., On model selection consistency of Lasso , 2006.
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. "Double/debiased machine learning for treatment and structural parameters. ", 2018.
Discussion
で，絶望の谷はどこにあるのでしょうか？