はじめに

回帰分析って知ってますか？

こんな質問されたらゾッとしますよね。知らないというわけでもないですが、きっとかなり深いレベルで知っていないとボコボコにされるでしょう。

R言語の lm() や Python の sklearn.LinearRegression() を呼び出して、推定した回帰係数を眺めたことがある、紙とペンでOLSによる係数を求めたことがある、OLS推定量がBLUEであることを確かめた、因果効果を回帰分析で推定した。。。どこまでやれば回帰分析を「知っている」と言って良いのでしょうか。

そんなことは何でも良いのです。とにかく回帰分析はシンプルながら、最も奥深いテーマの一つであり、多層的に理解する必要があります。

この記事では、回帰分析の理解をLevel 0（道具として使うだけ）からLevel 5（高次元データにおけるDMLなど）まで、6段階に分けて解説します。各レベルで「何が見えるようになるか」「どこに崖があるか」を丁寧に追っていきます。

思いの外分量のある記事になってしまったので、対象ごとにお勧めの章を整理しておきます。

回帰分析への理解は、ダニングクルーガー効果の絶望の谷を味わい続けるのにもってこいです。さて、本記事を通じて一度「絶望の谷」に落ちて、そこから一緒に「啓発の坂」を登っていきましょう。

https://hrmos.co/trend/talent-management/5397/ より引用。

Level 0: 道具としての回帰

最初の段階は、回帰をブラックボックスのツールとして使う段階です。

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)
print(model.coef_, model.score(X, y))

この程度の実装はAIで書けるレベルを通り越してもはやボイラープレート^[1]で、書けること自体に価値はありません。

回帰分析をこんな感じで捉えている人は多いでしょう：

回帰係数の符号と大きさを見て「学歴は効いてる」「経験年数も正に効いてる」と語ります。R²が0.7なら「全体的にまあまあ説明できてる」と判断する感じです^[2]。

この段階の限界は明らかで、以下の問いに答えられません。

• 係数が現実世界において何を意味するのか
• いつどんな時にその係数を信じていいのか
• 「効いてる」とはどういう意味か

ここから先のレベルで、これらの問いに段階的に答えられるようになっていきます。

Level 1: 古典的な統計的解釈

このレベルは簡単な統計学の教科書を一通り学び終えた水準です。回帰分析の理論は習得しているものの、深くは踏み込んでいない、統計検定でいうところの2級水準のイメージです。具体的には以下の事柄を理解している状態です。

最小二乗法の幾何学

線形回帰モデルを行列で書きます。

ここで \mathbf{y} \in \mathbb{R}^n、\mathbf{X} \in \mathbb{R}^{n \times p}（列フルランクを仮定）、\boldsymbol{\beta} \in \mathbb{R}^p。自然言語でいうと、 p 次元 のデータが n 件あって、スカラー値 y を回帰（予測）する問題設定です。

最小二乗推定量（OLS）は残差二乗和を最小化することで得られます（今後断りなしに推定したパラメータにはハット（^）をつけます）。

このOLS計算では予測誤差の二乗和を最小にする係数を求めていますが、幾何学的に理解することも重要です。

回帰分析の幾何学的理解。

このアニメーションを説明していきます。

空間設定

このアニメーションの空間は n 次元です。p ではないことに注意してください。

• 観測ベクトル y: n 個のサンプルの値を並べたベクトルです。これは n 次元空間 \mathbb{R}^n の中の「1点」を指します。
• データ行列 X: n \times p 行列であり、各列ベクトル x_1, x_2, \dots, x_p はすべて \mathbb{R}^n のベクトルです。

まずは黒矢印に着目してください。これは基底ベクトルです（x_1, x_2）。p 個存在します（今p=2）。基底ベクトルは行列 \mathbf{X} の列ベクトルで、各ベクトルはもちろん n 次元です。これらのベクトルが張る空間^[3]が、水色の平面です。別の言い方をすると列空間 M(X) です。線形回帰モデルは、「予測値はこの平面上のどこかに存在しなければならない」という制約を持っています。

念の為の補足ですが、これはあるデータ y_i に関するアニメーションではなく、n個のデータ全て対象にした回帰分析全体を表したアニメーションです。

データ生成プロセス

次にデータ生成プロセスを考えます。

• 真の平均 \eta（緑色の矢印）:
  自然界や社会における「真のルール」^[4]による値です。通常分析者は知る術がありません。

• モデル内成分 \eta_0（青色の矢印）:
  真の平均のうち、我々が用意した説明変数 x_1, x_2 で説明できる部分です。平面上に位置しています。

• バイアス \eta_1（オレンジ色の矢印）:
  真の平均のうち、用意した変数では説明しきれない「モデルの不備」による成分です。平面に対して垂直（直交）な成分として描かれています。本来は y の説明に必要だった説明変数 x_3, x_4 などを分析者が入手できなかったという状況です。この成分の分解を

  \eta = \eta_0 + \eta_1
  と記述します。いわゆるバイアスと呼ばれるものです。

• 観測値 y（赤色の矢印）:
  真の平均 \eta に、さらにノイズ（誤差 \epsilon）が加わった、私たちが実際に手にするデータです。データは観測誤差なども含めて、あらゆる原因で乗っかります。

最小二乗法による推定（直交射影）

私たちがやりたいことは、手元にある赤色の矢印 y から、尤もらしい平面上の点を特定することです。

• 推定値 \hat{y}（紫色の矢印）:
  観測値 y を平面 M(X) に向かって「垂直に」落とした点です。これを直交射影と呼びます。赤矢印の観測値 y のノイズによって、この紫の矢印は振り回されることになり、それをバリアンスと呼びます。
• なぜ垂直なのか？:
  平面上の点の中で y との距離（残差の二乗和）を最小にする点は、幾何学的に「垂直に下ろした足」になるからです。

幾何学的視点から見たバイアス・バリアンストレードオフ

説明変数が少なすぎると、モデルの表現力不足に陥ります。真の平均値と、モデルで表現できる平面との間の距離が開いてバイアスが大きくなります。一方、変数を増やしすぎてモデルが複雑になると、ノイズに反応しやすくなってしまいます。その分、平面自体が真の値に近付いてくれるはずなのでバイアスは小さくなります。

回帰分析の道具

この段階の理解をしている人は以下の事柄も理解していることでしょう。

検定と信頼区間

https://zenn.dev/zenkigen_tech/articles/8ffcd1b39e1e5f

最小二乗法そのものは、誤差項の分布に特定の仮定を置かなくても計算可能です。しかし、推定された係数が「偶然得られたものではないか」を統計的に検定するためには、誤差項に関する強い仮定が必要になります。

ここで、誤差項ベクトル \boldsymbol{\varepsilon} が互いに独立で、平均 \mathbf{0}、分散 \sigma^2 の正規分布に従うと仮定します。

この「正規性の仮定」を置くことで、推定量 \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} は線形変換の性質から、それ自身も多変量正規分布に従うことが導かれます。

この性質を起点に、回帰分析における代表的な2つの検定が構成されます。

1. 個別係数のt検定

「ある特定の説明変数が、目的変数に対して本当に影響を与えているか」を確認するための検定です。

通常、帰無仮説 H_0: \beta_j = 0（変数 x_j は影響を持たない）を立てます。
真の誤差分散 \sigma^2 は未知であるため、残差から求めた不偏分散 \hat{\sigma}^2 を用いて推定量の標準誤差 \mathrm{SE}(\hat{\beta}_j) を計算します。

この t 統計量は、自由度 n - k - 1 （n はサンプルサイズ、k は説明変数の数）の t分布に従います。算出されたp値が有意水準（例えば 0.05）を下回れば帰無仮説を棄却し、「その変数は意味がある可能性が高い」と考察できます。

2. モデル全体の有意性を評価する F検定

「構築した回帰モデル全体が、単なる定数モデル（平均値予測）よりも意味のある予測をしているか」を確認するための検定です。

帰無仮説として、定数項以外のすべての回帰係数がゼロである（H_0: \beta_1 = \beta_2 = \dots = \beta_k = 0）と仮定します。これは、決定係数 R^2 を用いて以下のような F 統計量として表現できます。

この統計量は、自由度 (k, n - k - 1) の F分布に従います。F値が十分に大きければ帰無仮説を棄却し、「少なくとも1つの説明変数は目的変数の説明に寄与している（モデル全体として意味がある）」のではないかと考えられます。複数係数の線形制約の同時検定（例：\beta_1 = \beta_2 かどうか）への拡張も、このF検定の枠組みで行われます。

自由度調整済み決定係数

決定係数 R^2：

ただし R^2 には、説明変数を増やせば必ず R^2 は上がる（か変わらない）という落とし穴があります。これを補正したのが自由度調整済み決定係数です：

n - k - 1 は残差の自由度です。説明変数の数 k が増えるほど、分母が小さくなります。そして分母が小さくなると、右側の項が大きくなり、結果として全体の {R}_{\text{adj}}^2 が引き下げられます。つまり、「新しい変数を追加して得られる精度の向上」が「変数を増やしたことによるペナルティ」を上回らない限り、この値は上がりません。

多重共線性

説明変数同士が強く相関していると、\mathbf{X}^\top\mathbf{X} の行列式が0に近づき、(\mathbf{X}^\top\mathbf{X})^{-1} の要素が極端に大きくなります。これが多重共線性です（いわゆるマルチコ^[5]）。

有名な診断指標は VIF（Variance Inflation Factor） です。

ここで R_j^2 は j 番目の説明変数を他の説明変数で回帰したときの決定係数です。VIF > 10 が要注意の目安です。

残差診断

最小二乗法（OLS）による推定が適切であるためには、誤差項に関して「線形性・等分散性・正規性・独立性」といった前提条件が満たされている必要があります^[7]。これらを視覚的に確認する作業が回帰診断です。

主に以下の4つのグラフを確認します。

• Residuals vs Fitted（残差 vs 予測値）：線形性と等分散性の確認。平均 0 の周りにランダムに分散していれば OK。データに順番やIDがあって、誤差に周期性があってもよろしくないです。
• Normal Q-Q：誤差項の正規性の確認。点が直線上に並んでいれば OK。
• Scale-Location：等分散性の確認。標準化残差の絶対値の平方根をプロットし、ばらつきに偏りがないかを見ます。
• Residuals vs Leverage：外れ値と影響力の確認。クックの距離（Cook's distance）が大きく、モデルに過度な影響を与えているデータ点がないかをチェックします。

回帰診断図。

Level 1 のまとめ

このレベルに達した分析者は、回帰分析を「ブラックボックスの関数呼び出し」ではなく、最小二乗法という最適化問題の解として、また n 次元空間における直交射影として捉えられえています。

• OLS推定量を行列で導出し、その幾何学的意味（列空間への射影）を説明できるか？
• t検定・F検定が「何を仮定して」「何を主張している」のかを区別できるか？
• R^2 と自由度調整済み R^2 の違い、多重共線性のVIF、残差診断の4プロットの意味を語れるか？

こうした道具立てを使いこなし、実務で統計的な分析が一通り回せるようになれば、Level 1 です。多くの統計学の教科書ではこの辺りまで書かれており、ある程度簡単な実務で統計的な処理をする分には問題がない水準と思います。

ただし、ここまでの理解は、データに対する仮定が満たされている場合に過ぎず、その仮定自体を問う段階にはまだ訪れておらず、先は長いです。

Level 2: ガウス=マルコフ定理と仮定の体系

ガウス=マルコフ定理の理論的背景と、その仮定が満たされなかった場合の対処法（処世術）を身につけていると Level 2 です。

ガウス=マルコフ定理

OLSが「良い推定量」である根拠は何か？その答えがガウス=マルコフ定理です。

BLUEは平たくいうと、不偏（平均的に間違わない）であるシンプルな計算式（線形）の中で、一番ブレが少なく信用できる推定量ということです。手元のデータで線形回帰する場合、いつものOLSを使うのが最も精度が良いということを言っています。

仮定（A1〜A4）：

1. 線形性：\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}
   目的変数 \mathbf{y} が、説明変数 \mathbf{X} とパラメータ \boldsymbol{\beta} の線形結合で表されるという仮定です。ここで重要なのは「パラメータ \boldsymbol{\beta} に対して線形である」ということであり、説明変数自体が x^2 や \log(x) のように非線形であっても構いません。また、これは「真のデータ生成プロセスがこの数式通りである（必要な変数が過不足なく含まれている）」という、モデルの正しい定式化自体を要求する強い仮定でもあります。
2. 外生性（厳密外生性）：\mathbb{E}[\boldsymbol{\varepsilon} \mid \mathbf{X}] = \mathbf{0}
   説明変数 \mathbf{X} がどんな値をとったとしても、誤差項 \boldsymbol{\varepsilon} の期待値はゼロである、という仮定です。直感的には「観測できないノイズ \boldsymbol{\varepsilon} の中に、説明変数 \mathbf{X} と連動して動くような隠れた要因（情報）が残っていない」ことを意味します。
   これが満たされない代表例が欠落変数バイアスです。例えば、賃金を予測するモデルで「個人の本来の能力」をデータ化できず誤差項に押し付けた場合、能力と学歴（説明変数）には相関があるため、この外生性の仮定が崩壊してしまいます。
3. 等分散性・無相関：\text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \sigma^2 \mathbf{I}
   行列 \sigma^2 \mathbf{I} （\mathbf{I} は単位行列）という美しい形は、2つのことを同時に仮定しています。
   • 等分散性：対角成分がすべて \sigma^2。つまり、説明変数の値によらず、誤差のばらつきの大きさが常に一定であること。例えば「企業の規模が大きくなるほど、利益の予測誤差も大きくなる」ようなデータではこの仮定に反します。
   • 無相関：非対角成分がすべて 0。あるサンプルの誤差が、別のサンプルの誤差に影響を与えないこと（自己相関がないこと）。時系列データ等ではしばしばこの仮定が破られます。
4. \mathbf{X} の列フルランク
   説明変数の中に「他の変数の完全な組み合わせで作れる変数」が存在しない、という数学的な仮定です（完全な多重共線性がないこと）。
   もし、身長をcmで表した変数とmで表した変数を同時にモデルに入れたり、カテゴリカル変数のダミー変数を「すべて」入れてしまったり（ダミー変数の罠）すると、各列が線形従属になります。その結果、行列 \mathbf{X}^\top\mathbf{X} の逆行列が存在しなくなり、計算式の上でOLS推定量 \hat{\boldsymbol{\beta}} を求めること自体が不可能になります。

BLUEであることの証明（導出）

なぜOLS推定量 \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} が最小分散を持つことを、任意の線形推定量と比較することで証明します。

任意の線形推定量を \tilde{\boldsymbol{\beta}} = \mathbf{C}\mathbf{y} と置きます（\mathbf{C} は確率変数ではない定数行列とします）。
ここで、\mathbf{C} を \mathbf{C} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top + \mathbf{D} と定義します（\mathbf{D} はOLSの係数行列からのズレを表す行列）。

1. 不偏性の条件
\tilde{\boldsymbol{\beta}} が不偏推定量であるためには、\mathbb{E}[\tilde{\boldsymbol{\beta}} \mid \mathbf{X}] = \boldsymbol{\beta} となる必要があります。

これが任意の \boldsymbol{\beta} で成り立つためには、\mathbf{C}\mathbf{X} = \mathbf{I} でなければなりません。
\mathbf{C} の定義を代入すると、

したがって、\tilde{\boldsymbol{\beta}} が不偏であるための必要十分条件は \mathbf{D}\mathbf{X} = \mathbf{0} となります。

2. 分散の比較
次に、\tilde{\boldsymbol{\beta}} の分散共分散行列を求めます。

ここで \mathbf{C}\mathbf{C}^\top を展開します。

不偏性の条件 \mathbf{D}\mathbf{X} = \mathbf{0}（およびその転置 \mathbf{X}^\top\mathbf{D}^\top = \mathbf{0}）を用いると式が整理されます。

したがって、\tilde{\boldsymbol{\beta}} の分散は以下のようになります。

任意の行列 \mathbf{D} に対して \mathbf{D}\mathbf{D}^\top は半正定値行列となるため、\text{Var}(\tilde{\boldsymbol{\beta}} \mid \mathbf{X}) \succeq \text{Var}(\hat{\boldsymbol{\beta}} \mid \mathbf{X}) が成立します。
以上により、OLS推定量 \hat{\boldsymbol{\beta}} は線形不偏推定量の中で最小の分散を持つ（BLUEである）ことが証明されました。

仮定に対する違反とその影響

Level 2の分析者には、各仮定が満たされなかった場合に、「どの性質が失われるのか」を明確に区別できるようになっていることが求められます。

実務上、A2（外生性）の仮定に対する違反が多く、深刻化しやすいです。これは点推定自体が真のパラメータから乖離してしまうためです（欠落変数バイアスや同時決定などの問題）。
一方、A3（等分散性）の仮定が満たされない場合、点推定の不偏性・一致性は維持されます。この場合、OLS推定量自体を捨てるのではなく、不均一分散に頑健な標準誤差（ロバスト標準誤差）を用いて標準誤差の計算式を修正するアプローチが一般的です。

1. 線形性の仮定が崩れた場合（モデルの誤特定）

データの真の姿が非線形であるにもかかわらず、単純な直線を当てはめてしまうと、係数の解釈が歪み、予測精度も低下します。この場合、OLSの枠組み（パラメータについて線形）を維持しつつ、変数の表現を工夫するか、あるいは別のアルゴリズムに頼るかの選択になります。

変数変換と多項式の追加

実務で最も頻繁に行われるアプローチです。散布図や残差プロットを確認し、関係性が直線でない場合は説明変数を変換してモデルに投入します。

• 対数変換（\log(x)）: 経済データ（賃金、売上、人口など）で頻出します。「x が1%変化したときに y が何%変化するか（弾力性）」という解釈が可能になり、極端に大きな外れ値の影響を抑える効果もあります。
• 多項式の追加（x^2, x^3）: 「年齢と賃金（中年でピーク）」「気温と電力消費（夏と冬にピーク）」のように、U字型や逆U字型の関係が想定される場合、説明変数に x だけでなく x^2 を同時に投入します。パラメータ \boldsymbol{\beta} に対しては依然として線形であるため、OLSがそのまま適用できます。

ラムゼイのRESET検定（Ramsey RESET Test）

「自分の作ったモデルが適切に線形性を捉えられているか（重要な非線形項を見落としていないか）」を統計的にチェックする処世術です。

1. まず通常のモデルで回帰し、予測値 \hat{\mathbf{y}} を求めます。
2. 次に、その予測値の2乗や3乗（\hat{\mathbf{y}}^2, \hat{\mathbf{y}}^3）を新たな説明変数として元のモデルに追加し、再度回帰を行います。
3. もし追加した項の係数が統計的に有意であれば、「モデルにはまだ表現しきれていない非線形な構造が残っている（誤特定がある）」と判断し、変数変換などを再検討します。

機械学習手法への切り替え

分析の目的が「係数から因果関係を解釈すること」ではなく、「未知のデータを高精度で予測すること」である場合の手段です。
ランダムフォレストや勾配ブースティング（LightGBMなど）の決定木ベースのアルゴリズムは、特徴エンジニアリングを頑張らなくとも、データから非線形性や変数間の相互作用を自動的に学習してくれますし、解釈性を除けばニューラルネットワークは心強い味方です。

2. 外生性の仮定が崩れた場合（内生性問題）

説明変数と誤差項が相関してしまう内生性は、推定量の不偏性・一致性を破綻させる最も深刻な問題です。これを解決するには、相関の根源（交絡や逆因果）を断ち切るアプローチが求められます。

操作変数法（IV）と2段階最小二乗法（2SLS）

内生性を持つ説明変数（X）の中から、「誤差項（見えない要因）と相関している部分」を取りいて回帰を行う分析方法です。

これには、以下の2条件を満たす 操作変数（Z） を見つけてくる必要があります。

1. 関連性: 操作変数 Z は、問題の説明変数 X と強い相関を持つ。
2. 除外制約: 操作変数 Z は、誤差項 \boldsymbol{\varepsilon} とは無相関であり、X を通じてのみ目的変数 Y に影響を与える。

【2SLSの手順】

• 第1段階: X を Z で回帰し、X の予測値 \hat{X} を求めます。この \hat{X} は Z によって作られているため、誤差項とは無相関になります。
• 第2段階: 本来の目的変数 Y を、生の X の代わりに \hat{X} で回帰します。これにより、内生性バイアスを取り除いた真の因果効果 \boldsymbol{\beta} を推定できます。

パネルデータ分析（固定効果モデル）

同一の個人や企業を、複数期間にわたって観測したデータ（パネルデータ）がある場合に使える処世術です。
欠落変数バイアスの多くは、「個人の持って生まれた能力」や「企業の立地」など、観測できないが時間的に変化しない要因（個体固有効果）によって引き起こされます。
固定効果モデルでは、各個体専用のダミー変数を投入する（あるいは個体ごとの平均値を引き算して消去する）ことで、これらの「見えない要因」を丸ごと吸収・排除します。これにより、時間とともに変化する変数の純粋な効果だけを抽出できます。

3. 等分散性・無相関の仮定が崩れた場合への処世術

この仮定が崩れても、OLS推定量の係数 \hat{\boldsymbol{\beta}} 自体は（不偏なので）使い物になります。問題は、OLSが BLUE（最小分散）ではなくなる という点です。これを効率性の喪失と呼びます。誤差項に不均一分散がある場合、誤差が小さいサンプルの情報には重みを置き、誤差が大きいサンプルの情報は軽めに扱った方が、よりブレの少ない（分散の小さい）推定が可能になります。OLSはすべてのデータを一律に扱うため、この情報の最適化ができず、結果として推定値が「たまたまの外れ値」に振り回されやすくなってしまうのです。また、この効率性の喪失だけでなく、「標準誤差の計算式」が狂うため、t検定の結果なども信頼できなくなります。

不均一分散下での効率性のロス

実際に、等分散性の仮定が崩れ、誤差項の分散共分散行列が \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \boldsymbol{\Omega} （ \boldsymbol{\Omega} \neq \sigma^2 \mathbf{I}）となった場合を確認してみます。

このとき、通常のOLS推定量 \hat{\boldsymbol{\beta}}_{\text{OLS}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} の「真の分散」は、以下のように計算されます^[8]。

一方、この不均一分散の構造 \boldsymbol{\Omega} を完全に把握し、それをデータに組み込んで最適化した推定量が 一般化最小二乗法（GLS） です。GLS推定量 \hat{\boldsymbol{\beta}}_{\text{GLS}} の分散は次のようにスッキリした形になります。

エイトケンの定理（ガウス=マルコフ定理の一般化）により、任意の線形推定量の中でGLSが最小分散を持つため、両者の分散共分散行列の間には必ず以下の大小関係が成立します。

どれくらい効率が悪くなるのか？

この分散の差（＝無駄にしてしまった効率性）の大きさは、説明変数 \mathbf{X} と誤差の分散の大きさが、どの程度連動しているかに依存して決まります。

例として「企業の規模（X）が大きいほど、利益の予測誤差の分散（\omega_i^2）も大きくなる」というデータを考えます。
最適な推定量であるGLSは、分散の逆数行列 \boldsymbol{\Omega}^{-1} を掛けることで「ノイズが大きすぎる大企業のデータは重視せず、ノイズの小さい中小企業のデータに重みを置いて直線を引く」という調整（加重最小二乗法）が行われます。
しかしOLSは、分散の大きさに限らず全てのデータを平等に扱うため、大企業データの巨大なノイズ（たまたまの大きな外れ値）に回帰直線が激しく引っ張られてしまいます。この連動が強いデータであるほど、OLSの分散はGLSよりも大きくなってしまいます（効率性の悪化）。

ロバスト標準誤差（不均一分散に頑健な標準誤差）

実務分析においてよく使う方法です。OLSのデフォルトの分散共分散行列の計算式 \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1} は、すべてのサンプルで分散 \sigma^2 が一定であることを前提としています。Huber-Whiteのロバスト標準誤差は、各サンプルの実際の残差の2乗（\hat{\varepsilon}_i^2）を利用して、不均一な分散の構造をデータから直接近似し、標準誤差を正しく補正します。多くの統計ソフト（Pythonのstatsmodelsなど）では、引数に cov_type='HC3' などを指定するだけで計算してくれます。

クラスター化ロバスト標準誤差

データの中に「同じクラスの生徒」「同じ都道府県の企業」「同一人物の複数年のデータ（パネルデータ）」といったグループ（クラスター）構造がある場合、同一クラスター内のサンプル同士は似たような未観測のショックを共有するため、誤差項に相関（空間的・時間的な自己相関）が生じます。

これを無視して通常の標準誤差を計算すると、実質的な情報量を過大評価してしまいます。たとえば「同じクラスの生徒40人」のデータは、完全に独立した40人分の情報を持っているわけではありません。これを独立だと仮定して分析すると、標準誤差が不当に小さく（t値が不当に大きく）計算され、本当は有意でないのに有意と判定してしまう危険性が非常に高くなります。これをモウルトン・バイアス（Moulton bias）と呼びます。

クラスター構造とサンドイッチ推定量

全データ数を N、クラスター数を G、各クラスター g に属するデータ数を N_g とします（g = 1, \dots, G）。
モデルをクラスターごとに分割して書くと、以下のようになります。

通常のOLSはすべての誤差項が無相関だと仮定しますが、クラスター構造がある場合、異なるクラスター間の誤差は無相関（独立）ですが、同じクラスター内の誤差には相関があると考えます。つまり、（データをクラスターで固めて並び替えれば）誤差項の全体の分散共分散行列 \boldsymbol{\Omega} は、クラスターごとのブロック対角行列になります。

この構造を許容した上で、正しい分散を推定するのがクラスター化ロバスト分散推定量（CRVE: Cluster-Robust Variance Estimator） です。通常のロバスト標準誤差（Huber-White推定量）の考え方を拡張し、個別の残差の2乗の代わりに、クラスターごとの残差ベクトル \hat{\boldsymbol{\varepsilon}}_g の外積を用いてサンドイッチの「具」^[9]の部分を構成します^[10]。

この式が意味するのは、クラスター内の任意の相関（不均一分散を含む）を許容しつつ、漸近的に正しい標準誤差を計算できるということです。ただし、この大数の法則が機能するためにはクラスターの数 G 自体が十分に大きいこと（一般的に30〜50以上が目安）が求められます。

【実務での処世術・実装方法】

実務の分析では非常によく使います。数式は複雑ですが、実装はオプション一つで解決してしまいます。

Python（statsmodels）の場合
fit() メソッドの引数に cov_type='cluster' を指定し、どの変数がクラスターを示すのかを渡します。

import statsmodels.formula.api as smf

# 'cluster_id'という列がクラスター（県や学校など）を表すとします
model = smf.ols('y ~ x1 + x2', data=df)

# クラスター化ロバスト標準誤差を指定して推定
result = model.fit(cov_type='cluster', cov_kwds={'groups': df['cluster_id']})

print(result.summary())

Rの場合
estimatr パッケージの lm_robust() を使うのが最も簡単です。

library(estimatr)

# clusters引数にクラスターのIDを指定するだけ
result <- lm_robust(y ~ x1 + x2, data = df, clusters = cluster_id)

summary(result)

どちらの言語でも、出力される係数の推定値（coef）自体は通常のOLSと全く同じ値になります。しかし、標準誤差（std err）がクラスター構造を加味して適切に大きめに補正され、それに伴いt値とp値がより慎重な値に変化していることが確認できるはずです。

一般化最小二乗法（GLS / FGLS）

ロバスト標準誤差が「OLSの計算結果をそのまま残し、標準誤差の計算式だけを後から補正する」という対症療法であるのに対し、GLS（一般化最小二乗法）は「データを根本から変換することで、OLSが再びBLUE（最小分散）となる性質を取り戻す」というアプローチです。

なぜデータを変換するのか？

不均一分散や自己相関が存在する場合、誤差項の分散共分散行列は \text{Var}(\boldsymbol{\varepsilon} \mid \mathbf{X}) = \boldsymbol{\Omega} となります（\boldsymbol{\Omega} は対角成分がバラバラだったり、非対角成分がゼロでなかったりする正定値対称行列です）。

ここで、行列の性質から \boldsymbol{\Omega} = \mathbf{P}\mathbf{P}^\top となるような正則行列 \mathbf{P} が存在します。計算をわかりやすくするため、この \mathbf{P} の逆行列を \boldsymbol{\Omega}^{-1/2} と表記しましょう（すなわち \boldsymbol{\Omega}^{-1/2}\boldsymbol{\Omega}(\boldsymbol{\Omega}^{-1/2})^\top = \mathbf{I} です）。

元の回帰モデル \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} の両辺に、左からこの \boldsymbol{\Omega}^{-1/2} を掛け合わせます。

これを新しい変数 \tilde{\mathbf{y}} = \tilde{\mathbf{X}}\boldsymbol{\beta} + \tilde{\boldsymbol{\varepsilon}} として置き換えます。すると、新しい誤差項 \tilde{\boldsymbol{\varepsilon}} の分散共分散行列は以下のようになります。

分散共分散行列が単位行列 \mathbf{I} になりました。これは変換後のデータにおいては、完全に等分散かつ無相関の仮定（A3）が満たされていることを意味します。
したがって、この変換されたデータに対して通常のOLSを適用すれば、ガウス=マルコフ定理が再び成立し、BLUE（最良線形不偏推定量）が得られます。

この推定量 \hat{\boldsymbol{\beta}}_{\text{GLS}} は以下のように導出されます。

この式の直観的な理解としては、\boldsymbol{\Omega}^{-1} を間に挟むことで、誤差（ノイズ）が大きいサンプルの情報は割り引いて軽く扱い、誤差が小さいサンプルの情報を重視して回帰直線を引くという加重最小二乗法（WLS）を行っています。これにより、OLSよりも圧倒的に効率的（分散が小さい）な推定が可能になります。

実行可能一般化最小二乗法（FGLS）による実務的アプローチ

理論上は完璧なGLSですが、実務においては真の分散共分散行列 \boldsymbol{\Omega} を知る由がないという根本的な問題があります。サンプルサイズが N のとき、推定すべき \boldsymbol{\Omega} の要素数は N(N+1)/2 個もあり、計算不可能です。

そこで実務では、\boldsymbol{\Omega} に「何らかの構造」を仮定し、データから \boldsymbol{\Omega} の推定値 \hat{\boldsymbol{\Omega}} を作って代入する 実行可能一般化最小二乗法（FGLS: Feasible GLS） が用いられます。

【FGLSの一般的なステップ】

1. とりあえずOLSを実行: 通常のOLSで \mathbf{y} を \mathbf{X} に回帰し、とりあえずの残差 \hat{\boldsymbol{\varepsilon}} を求めます。
2. 分散の構造をモデリング: 残差の2乗 \hat{\varepsilon}_i^2 が「何によって大きくなるのか」をモデル化します。例えば、ある変数 x_{1i} が大きくなるほど分散が広がると想定される場合、\log(\hat{\varepsilon}_i^2) = \alpha_0 + \alpha_1 x_{1i} + v_i のような補助回帰を行います。
3. 重みの作成: 補助回帰の予測値から、各サンプルの分散の推定値 \hat{\omega}_i^2 を計算し、これを用いて推定された分散共分散行列 \hat{\boldsymbol{\Omega}} を構築します。
4. FGLSの実行: 推定した \hat{\boldsymbol{\Omega}} をGLSの公式に代入し、最終的な係数 \hat{\boldsymbol{\beta}}_{\text{FGLS}} = (\mathbf{X}^\top\hat{\boldsymbol{\Omega}}^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\hat{\boldsymbol{\Omega}}^{-1}\mathbf{y} を計算します。

大標本においては、分散構造のモデルが正しい場合、\hat{\boldsymbol{\Omega}} を用いたFGLS推定量は、真の \boldsymbol{\Omega} を用いたGLSと同じ漸近分布を持って、漸近的に効率的になることが証明されています。ただし、小標本では逆にOLSより精度が悪くなるリスクもあるため、データサイズと分散の構造に対する理解が求められる難しいアプローチともいえます。

def fgls_diagonal(y, X, z_vars=None, max_iter=1, tol=1e-6):
    if z_vars is None:
        z_vars = X
    
    beta_prev = None
    for it in range(max_iter):
        # Step 1: とりあえずOLSを実行し、残差を得る
        if it == 0:
            resid = sm.OLS(y, X).fit().resid
        else:
            resid = y - X @ beta_prev
        
        # Step 2: 分散の構造をモデリング
        #         log(残差²) を z_vars に補助回帰
        log_e2 = np.log(resid**2 + 1e-12)
        aux_res = sm.OLS(log_e2, z_vars).fit()
        
        # Step 3: 重みの作成
        #         補助回帰の予測値から ω̂_i² を計算 → Ω̂ を構築（対角なので重みベクトルで表現）
        omega_hat = np.exp(aux_res.fittedvalues)
        weights = 1.0 / omega_hat
        
        # Step 4: FGLSの実行
        #         β̂_FGLS = (Xᵀ Ω̂⁻¹ X)⁻¹ Xᵀ Ω̂⁻¹ y を解く
        W_sqrt = np.sqrt(weights)
        X_w = X * W_sqrt[:, None]
        y_w = y * W_sqrt
        beta_new = np.linalg.lstsq(X_w, y_w, rcond=None)[0]
        
        if beta_prev is not None and np.max(np.abs(beta_new - beta_prev)) < tol:
            break
        beta_prev = beta_new
    
    return beta_new, omega_hat

4. 列フルランクの仮定が崩れた場合（多重共線性）への処世術

完全な多重共線性は、データに由来する避けられない問題というよりも、分析者によるミスです。

冗長なダミー変数

完全な多重共線性の代表例として、ダミー変数の種類の数だけダミー変数を用意してしまうものがあります。
例えば、「春・夏・秋・冬」という4つの季節カテゴリがあるとき、4つすべてのダミー変数をモデルに入れると、それらの和が常に「1（定数項と同じ）」になり、完全な線形従属が発生して逆行列が計算できなくなります。

これの対処は簡単で、カテゴリが N 個ある場合、モデルに投入するダミー変数は必ず N-1 個にします。除外した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までの議論は、不偏推定量（\mathbb{E}[\hat{\beta}] = \beta）のクラスにおいて、いかに分散を最小化するかという最良線形不偏推定量（BLUE）に焦点を当ててきました。しかし、実務的な予測においては、不偏性に固執することが必ずしも全体の予測誤差を最小化するとは限りません。

Level 3では、回帰の本質を「条件付き期待値の推定」と再定義し、バイアスとバリアンスのトレードオフを能動的に制御する手法群を扱います。

回帰の上位概念

統計的学習の枠組みにおいて、回帰問題は「入力 X が与えられたときのアウトカム Y の最良の予測値 f(X) を求めること」と考えられます。ここで、予測誤差の評価指標として平均二乗誤差（MSE）を採用すると、以下の最小化問題に帰着します。

この解は、理論的には条件付き期待値 f(x) = \mathbb{E}[Y \mid X = x] となります。

第1項はデータの本質的なノイズであり、g(X) に依存しません。第2項を最小化するには g(x) = \mathbb{E}[Y \mid X = x] とすればよいことがわかります。線形回帰は、この f(X) に \mathbf{X}\boldsymbol{\beta} という線形関数族の制約を置いた特殊ケースに相当します。

正則化：Ridge、Lasso

高次元データ（p \approx n）や多重共線性がある場合、行列 \mathbf{X}^\top\mathbf{X} の固有値が非常に小さくなり、その逆行列を用いるOLS推定量 \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y} の分散は爆発的に増大します。

正則化（Regularization）は、目的関数に残差平方和（RSS）以外のペナルティ項を加えることで、この分散を抑制します。

Ridge回帰（L2正則化）

係数の二乗和をペナルティとして加えます。

正規方程式を解くと、以下の解が導かれます。

ここで \lambda \mathbf{I} が加わることにより、\mathbf{X}^\top\mathbf{X} が特異に近い場合でも逆行列の計算が安定します。数理的には、推定量を 0 の方向へ縮小（Shrinkage）させる効果を持ち、バリアンスを小さくします（が、もちろんバイアスは増えます）。

Lasso回帰（L1正則化）

係数の絶対値の和をペナルティとして加えます。

L1ノルムの幾何学的性質により、多くの係数が正確に 0 になる「スパース性」が生じます。これは予測精度の向上と同時に、自動的な変数選択としての機能を果たします。なお、Ridgeと異なり解析的に解けないので数値的に解くことになります。

正則化の解釈：ベイズ推定との関連

正則化はベイズ統計学における最大事後確率（MAP）推定として解釈できます。

• Ridge: 係数 \boldsymbol{\beta} に平均 0 のガウス分布を事前分布として仮定した場合のMAP推定。
• Lasso: 係数 \boldsymbol{\beta} にラプラス分布（二重指数分布）を事前分布として仮定した場合のMAP推定。

線形回帰 y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i, \varepsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2) において、事後分布はベイズの定理より、

MAP 推定は事後分布を最大化する点なので、対数をとって、

尤度項はガウス誤差の仮定から、

したがって MAP 推定は次の最小化問題と等価になります。

Ridge の場合：\beta_j \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \tau^2) を仮定すると、

となり、\lambda = \sigma^2 / \tau^2 とおけば Ridge 回帰の目的関数

と一致します。

Lasso の場合：\beta_j \overset{\text{i.i.d.}}{\sim} \text{Laplace}(0, b)、すなわち p(\beta_j) = \frac{1}{2b} \exp(-|\beta_j| / b) を仮定すると、

となり、\lambda = 2\sigma^2 / b とおけば Lasso の目的関数

と一致します。

事前分布の裾の重さと原点での尖りが、推定量の性質を決めています。ガウス事前は原点で滑らかなので係数を 0 方向に「縮小」するだけですが、ラプラス事前は原点で微分不可能な尖りを持つため、解が厳密に \beta_j = 0 となる点に張り付き、スパース性が生じます。

バイアス-バリアンス分解

正則化がなぜ有効かを理解するための概念が、バイアス-バリアンス分解です。特定の点 x_0 における期待予測誤差は、以下の3つの成分に分解できます。

OLSは不偏推定量であるため、バイアスは 0 ですが、モデルが複雑であったりデータが少なかったりするとバリアンスが非常に大きくなります。
正則化は、\lambda を調整することで意図的にバイアスを導入し、それ以上にバリアンスを削減する手法です。その結果、全体のMSE（\text{Bias}^2 + \text{Variance}）が最小化されるポイントを探るのが Level 3 です。

非線形への一般化

線形モデルの枠組みを保ちつつ、データの非線形な関係を捉えるための拡張手法です。

多項式回帰とスプライン回帰

説明変数を \phi(x) = (x, x^2, \ldots, x^k) と高次化する多項式回帰は、特定の領域での極端な挙動に弱いです（ルンゲの現象）。
これに対し、スプライン回帰は領域を複数の区間に分割し、接続点（ノット）で連続かつ滑らか（一般に2次微分まで連続）になるよう低次多項式を当てはめます。

多項式回帰と自然キュービックスプラインの比較。高次の多項式回帰（赤）は学習データの範囲外で予測が極端に暴れる（ルンゲ現象）のに対し、自然スプライン（緑）はノット間で滑らかに結合しつつ、境界外では線形に保たれるためバリアンスが抑えられています。

特にB-スプラインや自然スプラインは、計算の安定性と解釈性のバランスに優れており、非線形なトレンドを捉えるための標準的な手法です。

一般化加法モデル（GAM）

複数の説明変数がある場合、各変数ごとの非線形効果を足し合わせた形でモデル化します。

各 f_j は平滑化スプラインなどの非パラメトリックな手法で推定されます。「変数間の加法性」という制約を置くことで、高次元データにおいても非線形性を維持しつつ、次元の呪いを回避し、各変数の寄与を視覚的に解釈することを可能にします。

一般化加法モデル（GAM）による成分プロットの概念図。各変数の非線形な影響（U字型、対数型など）を個別の関数 f_j(x_j) として抽出し、視覚的に解釈できます。X軸上のラグプロット（短い縦線）はデータの密度を表し、データがまばらな領域での極端な推定が起きていないか確認できます。

カーネル回帰

特定の関数形を仮定せず、予測したい点 x の近傍にあるデータ点 (x_i, y_i) の重み付き平均をとる手法です。Nadaraya-Watson推定量は次のように定義されます。

ここで K_h はカーネル関数（窓関数）であり、x に近い点ほど大きな重みを付与します。これは条件付き期待値 \mathbb{E}[Y|X] を直接的に推定しようとするノンパラメトリック手法の代表例です。

Level 3のまとめ

このレベルに達した分析者は、不偏性に捉われず、予測の高性能化を目指せます。回帰の目的が「未知のデータをうまく予測する」ことに力点を置き、あえてバイアスを導入してバリアンスを削るというトレードオフが視野に入ってきます。

• 回帰の上位概念が「条件付き期待値 \mathbb{E}[Y \mid X] の推定」であることを理解し、線形回帰がその一特殊ケースだと位置付けられるか？
• Ridge・Lassoを単なる「過学習対策」ではなく、ベイズ的にはガウス事前・ラプラス事前を置いたMAP推定として説明できるか？
• バイアス-バリアンス分解の各項を式で書き下し、\lambda を動かしたときの各項の挙動を語れるか？
• スプライン・GAM・カーネル回帰など、線形性の枠を保ったまま非線形性を取り込む手法を、目的に応じて使い分けられるか？

（自由度調整済み）決定係数 R^2 が高いか低いかに一喜一憂せず、モデルの複雑さをどこで制御し、いかにして汎化誤差を小さくするかという予測のトレードオフと向き合えれば、Level 3です。

さて、ここまでの議論は「相関」の話であり、X を動かしたら Y がどうなるかという「因果」の話はまだ登場していません^[11]。

Level 4: 因果推論へ

Level 3までで、私たちは「手元のデータに対して、予測するための妥当な線を引く方法」を習得しました。次の壁は、「回帰係数が有意であること、予測性能が高いこと」と「変数 X を動かせば Y が変わること」は、全く別の話だということです。

Level4では、予測モデルとしての回帰を、因果推論のためのツールとして見ていきます。

観測データの限界と内生性

私たちが手にするデータの多くは、ランダム化比較試験（RCT）によって得られたものではなく、単なる観測データです。観測データにおいて、Level2で学んだガウス＝マルコフの仮定 A2（外生性：\mathbb{E}[\varepsilon \mid \mathbf{X}] = 0）が満たされないことが多々あります。

この仮定が崩れた状態を 内生性（Endogeneity） があると言います。内生性があるとき、OLS推定量は不偏性も一致性も失い、因果推論が適切にできなくなります。

内生性を引き起こす3パターン

1. 欠落変数バイアス（OVB）: X と Y の両方に影響を与える「共通の要因」をモデルに入れられていない。
2. 同時決定（逆の因果）: X が Y を決めるだけでなく、Y も X を決めている（例：価格と需要）。
3. 測定誤差: 説明変数 X 自体に誤差が含まれており、それが誤差項と相関してしまう。

欠落変数バイアス（OVB）

まずは説明変数が欠落した場合（必要な変数を入力できていない場合）、推定した回帰係数が信頼できなくなることを確認します。

真のモデルが以下の2変数モデルであるとします。

ここで、分析者が x_2 の存在を知らず（あるいはデータがなく）、x_1 だけで回帰を行ったとします。

このとき、推定される係数 \hat{\alpha}_1 は、以下のようになります。

ここで、手元のデータ y は真のモデル y = \beta_1 x_1 + \beta_2 x_2 + \varepsilon によって生み出されているので、これを \hat{\alpha}_1 に代入します。

最後に、期待値はをとると以下のようになります（ここで「誤差項 \varepsilon は x_1 とは無相関である」という外生性の仮定から、\mathbb{E}[\text{Cov}(x_1, \varepsilon)] = 0 を使っています）。

右辺の第2項が欠落変数バイアスです。

• \beta_2 \neq 0（入れ忘れた変数 x_2 が y に影響を与える）
• \text{Cov}(x_1, x_2) \neq 0（入れ忘れた変数 x_2 が x_1 と相関している）

この2条件が揃ったとき、\hat{\alpha}_1 は真の因果効果 \beta_1 から離れていきます。例えば、「教育が賃金に与える効果」を測る際、「能力」を入れ忘れると、教育の係数は「純粋な教育の効果」＋「能力が高い人ほど教育を受けるという相関」が混ざったものになってしまうのです。

Frisch-Waugh-Lovell (FWL) 定理

次に、回帰分析において変数を追加する（コントロールする）と、因果推論に近づくのかを見ていきます。FWL定理 です。

FWL定理
重回帰モデル y = \mathbf{X}_1 \boldsymbol{\beta}_1 + \mathbf{X}_2 \boldsymbol{\beta}_2 + \varepsilon における \hat{\boldsymbol{\beta}}_1 は、以下の3ステップで得られる推定値と厳密に一致する。

1. y を \mathbf{X}_2 で回帰し、その残差 \tilde{y} を求める。
2. \mathbf{X}_1 を \mathbf{X}_2 で回帰し、その残差 \tilde{\mathbf{X}}_1 を求める。
3. \tilde{y} を \tilde{\mathbf{X}}_1 で回帰する。

この定理の主張は、「コントロールする」とは、悪影響を及ぼす変数の影響を事前に削除することです。

このプロセスを、情報（分散）の重なりを示すベン図で視覚的に理解してみましょう。

1. Step 1 & 2の意味（フィルター）
   y と x_1 の中には、それぞれ交絡因子 x_2 によって作られた「不純物（x_2 と相関している部分）」が含まれています。ここで x_2 で回帰して残差をとるという行為は、「x_2 の（線形射影）成分を除くフィルター」にかけることと同じです。これにより、x_2 の（線形成分の）影響を取り除いた部分（\tilde{y} と \tilde{x}_1）が抽出されます。
2. Step 3の意味（純粋な相関）
   重回帰分析における x_1 の係数とは、この「x_2 の影響を完全に削ぎ落とした純粋な部分」の相関を見ていることになります。

FWL定理を知っていると、重回帰分析をすれば因果関係がわかるという安易な考えを持たなくなるでしょう。私たちが本当に知りたい因果を抽出するためには、「モデルに入れていない未観測の変数」の影響までも、残差から完璧に削ぎ落とさなければならないのですが、データが存在しない以上それは不可能です。

だからこそ因果推論では、単なる回帰分析を諦め、「自然実験」や「操作変数」といった特殊なアプローチ（次節の識別戦略）を用いて、この「純粋な部分」を頑張って抽出しようと試みるのです。

識別戦略（Identification Strategy）

ここで紹介する戦略はすべて、擬似的なランダム化を作り出すための工夫です。

1. 操作変数法（IV / 2SLS）

X とは相関するが、誤差項 \varepsilon とは相関しない第3の変数 Z（操作変数）を利用します。「Z を通じて動いた X の変動だけ」を取り出すことで、内生性をクレンジングします。

• 2段階最小二乗法 (2SLS): 第1段階で \hat{X} = \text{Proj}(X \mid Z) を作り、第2段階で y を \hat{X} に回帰します。

2. 固定効果モデル（Fixed Effects）

パネルデータを用いて、個体固有の「変わらない特性（性格、企業の社風など）」を、個体ごとの平均からの差分をとることで取り除きます。

これにより、観測不可能な時不変の交絡因子によるバイアスを取り除きます。

3. 差分の差分法（DID）

「政策が行われたグループ」と「行われなかったグループ」の、前後の変化の「差」を比較します。

共通のトレンドさえ仮定（平行トレンド仮定）できれば、説得力のある方法です。

4. 回帰不連続デザイン（RDD）

「テストが60点以上なら合格」といった閾値（カットオフ）を利用します。59点の人と61点の人は「ほぼ同じ」であるという仮定に基づき、カットオフ付近での y のジャンプを因果効果と見なします。

Level 4のまとめ

このレベルに達した分析者は、出力された回帰係数を安易に「効果」と解釈しません。

• この係数は内生性によって過大評価されていないか？
• どの変数をコントロールすべきで、どの変数を入れてはいけないのか（合流点バイアスなど）？
• そもそもこの効果は識別可能なのか？

こうした因果的な問いを持ち、バイアスを適切に取り除けられれば（取り除けないことを理解できていれば）、Level 4です。

Level 5: 超高次元

Level 3 では、Lasso や Ridge を「予測精度を上げるための正則化」として見ました。
Level 5 では一歩進んで、p が非常に大きい状況下で、どのように推論を成立させるかを考えます。

p > n での問題

線形回帰モデル

を考えます。ここで p > n のとき、行列 \mathbf X のランクは高々 n なので、

です。したがって \mathbf X^\top \mathbf{X} は特異となり、通常の OLS 公式

は使えません。

使えないというのは、最小二乗問題の解が存在しないのではなく、解が一意に定まらないのです。 実際、\hat{\boldsymbol{\beta}}_0 が最小二乗解なら、任意の \mathbf v \in \mathrm{Null}(\mathbf X) ^[12]に対して

なので、\hat{\boldsymbol{\beta}}_0 + \mathbf v も同じ予測値を与えます。 つまり、同じ \hat{\mathbf y} を与える係数ベクトルが無数に存在するのです。

このとき必要になるのが、「その中でどの解をもっともらしいと見なすか」という追加の最適化問題です。そこで重要になってくるのがスパース性（sparsity） です。

以後、必ずしも p > n とせず、単に p が非常に大きく、スパース性がある状況を考えます。

スパース性と Lasso

高次元では、説明変数の総数 p は巨大でも、実際に効いている変数の数 s は小さい、つまり

であると仮定することが多いです。これを スパース性 と呼びます。
ここで、 \boldsymbol{\beta}^\astを真の係数ベクトルで、意味のない変数にかかる係数は0になっています。また、\mathrm{supp}(\boldsymbol{\beta}^\ast)は\boldsymbol{\beta}^\astの中で0ではない成分の番号の集合を表します。

この仮定のもとで代表的なのが Lasso です。

L^1 ペナルティ(\|\boldsymbol{\beta}\|_1 = \sum_j |\beta_j| は、多くの係数をちょうど 0 に押し込む性質があり、無数にある最小二乗解の中から、疎なものを選ぶという役割を果たします。

Level 3 ではこれを「予測のための正則化」と見ましたが、Level 5 ではさらに一歩進めて、どの条件のもとで、Lasso が真の \boldsymbol{\beta}^\ast に近づくのかを考えます。

1. 変数推定と変数選択は別物

Lasso について議論するとき、最初に区別すべきは次の 2 つです。

• (A) 推定された係数 \hat\beta が真の係数 \beta^\ast に近いか?
• (B) 推定されたサポート \hat S = \operatorname{supp}(\hat\beta) が真のサポート S = \operatorname{supp}(\beta^\ast) と一致するか?

(A) は予測・推定誤差の問題、(B) は変数選択の一貫性 (selection consistency) の問題で、要求される条件のレベルが全く違います。

1-1. 推定誤差が小さくなるための条件: Restricted Eigenvalue 条件^[13]

Lasso の解 \hat\beta_{\mathrm{Lasso}} は、最適化問題

を解いて得られます。p > n の高次元では、グラム行列 \hat\Sigma = \mathbf X^\top \mathbf X / n が退化していて最小固有値は 0 になるため、解が無数にある状態です。この中無数にある解の中から、何かしらの基準で最適なものを選ばねばなりません。

Lasso の最適性条件 (KKT 条件)

Lasso の解 \hat\beta_{\mathrm{Lasso}} がどのような点であるかを記述する条件として、KKT 条件 (Karush–Kuhn–Tucker 条件) を考えます。

【KKT 条件とは】

凸最適化問題 \min_x f(x) の最適解 \hat x では、目的関数の勾配が 0 になります(\nabla f(\hat x) = 0)。これが滑らかな場合の最適性条件です。KKT 条件はこれを拡張したもので、(i) 制約付き最適化や (ii) 微分不可能な目的関数にも対応できる最適性の必要十分条件です。

Lasso の目的関数

は、\ell_1 ペナルティ項 \|\beta\|_1 = \sum_j |\beta_j| のせいで、いずれかの \beta_j = 0 となる点では微分不可能です。そのため通常の「勾配 = 0」の代わりに、劣微分 (subdifferential) を使った KKT 条件を考えます。j 番目の成分について劣微分を取り、\hat\beta = \hat\beta_{\mathrm{Lasso}} で 0 が劣微分集合に含まれることを要求すると、各 j = 1, \ldots, p について

が成り立ちます。ここで \mathbf x_j \in \mathbb R^n は設計行列 \mathbf X の j 列目 (特徴量 j の観測値ベクトル)、\mathbf y - \mathbf X \hat\beta は Lasso 解の残差です。

この式の意味を考えます。|s_j| \le 1 なので、j がどの変数であっても

が成立します。\frac{1}{n}\mathbf x_j^\top (\mathbf y - \mathbf X \hat\beta) は特徴量 \mathbf x_j と残差ベクトルとの (スケールされた) 内積であり、特徴量を標準化していれば残差と特徴量の標本相関に対応する量です。つまり KKT 条件は「Lasso 解 \hat\beta の残差は、どの特徴量とも相関の大きさが \lambda 以下になっている」ということを述べています。これは「Lasso が残差をできるだけ各特徴量と無相関にするように \hat\beta を選ぶ」という Lasso のアルゴリズム的解釈とも合致します。

コーン制約の導出

次に、誤差ベクトル \boldsymbol\delta = \hat\beta - \beta^\ast がどのような領域に存在しうるかを調べます。

\hat\beta が最適解であることから、F(\hat\beta) \le F(\beta^\ast)、すなわち

が成り立ちます。\mathbf y = \mathbf X \beta^\ast + \boldsymbol\varepsilon を代入して両辺を展開・整理すると、

が得られます。チューニングパラメータ \lambda を、ノイズの寄与 \frac{1}{n}\|\mathbf X^\top \boldsymbol\varepsilon\|_\infty よりも十分大きく取る (具体的には \lambda \ge \frac{2}{n}\|\mathbf X^\top \boldsymbol\varepsilon\|_\infty となるよう取る) と、ヘルダーの不等式と \ell_1 ノルムの分解 \|\beta\|_1 = \|\beta_S\|_1 + \|\beta_{S^c}\|_1 を使って

という basic inequality が導けます。

左辺は \|\mathbf X \boldsymbol\delta\|_2^2 \ge 0 なので非負です。したがって右辺も非負でなければならず、

を得ます。これが、本記事で扱う コーン制約 \mathcal C(S, 3) そのものです^[14]。

コーン制約の意味

この不等式は、誤差ベクトル \boldsymbol\delta = \hat\beta - \beta^\ast について

を主張しており、「Lasso の誤差は、真サポート方向の大きさの 高々 3 倍までしか、非サポート方向に広がらない」ことを意味します。これは \ell_1 ペナルティが「非サポート方向の係数を 0 に引っ張る」効果として直感的にも理解できます。

この結果は Lasso の誤差解析にとって重要です。p > n では設計行列 \mathbf X は \mathbb R^p 全体で \|\mathbf X \boldsymbol\delta\|_2^2 \gtrsim \|\boldsymbol\delta\|_2^2 を満たすことは不可能です(退化しているため)。しかし、誤差ベクトルがコーン \mathcal C(S, 3) の中にしか存在しない以上、コーンの内側で \|\mathbf X \boldsymbol\delta\|_2^2 \gtrsim \|\boldsymbol\delta\|_2^2 が成立すれば誤差解析には十分です。これを定式化したのが Restricted Eigenvalue 条件です^[15]。

Restricted Eigenvalue (RE) 条件

真のサポートを S、|S| = s とする。ある定数 \kappa > 0 が存在して、コーン

上のすべての \boldsymbol\delta \ne 0 で

が成り立ちます。コーンの広さを決める定数として、ここでは 3 を採用しています(よく使われる値で、5 を使う流儀もあります)。

1-2. RE 条件のもとでの誤差レート

RE 条件と、適切なチューニングパラメータ

のもとで、Lasso は典型的に

というレートで真の係数に近づくことが示されます。これはサンプル数の対数オーダーで次元 p を許容できるということです。

ここで強調すべきは、RE 条件は「Lasso の推定誤差を抑える条件」であり、変数選択の完全一致を保証する条件ではないということです。

1-3. 真の変数集合を当てるための条件: Irrepresentable Condition

\hat S = \mathrm{supp}(\hat\beta_{\mathrm{Lasso}}) が真のサポート S ^[16]を正確に回復する (selection consistency) ためには、より強い条件が必要です。代表的なのが Irrepresentable Condition (IC)^[17] です。グラム行列を \boldsymbol\Sigma = \mathbf X^\top \mathbf X / n と書き、サポート S と非サポート S^c への分割を考えると、

を要求します。直感的な理解としては、「本当は不要な変数が、本当に必要な変数を強く模倣できないこと」を要求している感じです。\boldsymbol\Sigma_{S^c S}\boldsymbol\Sigma_{SS}^{-1} は「サポート外の変数を、サポート内の変数の線形結合でどれだけ表現できるか」を測る量で、これが 1 に近づくほど Lasso は機能しにくくなります。

さらに、信号強度も十分大きくなければいけません。具体的には、

のような beta-min 条件 が必要です。これは「真に非ゼロな係数は、ノイズに紛れない程度に大きい」という条件です。

1-4. 強さの順序

二つの条件には強弱の関係があります。

つまり、相関が中程度まで強くなると、予測と推定はうまくいくが、変数選択は失敗するというレジームに入ります。これが「Lasso は予測には強いが、選ばれた変数の解釈には注意」という実務の鉄則の理論的背景です。

1-5. 数値実験で確認する

理屈だけだとピンとこないので、実際に確認してみましょう。equicorrelation 構造（全特徴量が一様に相関 \rho を持つ）の設計行列で、\rho を動かしながら推定誤差とサポート回復率を測定します。

import numpy as np
from sklearn.linear_model import Lasso

def simulate(rho, n=300, p=60, s=5, n_trials=80):
    rng = np.random.default_rng(0)
    # equicorrelation: 対角 1, それ以外 rho
    cov = np.full((p, p), rho)
    np.fill_diagonal(cov, 1.0)
    L = np.linalg.cholesky(cov + 1e-8 * np.eye(p))

    # 真の係数: 最初の s 個だけが 1.5
    beta_star = np.zeros(p)
    beta_star[:s] = 1.5
    support_true = set(np.where(beta_star != 0)[0])

    l2_errors, exact_match = [], []
    for t in range(n_trials):
        Z = rng.normal(size=(n, p))
        X = Z @ L.T
        X = (X - X.mean(axis=0)) / X.std(axis=0)
        y = X @ beta_star + rng.normal(scale=1.0, size=n)

        lam = 1.5 * np.sqrt(2 * np.log(p) / n)
        model = Lasso(alpha=lam, max_iter=30000).fit(X, y)
        l2_errors.append(np.linalg.norm(model.coef_ - beta_star))
        support_hat = set(np.where(np.abs(model.coef_) > 1e-6)[0])
        exact_match.append(int(support_hat == support_true))
    return np.mean(l2_errors), np.mean(exact_match)

複数の \rho で実行して結果をプロットすると、次のようになります。

相関 \rho を動かしたときの推定誤差 (青) とサポート完全回復率 (赤) の挙動。横軸は equicorrelation の相関 \rho、左縦軸は \|\hat\beta - \beta^\ast\|_2、右縦軸はサポート完全一致 \hat S = S となった試行の割合。背景色で 3 つのレジームに分けています。

この実験結果には 3 つのレジームがあり、それぞれ理論と対応しています。

(A) 低相関領域 (\rho \lesssim 0.1): RE も IC も成立
両方の指標が良好です。サポートは100% 完全回復できています。推定誤差は中程度の値ですが、これは\lambdaによる縮小バイアスが残っているためです。

(B) 中相関領域 (0.1 \lesssim \rho \lesssim 0.6): RE は成立、IC は不成立
ここが最も興味深いレジームです。サポート完全回復率は急落し、\rho \approx 0.4 で既にほぼ 0 になります — IC が破れて、相関の強い偽の変数が混入するからです。ところが推定誤差はむしろ減少しており、\rho \approx 0.4 あたりで最小を取ります^[18]。これが「選択は失敗しているが推定は成功している」というレジームで、Lasso と多重共線性に関する最も実務的に重要な現象です。

(C) 強相関領域 (\rho \gtrsim 0.6): RE も不成立
ここでは推定誤差も増加に転じ、\rho \to 1 では発散していきます。equicorrelation の相関行列の最小固有値は 1 - \rho なので、\rho \to 1 で設計行列は退化し、RE の定数 \kappa も小さくなって理論保証そのものが緩くなります。両方の保証が同時に弱まる領域です。

1-6. 実務的教訓

この実験から学べることをまとめます。

• 「Lasso は変数選択をしてくれる」は誤解: 変数間に相関があると IC が成り立たず、サポート回復は信頼できなくなります。
• 予測・推定だけが目的なら Lasso はそこそこ頑健: 中程度の相関までは推定誤差は大きく崩れません。
• 極端な多重共線性は両方を壊す: \rho \to 1 では RE まで成立しなくなるため、Elastic Net や次元削減など別の手段を検討すべきです。

この区別を理解していないと、「Lasso で選ばれた変数は重要な変数」などと誤った主張をする恐れがあります。

2. Post-Selection Inference の罠

Lasso に関して、実務でついやってしまいがちな罠があります。

Lasso で変数選択して、選ばれた変数だけで OLS を回し、普通の標準誤差で有意性検定する

これはそのままではNGです。

2-1. なぜダメなのか

理由は単純で、同じデータを使って「選ぶ」ことと「検定する」ことを両方やっているからです。選択イベント自体が \mathbf y に依存している以上、その後の OLS 推論で「最初から変数が固定されていた」かのように扱うと、選択の不確実性を無視してしまいます。

直感的には、「たまたま大きく見えた係数だけを拾って、そのあとで“やっぱり大きいですね”と検定している」状態です。これはいわゆる winner's curse ^[19]であり、p 値は過小評価され、信頼区間の被覆率もいい加減になります。

数式で書くと、通常の信頼区間 CI_j^{\text{naive}} について一般には

です。つまり、「選ばれた後」の条件付きでは、通常の 95% 信頼区間は 95% を保たないのです。

2-2. どう対処するか

大きく次の 3 つの方法があります。

上2つはイメージしやすいと思うので、次節では3 つ目のデバイアスド Lasso を詳説します。

3. デバイアスド Lasso

Lasso は L^1 ペナルティのために係数を 0 の方向に縮めます。この縮小は予測には有利ですが、パラメータ推定には不利です。そこで、縮小による一次のバイアスを打ち消す補正を加えるのがデバイアスド Lasso のアイデアです。

3-1. 補正

デバイアスド Lasso の推定量は

と書けます。ここで \hat{\boldsymbol\Theta} は \mathbf X^\top \mathbf X / n の逆行列の近似で、典型的には nodewise Lasso で構成します。

理解を深めるために展開すると:

第 1 項と第 2 項を合わせると \bigl(I - \hat{\boldsymbol\Theta} \hat\Sigma\bigr)(\hat\beta_{\mathrm{Lasso}} - \beta^\ast) となり、\hat{\boldsymbol\Theta} \approx \hat\Sigma^{-1} ならこれは小さくなります。残る第 3 項 \hat{\boldsymbol\Theta} \mathbf X^\top \boldsymbol\epsilon / n は中心極限定理で正規分布に近づきます。

3-2. 漸近正規性

適切な正則条件のもとで、各係数 j について

が成り立ちます。ここで \hat\tau_j は対応する漸近標準誤差です。これにより、高次元でも個別係数について信頼区間や仮説検定ができるようになります。

もちろん、これは無条件に成立する魔法ではありません。設計行列の条件、真の係数の疎性、さらには \boldsymbol\Theta 側の疎性など、いくつかの仮定が必要です。それでも、「高次元で正則化しつつ推論する」という方向性を切り拓いた重要な考え方です。

3-3. 実装と被覆率の確認

scikit-learn は直接サポートしていませんが、nodewise Lasso を使って自前で実装できます。

import numpy as np
from sklearn.linear_model import Lasso
from scipy import stats

def debiased_lasso_ci(X, y, j, lam, alpha=0.05):
    """j 番目の係数についてデバイアスド Lasso の 95% 信頼区間を返す"""
    n, p = X.shape
    # 通常の Lasso
    lasso = Lasso(alpha=lam, max_iter=30000).fit(X, y)
    beta_hat = lasso.coef_
    residual = y - X @ beta_hat
    sigma_hat = np.sqrt(np.sum(residual**2) / max(n - np.sum(beta_hat != 0), 1))

    # nodewise Lasso: X_j を他の列で回帰して "逆共分散の j 行" を作る
    X_minus_j = np.delete(X, j, axis=1)
    x_j = X[:, j]
    lam_node = np.sqrt(2 * np.log(p) / n)
    nodewise = Lasso(alpha=lam_node, max_iter=30000).fit(X_minus_j, x_j)
    gamma = nodewise.coef_
    z_j = x_j - X_minus_j @ gamma            # X_j の "予測できない部分"
    tau_sq = (x_j @ z_j) / n

    # デバイアス補正項
    correction_factor = z_j / (n * tau_sq)
    correction = correction_factor @ residual
    beta_db = beta_hat[j] + correction

    # 漸近分散から CI
    se = sigma_hat * np.sqrt(np.sum(z_j**2)) / (n * abs(tau_sq))
    z_alpha = stats.norm.ppf(1 - alpha / 2)
    return beta_db, beta_db - z_alpha * se, beta_db + z_alpha * se

Naive (Lasso → OLS → 通常の CI) と Debiased Lasso で、それぞれ 95% 信頼区間の被覆率と平均幅を比較した結果が次の図です。

左: 95% 信頼区間の被覆率。点線が公称値 95%、これに近いほど良い結果です。Naive (赤) はどの \rho でも公称95%を下回り、\rho が大きくなるにつれて 90% 付近まで低下します。一方 Debiased Lasso (緑) は \rho \le 0.6 で公称95% にほぼ届いており、\rho = 0.8 でも 92% と Naive より明確に良い被覆を保ちます。右: 信頼区間の平均幅。両者はほぼ同程度の幅で、Debiased は Naive より少しだけ広い程度です。つまり Debiased Lasso は「ほぼ同じ効率で、より正しい被覆」を実現していると言えます。

メッセージを 3 つに整理すると:

1. Naive は失敗: 95% を一貫して下回り、相関とともに悪化していることが確認できます。これが正に Post-Selection Inference の罠です。
2. Debiased Lasso は理論に近い: 中相関までは95% に近い被覆を保ちます。
3. 極端な相関には限界がある: 極端な高相関 (\rho = 0.8) ではデバイアスドの仮定 (\Theta の疎性など) も次第に弱くなり、被覆率は徐々に低下します。それでも Naive よりは常に良いという結果になっています。

4. Double / Debiased Machine Learning (DML)

最後に、デバイアスドの考え方を機械学習による因果効果推定にまで拡張した DML を紹介します。これは Chernozhukov et al. (2018) ^[20] で体系化された枠組みです。

4-1. 部分線形モデル

シンプルな部分線形モデルを考えます。

• Y: アウトカム
• D: 処置変数 (treatment)
• X: 高次元の共変量 (confounders)
• \theta_0: 知りたい因果効果
• g_0, m_0: 複雑な局外 (nuisance) パラメータ
• U, V: ノイズ (E[U|X, D] = 0, E[V|X] = 0)

ここで g_0(X) は、処置 D の効果を除いた後に X が Y に与える成分です。
一方、DML の partialling-out 表現では、次の2つの条件付き期待値を局外関数として推定します。

このとき、

であり、一般には g_0(X) \neq \ell_0(X) です。

\theta_0 が興味の対象で、g_0, m_0 は推論したくないけれど無視もできない、という構造です。

したがって、

となります。
つまり、Y と D の両方から X で説明できる部分を取り除き、残差同士を回帰すれば \theta_0 を推定できます。

4-2. なぜ素朴に ML を入れるだけではダメなのか

g_0(X) や m_0(X) をランダムフォレストやブースティングで柔軟に推定したくなります。しかし、単純にそれらの予測値を代入して

のように \theta_0 を推定すると、

• \hat g の正則化バイアス (ML 推定器は予測精度を上げるため一般にバイアスを含む)
• 同じデータで \hat g を学習・適用することによる過学習バイアス

が \hat\theta にそのまま伝播してしまいます。つまり、予測器としては優秀でも、推論対象 \theta_0 の \sqrt n一致性や漸近正規性が成立しなくなります。

4-3. DML の核心: 直交化と cross-fitting

DML はこの問題を 2 つのアイデアで解決します。

(1) Neyman 直交化

局外パラメータ \eta = (g, m) の推定誤差に一次のスコア関数を作ります。部分線形モデルでは

を使います。重要な性質は、\eta = (g_0, m_0) のまわりで偏微分すると

となること(Neyman 直交性)です。つまり g, m が少しズレても \theta の推定への影響が 一次では消える ように設計されています。

(2) Cross-fitting

データを K 個に分割し、ある fold の \hat g, \hat m は 別 fold で学習します。これにより、「同じデータで学習した複雑な予測器を、同じデータにそのまま当てる」ことによる過学習バイアスを取り除きます。

このとき、部分線形モデルの DML 推定量は実質的に

となります。これは見方を変えると、

X で説明できる部分を Y と D の両方から取り除き、残差どうしを回帰している

ことに相当します。この意味で DML は、Frisch–Waugh–Lovell (FWL) 定理の機械学習版と見なせます。

4-4. 理論的保証

十分条件の一つとして、局外パラメータの推定誤差が

を満たせば、\hat\theta_{\mathrm{DML}} は \sqrt n-一致かつ漸近正規になります。各推定器がそれぞれ n^{-1/4} 程度の収束レートを持てばこの条件は満たされるので、DML はかなり柔軟な機械学習器を許容します。これは線形モデルが必要とする o_p(1) レベルの一致性よりは厳しいですが、\sqrt n-一致性よりはずっと緩いという絶妙なバランスです。

4-5. 実装で確認する

実装は驚くほどシンプルで、本質は数十行で書けます。

import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import KFold

def naive_plugin(X, D, Y):
    """素朴な plug-in: 同じデータで g_hat を学習して Y から引く"""
    rf = RandomForestRegressor(n_estimators=50, max_depth=8, random_state=0).fit(X, Y)
    g_hat = rf.predict(X)               # ← 同じデータで予測 (過学習バイアス)
    return np.sum(D * (Y - g_hat)) / np.sum(D**2)

def dml_estimate(X, D, Y, K=5, seed=0):
    """cross-fitting + Neyman 直交スコア"""
    n = len(Y)
    kf = KFold(n_splits=K, shuffle=True, random_state=seed)

    Y_resid = np.zeros(n)
    D_resid = np.zeros(n)

    for train, test in kf.split(X):
        ell = RandomForestRegressor(
            n_estimators=50,
            max_depth=8,
            random_state=seed
        ).fit(X[train], Y[train])

        m = RandomForestRegressor(
            n_estimators=50,
            max_depth=8,
            random_state=seed
        ).fit(X[train], D[train])

        Y_resid[test] = Y[test] - ell.predict(X[test])
        D_resid[test] = D[test] - m.predict(X[test])

    return np.sum(D_resid * Y_resid) / np.sum(D_resid ** 2)

データ生成過程として g_0, m_0 に非線形な真の関数を入れ、Naive と DML で \theta_0 = 2.0 をどれだけ正確に推定できるか比較します。

Naive plug-in (赤) は \theta_0 = 2.0 から大きくバイアスを持つ (実験では \hat\theta \approx 0.43, バイアス -1.57) のに対し、DML (緑) は真値のすぐ近くに集中して分布する (実験では \hat\theta \approx 1.90, バイアス -0.10)。これが Neyman 直交化と cross-fitting の威力です。

差は劇的です。Naive plug-in は真値から大きく外れた値に集中しており、いくらサンプル数を増やしてもこのバイアスは消えません。一方 DML は真値の周りに集中しており、\sqrt n で詰めていける期待があります。

4-6. DML で何が嬉しいのか

DML は、

• 柔軟な機械学習で局外パラメータ g_0, m_0 を扱える（線形性などの強い仮定が要らない）
• 興味のあるパラメータ \theta_0 については 厳密な統計推論 （漸近正規性、信頼区間）ができる

というところがとても嬉しいです。これは因果推論を機械学習時代にアップデートする上で重要な枠組みです。R/Python ともに DoubleML というしっかりした実装パッケージがあるので、実務ではそちらの利用も検討してください。

https://docs.doubleml.org/

Level 5のまとめ

このレベルに達した分析者は、p が非常に大きい現代的なデータ状況において、「正則化すれば変数選択もできて推論もできる」といった淡い期待を持ちません。Lasso が何を保証し、何を保証しないのかを、条件のレベルで区別できることが求められます。

• 推定誤差を抑える条件（Restricted Eigenvalue）と、真の変数集合を当てる条件（Irrepresentable Condition + beta-min）の強弱関係を理解しているか？
• 「Lassoで変数を選んでからOLSで検定する」というナイーブな手続きがなぜダメなのか、Post-Selection Inferenceの罠を説明できるか？
• デバイアスドLassoが、L^1 ペナルティによる縮小バイアスを解析的に補正することで、高次元でも漸近正規性を回復する仕組みを語れるか？
• Double/Debiased Machine Learningが、Neyman直交化とcross-fittingという2つの仕掛けで、「機械学習の柔軟性」と「\sqrt{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.