賃金の決定要因を分析するために、就業者のデータで回帰分析を行ったとします。教育年数、経験年数、職種などを説明変数にして賃金を被説明変数とする回帰モデルです。しかし、ここに重大な落とし穴があります。このサンプルは母集団からの「ランダムサンプル」ではないのです。
就業するかどうかの意思決定自体が、賃金に影響する変数（学歴、スキル、育児の有無、配偶者の所得など）に依存しています。例えば、高い賃金を得られる見込みがある人ほど就業する傾向があります。このような場合、就業者のサンプルは母集団よりも「潜在的に高い賃金を得られる人」に偏っています。この偏りを無視して回帰分析を行うと、回帰係数の推定にバイアスが生じます。
このように「サンプルに入るかどうか」が分析対象の変数と系統的に関連する問題を選択バイアス（selection bias）と呼びます。ジェームズ・ヘックマン（James Heckman）はこの問題に対する画期的な補正法を提案し、2000年のノーベル経済学賞を受賞しました。
選択バイアスとヘックマン補正は以下の幅広い分野で重要です。
- 労働経済学: 賃金関数の推定（就業選択を考慮）。ヘックマン自身が女性の労働供給を分析した古典的研究
- 医学研究: 治療の効果分析（治療を受けるかどうかの選択が予後と相関）
- マーケティング: 商品購入者のみを対象にした満足度分析（購入選択バイアス）
- 教育: 大学進学者のみを対象にした教育効果の分析
- 金融: 融資承認者のみのデフォルト率分析（審査基準による選択）
本記事の内容
- 選択バイアスの数理的定式化と直感的理解
- 切断正規分布と条件付き期待値の導出
- 逆ミルズ比（Inverse Mills Ratio）の定義と性質
- ヘックマン2段階推定法の理論
- 識別条件と除外変数の役割
- 最尤推定法との比較
- トービットモデルとの関係
- 因果推論の枠組みからの再解釈
- Pythonシミュレーションと推定結果の可視化
前提知識
この記事を読む前に、以下の知識があると理解が深まります。
- 正規分布 — 正規分布の性質
- 交絡バイアスとその制御法 — バイアスの基本概念
- 最尤推定 — 最尤推定の基礎
選択バイアスの直感的理解
日常的なアナロジー
選択バイアスの直感をつかむために、大学のGPAと入学試験スコアの関係を考えましょう。
ある名門大学で「入学試験のスコアとGPAの間に負の相関がある」というデータが得られたとします。一見不思議ですが、これは選択バイアスの典型例です。
入学許可を得た学生は、「試験スコアが高い」か「推薦状が優れている」か、その両方の条件を満たしています。試験スコアがやや低い学生が入学しているとすれば、それは推薦状が特に優れていた（= GPA上昇に寄与する他の能力が高い）ためです。逆に、試験スコアが高い学生の中には、推薦状はそこそこで入学した人もいます。この選択プロセスにより、入学者のサンプル内では試験スコアとGPAに負の相関が生じるのです。
母集団全体（受験者全体）で見れば試験スコアとGPAは正の相関を持つかもしれないのに、選択されたサンプル（入学者）では相関が逆転する — これが選択バイアスの恐ろしさです。
3つのタイプの選択バイアス
選択バイアスは、データがサンプルに入るプロセスに応じていくつかのタイプに分類されます。
タイプ1: サンプルセレクション（Sample Selection）
最も一般的なタイプです。アウトカムが観測されるかどうかが、アウトカム自体と相関する変数に依存します。例: 就業者のみの賃金データ、生存者のみの健康データ。
タイプ2: 自己選択（Self-Selection）
個体が自らの判断で処置を選択する場合に生じます。例: 職業訓練プログラムに自主的に参加する人は、元々モチベーションが高い。
タイプ3: 生存者バイアス（Survivorship Bias）
分析時点まで「生存」したユニットのみがサンプルに含まれる場合に生じます。例: 長期間存続している企業のみの分析、戦争から帰還した航空機の損傷パターンの分析。
本記事では、主にタイプ1のサンプルセレクションを数学的に扱います。
数学的な定式化に進みましょう。
選択バイアスの数学的定式化
サンプル選択モデル
2つの潜在変数を含むモデルを考えます。
結果方程式（outcome equation）:
$$ \begin{equation} Y_i^* = \bm{X}_i^\top \bm{\beta} + \epsilon_i \end{equation} $$
選択方程式（selection equation）:
$$ \begin{equation} S_i^* = \bm{Z}_i^\top \bm{\gamma} + \nu_i \end{equation} $$
ここで: – $Y_i^*$: 潜在的なアウトカム（例: 潜在賃金） – $S_i^*$: 選択の潜在変数（例: 就業意欲の指標） – $\bm{X}_i$: 結果方程式の説明変数 – $\bm{Z}_i$: 選択方程式の説明変数（$\bm{X}_i$ を含みうる） – $\epsilon_i$, $\nu_i$: 誤差項
$Y_i^*$ が観測されるのは $S_i^* > 0$（すなわち選択された場合 $S_i = 1$）のときだけです。
$$ Y_i = \begin{cases} Y_i^* & \text{if } S_i^* > 0 \quad (S_i = 1) \\ \text{unobserved} & \text{if } S_i^* \leq 0 \quad (S_i = 0) \end{cases} $$
二変量正規性の仮定
誤差項が二変量正規分布に従うと仮定します。
$$ \begin{equation} \begin{pmatrix} \epsilon_i \\ \nu_i \end{pmatrix} \sim N\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma^2 & \rho\sigma \\ \rho\sigma & 1 \end{pmatrix}\right) \end{equation} $$
$\nu_i$ の分散は識別のために1に正規化しています（プロビットモデルの標準的な正規化）。$\rho = \text{Corr}(\epsilon_i, \nu_i)$ は2つの誤差項の相関係数であり、選択バイアスの大きさと方向を決定する鍵となるパラメータです。
選択バイアスの発生
$S_i = 1$（選択された）サンプルのみで $Y_i$ を $\bm{X}_i$ に回帰した場合の条件付き期待値を計算します。
$$ \begin{align} E[Y_i \mid S_i = 1, \bm{X}_i, \bm{Z}_i] &= \bm{X}_i^\top \bm{\beta} + E[\epsilon_i \mid \nu_i > -\bm{Z}_i^\top \bm{\gamma}] \\ &= \bm{X}_i^\top \bm{\beta} + \underbrace{E[\epsilon_i \mid \nu_i > -\bm{Z}_i^\top \bm{\gamma}]}_{\text{選択バイアス項} \neq 0} \end{align} $$
$\rho \neq 0$（$\epsilon_i$ と $\nu_i$ が相関）の場合、$E[\epsilon_i \mid \nu_i > -\bm{Z}_i^\top \bm{\gamma}] \neq 0$ です。つまり、選択されたサンプルの誤差項の条件付き期待値はゼロではなく、OLS推定にバイアスが生じます。
このバイアスは欠落変数バイアス（omitted variable bias）と見なすこともできます。選択バイアス項は $\bm{Z}_i$ の関数であり、これを省略して回帰を行うことで、$\bm{X}_i$ と相関する部分がバイアスとして現れるのです。
では、この選択バイアス項を具体的に計算するために、切断正規分布の性質を利用します。
切断正規分布と逆ミルズ比
切断正規分布の条件付き期待値
二変量正規分布の条件付き期待値の公式を使います。$(X, Y)$ が二変量正規分布に従い、$\text{Corr}(X, Y) = \rho$、$\text{Var}(X) = \sigma_X^2$、$\text{Var}(Y) = \sigma_Y^2$ のとき、
$$ E[X \mid Y > c] = \mu_X + \rho \frac{\sigma_X}{\sigma_Y} \cdot \frac{\phi\left(\frac{c – \mu_Y}{\sigma_Y}\right)}{\Phi\left(\frac{\mu_Y – c}{\sigma_Y}\right)} $$
これを $(\epsilon, \nu)$ に適用します。$\mu_\epsilon = \mu_\nu = 0$、$\sigma_\epsilon = \sigma$、$\sigma_\nu = 1$ として、$\nu > -\bm{Z}_i^\top \bm{\gamma}$ の条件を代入すると、
$$ \begin{align} E[\epsilon_i \mid \nu_i > -\bm{Z}_i^\top \bm{\gamma}] &= \rho \sigma \cdot \frac{\phi(-(-\bm{Z}_i^\top \bm{\gamma}))}{\Phi(-(-\bm{Z}_i^\top \bm{\gamma}))} \\ &= \rho \sigma \cdot \frac{\phi(\bm{Z}_i^\top \bm{\gamma})}{\Phi(\bm{Z}_i^\top \bm{\gamma})} \end{align} $$
最後の行では、$\phi$ の対称性（$\phi(-x) = \phi(x)$）を使いました。
逆ミルズ比の定義
上の式に現れる比を逆ミルズ比（Inverse Mills Ratio, IMR）と呼び、$\lambda(\cdot)$ で表記します。
$$ \begin{equation} \lambda(a) = \frac{\phi(a)}{\Phi(a)} \end{equation} $$
ここで $\phi(\cdot)$ は標準正規分布の密度関数、$\Phi(\cdot)$ は累積分布関数です。
逆ミルズ比の性質
逆ミルズ比 $\lambda(a)$ は以下の重要な性質を持ちます。
- 常に正: $\lambda(a) > 0$ for all $a$
- 単調減少: $a$ が増加すると $\lambda(a)$ は減少する
- 漸近的振る舞い: $a \to -\infty$ のとき $\lambda(a) \to +\infty$、$a \to +\infty$ のとき $\lambda(a) \to 0$
- 微分: $\lambda'(a) = -\lambda(a)[a + \lambda(a)]$
性質2と3の直感的理解: $a = \bm{Z}_i^\top \bm{\gamma}$ が大きい場合、選択確率 $\Phi(a)$ が1に近く、ほぼ全員が選択されるため、選択条件付けの影響は小さくなります。逆に $a$ が小さい（選択確率が低い）場合、選択されたサンプルは強くフィルタリングされているため、選択バイアスの影響が大きくなります。
条件付き期待値の導出の詳細
逆ミルズ比が出現する理由を、より詳しく導出しましょう。二変量正規分布 $(\epsilon, \nu)$ において、$\nu > -a$（$a = \bm{Z}_i^\top \bm{\gamma}$）という条件のもとでの $\epsilon$ の条件付き期待値を計算します。
まず、$\nu$ が与えられたときの $\epsilon$ の条件付き分布は、
$$ \begin{equation} \epsilon \mid \nu \sim N\left(\rho\sigma \cdot \nu, \, \sigma^2(1 – \rho^2)\right) \end{equation} $$
これは二変量正規分布の基本的な性質です。次に、$\nu > -a$ の条件のもとでの期待値を繰り返し期待値の法則で計算します。
$$ \begin{equation} E[\epsilon \mid \nu > -a] = E[E[\epsilon \mid \nu] \mid \nu > -a] = E[\rho\sigma\nu \mid \nu > -a] = \rho\sigma \cdot E[\nu \mid \nu > -a] \end{equation} $$
$\nu \sim N(0, 1)$ の切断正規分布の期待値は、
$$ \begin{equation} E[\nu \mid \nu > -a] = \frac{\phi(-a)}{\Phi(a)} = \frac{\phi(a)}{\Phi(a)} = \lambda(a) \end{equation} $$
最後の等号では $\phi$ の対称性 $\phi(-a) = \phi(a)$ を使いました。
したがって、
$$ \begin{equation} E[\epsilon \mid \nu > -a] = \rho\sigma \cdot \lambda(a) \end{equation} $$
が得られます。この導出から、選択バイアスの大きさが $\rho$（誤差の相関）、$\sigma$（結果方程式の誤差の標準偏差）、$\lambda(a)$（選択の厳しさ）の3つの要素で決まることが明確になります。
切断正規分布の条件付き分散
条件付き期待値だけでなく、条件付き分散も重要です。$\nu > -a$ の条件のもとでの $\epsilon$ の条件付き分散は、
$$ \begin{equation} \text{Var}(\epsilon \mid \nu > -a) = \sigma^2 \left[1 – \rho^2 \delta(a)\right] \end{equation} $$
ここで $\delta(a) = \lambda(a)[\lambda(a) – a]$ です。$0 < \delta(a) < 1$ であるため、選択条件付けにより分散は常に縮小します。直感的には、選択された個体は母集団全体よりも均質（アウトカムのばらつきが小さい）であるため、分散が小さくなるのです。
この不均一分散の存在がヘックマン2段階推定法の標準誤差修正を必要とする理由です。
これらの準備のもとで、ヘックマン補正の手順を導出できます。
ヘックマン2段階推定法
補正された結果方程式
選択バイアス項を代入すると、選択されたサンプルの条件付き期待値は次のようになります。
$$ \begin{equation} E[Y_i \mid S_i = 1, \bm{X}_i, \bm{Z}_i] = \bm{X}_i^\top \bm{\beta} + \rho\sigma \cdot \lambda(\bm{Z}_i^\top \bm{\gamma}) \end{equation} $$
これは、元の回帰方程式に逆ミルズ比 $\lambda(\bm{Z}_i^\top \bm{\gamma})$ を追加変数として含めることで、選択バイアスを補正できることを示しています。
しかし、$\bm{\gamma}$ は未知であるため、まず選択方程式から $\bm{\gamma}$ を推定する必要があります。これがヘックマン2段階推定法の「2段階」の意味です。
第1段階: プロビットモデルの推定
全サンプル（$S = 1$ と $S = 0$ の両方）を使って、選択方程式のプロビットモデルを最尤推定します。
$$ P(S_i = 1 \mid \bm{Z}_i) = \Phi(\bm{Z}_i^\top \bm{\gamma}) $$
プロビットモデルの対数尤度は次の通りです。
$$ \ell(\bm{\gamma}) = \sum_{i=1}^{n} [S_i \ln \Phi(\bm{Z}_i^\top \bm{\gamma}) + (1 – S_i) \ln(1 – \Phi(\bm{Z}_i^\top \bm{\gamma}))] $$
この最大化により $\hat{\bm{\gamma}}$ を得ます。
第2段階: 逆ミルズ比を追加したOLS
推定された $\hat{\bm{\gamma}}$ を使って逆ミルズ比を計算します。
$$ \hat{\lambda}_i = \lambda(\bm{Z}_i^\top \hat{\bm{\gamma}}) = \frac{\phi(\bm{Z}_i^\top \hat{\bm{\gamma}})}{\Phi(\bm{Z}_i^\top \hat{\bm{\gamma}})} $$
$S_i = 1$（選択された）サンプルのみで、次の回帰をOLSで推定します。
$$ \begin{equation} Y_i = \bm{X}_i^\top \bm{\beta} + \delta \hat{\lambda}_i + u_i \end{equation} $$
ここで $\delta = \rho\sigma$ であり、$u_i$ は新しい誤差項です。$\hat{\bm{\beta}}$ が選択バイアスを補正した結果方程式の回帰係数の推定量です。
選択バイアスの検定
$\hat{\delta}$ の有意性検定は、選択バイアスの存在の検定に対応します。
- $H_0: \delta = 0$（選択バイアスなし = $\rho = 0$）
- $H_1: \delta \neq 0$（選択バイアスあり = $\rho \neq 0$）
$\hat{\delta}$ が統計的に有意であれば、選択バイアスが存在する証拠であり、ヘックマン補正が必要であることを意味します。
標準誤差の修正
ヘックマン2段階推定法の第2段階のOLS標準誤差は、2つの理由で正しくありません。
- 逆ミルズ比が推定量: $\hat{\lambda}_i$ は第1段階の推定値 $\hat{\bm{\gamma}}$ に依存しているため、通常のOLS標準誤差は第1段階の推定誤差を考慮していません。
- 不均一分散: $u_i$ の分散は $\bm{Z}_i$ に依存するため、均一分散の仮定が成り立ちません。
修正された標準誤差は、$\hat{\bm{\gamma}}$ の推定誤差を考慮したサンドイッチ推定量で計算できます。具体的な公式はやや複雑ですが、多くの統計ソフトウェアがこの修正を自動的に行います。
識別条件と除外変数
除外変数の役割
ヘックマン2段階推定法が適切に機能するためには、選択方程式に含まれるが結果方程式には含まれない変数——除外変数（exclusion restriction）——が存在することが理論的に望ましいです。
除外変数 $Z_{\text{excl}}$ は以下の2つの条件を満たす必要があります。
- $Z_{\text{excl}}$ は選択 $S$ に影響する（$Z_{\text{excl}} \to S$）
- $Z_{\text{excl}}$ はアウトカム $Y^*$ に直接影響しない（$Z_{\text{excl}} \not\to Y^*$）
例えば、賃金関数の推定において: – 学歴、経験年数: $\bm{X}$ と $\bm{Z}$ の両方に含まれる（賃金にも就業にも影響） – 配偶者の所得、子供の数: $\bm{Z}$ にのみ含まれる除外変数の候補（就業決定に影響するが、潜在賃金には直接影響しないと仮定）
除外変数がない場合の識別
技術的には、除外変数がなくても（$\bm{Z} = \bm{X}$）ヘックマンモデルは識別可能です。なぜなら、逆ミルズ比 $\lambda(\bm{X}_i^\top \bm{\gamma})$ は $\bm{X}_i$ の非線形関数であり、この非線形性が $\bm{\beta}$ と $\delta$ の識別を提供するからです。
しかし、この非線形性のみに依存する識別は非常に脆弱です。逆ミルズ比は $\bm{X}_i^\top \bm{\gamma}$ の中間的な範囲ではほぼ線形に近い振る舞いをするため、$\bm{X}_i$ の線形関数と区別がつきにくくなります。結果として、多重共線性の問題が生じ、推定値が不安定になります。
したがって、実践上は少なくとも1つの除外変数を含めることが強く推奨されます。
除外変数の選択は、操作変数法における操作変数の選択と同様の困難さがあります。「選択に影響するが結果に直接影響しない」という条件はドメイン知識に基づいて判断する必要があり、データから検証することは一般に困難です。
最尤推定法との比較
完全情報最尤推定（FIML）
ヘックマン2段階推定法の代替として、完全情報最尤推定法（Full Information Maximum Likelihood, FIML）があります。FIMLは選択方程式と結果方程式を同時に推定します。
観測されたデータの対数尤度は次の通りです。
$$ \begin{align} \ell(\bm{\beta}, \bm{\gamma}, \sigma, \rho) &= \sum_{i: S_i = 1} \ln f(Y_i, S_i = 1 \mid \bm{X}_i, \bm{Z}_i) + \sum_{i: S_i = 0} \ln P(S_i = 0 \mid \bm{Z}_i) \end{align} $$
二変量正規性の仮定のもとで、$S_i = 1$ の観測の同時密度は次のように書けます。
$$ f(Y_i, S_i = 1 \mid \bm{X}_i, \bm{Z}_i) = \frac{1}{\sigma}\phi\left(\frac{Y_i – \bm{X}_i^\top\bm{\beta}}{\sigma}\right) \Phi\left(\frac{\bm{Z}_i^\top\bm{\gamma} + \rho(Y_i – \bm{X}_i^\top\bm{\beta})/\sigma}{\sqrt{1 – \rho^2}}\right) $$
2段階推定法とFIMLの比較
| 特徴 | 2段階推定法 | FIML |
|---|---|---|
| 効率性 | やや非効率 | 漸近的に効率的 |
| 計算の容易さ | 容易（プロビット + OLS） | 数値最適化が必要 |
| 頑健性 | 正規性仮定に敏感 | 正規性仮定に敏感 |
| 標準誤差 | 修正が必要 | 自然に正しい |
| 収束 | 常に収束 | 収束しない場合がある |
2段階推定法は計算が容易で常に収束する利点がありますが、FIMLよりも推定効率が低くなります。一方、FIMLは漸近的に効率的ですが、数値最適化の収束が保証されないことがあります。
半パラメトリック手法
ヘックマンモデルの正規性仮定への懸念
ヘックマン2段階推定法とFIMLはともに二変量正規性の仮定に依存します。この仮定が破れた場合、推定量にバイアスが生じます。実際の経済データでは、賃金分布が右に歪んでいたり、選択メカニズムが非対称であったりすることがあり、正規性の仮定が疑わしいケースは少なくありません。
セミパラメトリック推定法
正規性仮定を緩和するセミパラメトリック推定法がいくつか提案されています。
Powell（1987）の対称切断回帰推定量（symmetrically trimmed least squares, STLS）: トービットモデルの文脈で正規性仮定を緩和した推定量です。
Newey（1988）のシリーズ推定法: 逆ミルズ比の代わりに、選択確率のべき級数展開を追加変数として使います。選択方程式をノンパラメトリックに推定し、$\hat{P}(S_i = 1 \mid \bm{Z}_i)$ のべき級数
$$ \begin{equation} Y_i = \bm{X}_i^\top \bm{\beta} + \sum_{j=1}^{J} \delta_j \hat{P}_i^j + u_i \end{equation} $$
をOLSで推定します。逆ミルズ比という特定の関数形を仮定せず、一般的な制御関数をべき級数で近似するのです。
Heckman and Honoré（1990）のセミパラメトリック識別: 除外変数の存在と結合分布の尾の性質に関する弱い条件のもとで、選択モデルが識別可能であることを示しました。
Klein and Spady（1993）のセミパラメトリック推定: 選択方程式をセミパラメトリック単一指数モデルで推定する方法です。
これらのセミパラメトリック手法は正規性仮定への依存を軽減しますが、一般に推定の効率が落ちる（分散が大きくなる）というトレードオフがあります。
実用上のアドバイス
実務的には、以下のアプローチが推奨されます。
- まずヘックマン2段階推定法で分析を行う
- 正規性の仮定が疑わしい場合は、半パラメトリック手法でも推定し、結果が大きく変わらないかを確認する（感度分析）
- 結果が正規性仮定に敏感である場合は、半パラメトリック手法の結果を主たる推定として報告する
トービットモデルとの関係
トービットモデル
トービットモデル（Tobit model）は、ヘックマンモデルの特殊ケースです。結果方程式と選択方程式が同一の潜在変数で決定される場合に対応します。
$$ Y_i = \begin{cases} Y_i^* = \bm{X}_i^\top \bm{\beta} + \epsilon_i & \text{if } Y_i^* > 0 \\ 0 & \text{if } Y_i^* \leq 0 \end{cases} $$
トービットモデルでは、$S_i^* = Y_i^*$（選択と結果が同一の潜在変数）であるため、$\rho = 1$ です。これはヘックマンモデルの「選択方程式 = 結果方程式」という極端なケースです。
例えば、家計の耐久消費財への支出額の分析で、支出がゼロの家計が多い場合、「支出するかどうか」と「いくら支出するか」が同一の潜在変数（消費意欲）で決まると仮定するのがトービットモデルです。
ヘックマンモデルとの違い
ヘックマンモデルはトービットモデルよりも一般的です。
| 特徴 | トービットモデル | ヘックマンモデル |
|---|---|---|
| 選択メカニズム | $Y^* > 0$ で観測 | $S^* > 0$ で観測 |
| $\rho$ | $\rho = 1$（固定） | $\rho$ は推定対象 |
| 柔軟性 | 低い | 高い |
| 除外変数 | 不要 | 望ましい |
選択バイアスの実例と応用
賃金関数の推定（ヘックマンの原論文）
ヘックマンの1979年の原論文は、既婚女性の労働供給と賃金関数を分析したものです。当時の米国では既婚女性の就業率は約50%であり、賃金データは就業者からしか得られませんでした。
結果方程式（賃金関数）は教育年数、経験年数、経験の二乗を説明変数とし、選択方程式（就業決定）は賃金関数の変数に加えて、配偶者の所得、子供の数、子供の年齢を除外変数として含めました。ヘックマン補正なしのOLSでは教育の収益率が過小推定される傾向が見られ、選択バイアスの実証的重要性を示す先駆的な研究となりました。
金融における信用リスク評価
銀行の融資審査において、デフォルト率の分析に選択バイアスが生じます。デフォルトの有無が観測されるのは審査に合格して融資を受けた企業のみであり、審査で不合格になった企業のデフォルト率は分かりません。審査基準が信用力（デフォルトリスクに関連）に基づいているため、融資を受けた企業は平均的に信用力が高く、デフォルト率が母集団全体よりも低くなります。
この場合、除外変数としては「過去の取引関係の有無」「紹介者の有無」など、融資審査には影響するがデフォルト確率には直接影響しない変数が候補となります。
臨床試験における脱落問題
臨床試験で患者が途中で脱落（dropout）する場合も選択バイアスの一種です。副作用が強い患者や効果が見られない患者ほど脱落しやすく、試験を完了した患者のみで分析すると、治療効果が過大評価される可能性があります。この問題に対しては、ヘックマン型のモデルの他に、パターン混合モデル（pattern mixture model）や共有パラメータモデル（shared parameter model）が用いられます。
これらの実例から、選択バイアスが分野を超えて広く存在する問題であることが分かります。次に、因果推論の枠組みから選択バイアスを再解釈します。
因果推論の枠組みからの再解釈
選択バイアスとDAG
DAGと因果グラフの観点から、選択バイアスはコライダーへの条件付け（conditioning on a collider）として理解できます。
選択変数 $S$ がアウトカム $Y$ と説明変数 $X$ の共通の結果（コライダー）である場合、$S = 1$ で条件付ける（選択されたサンプルのみを分析する）と、$X$ と $Y$ の間に擬似的な関連が生じます。
$$ X \rightarrow Y, \quad X \rightarrow S, \quad \epsilon_Y \rightarrow Y \rightarrow S $$
$S$ に条件付けると、$X$ と $\epsilon_Y$ の間にバックドアパスが開き、$X$ の係数にバイアスが生じます。
ヘックマン補正の因果的解釈
ヘックマン補正は、この「コライダーへの条件付け」によるバイアスを逆ミルズ比で制御しています。逆ミルズ比は選択プロセスの情報を要約しており、これを回帰に含めることで、コライダーバイアスを除去します。
この解釈は、ヘックマン補正が単なる統計的テクニックではなく、因果推論の原理に根ざした手法であることを示しています。
ここまでの理論を、Pythonシミュレーションで確認しましょう。
Pythonシミュレーション
シミュレーション1: ヘックマン2段階推定法の実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.linear_model import LinearRegression
np.random.seed(42)
n = 5000
# --- データ生成 ---
X = np.random.normal(0, 1, n) # 結果と選択に共通の変数
Z_excl = np.random.normal(0, 1, n) # 除外変数（選択のみに影響）
# 誤差項（相関あり: rho = 0.6）
rho = 0.6
sigma = 2.0
cov_matrix = [[sigma**2, rho * sigma], [rho * sigma, 1]]
errors = np.random.multivariate_normal([0, 0], cov_matrix, n)
epsilon = errors[:, 0]
nu = errors[:, 1]
# 結果方程式（真のパラメータ）
beta_0, beta_1 = 3.0, 2.0
Y_star = beta_0 + beta_1 * X + epsilon
# 選択方程式
gamma_0, gamma_1, gamma_2 = -0.5, 0.3, 0.8
S_star = gamma_0 + gamma_1 * X + gamma_2 * Z_excl + nu
S = (S_star > 0).astype(int)
print(f"選択率: {S.mean():.3f}")
print(f"全サンプルのY平均: {Y_star.mean():.4f}")
print(f"選択されたサンプルのY平均: {Y_star[S==1].mean():.4f}")
print(f"選択バイアス（平均値の差）: {Y_star[S==1].mean() - Y_star.mean():.4f}")
# --- 推定 ---
# (a) ナイーブOLS（選択サンプルのみ）
X_sel = X[S == 1]
Y_sel = Y_star[S == 1]
ols_naive = LinearRegression().fit(X_sel.reshape(-1, 1), Y_sel)
beta_naive = [ols_naive.intercept_, ols_naive.coef_[0]]
# (b) ヘックマン2段階推定
# 第1段階: プロビットの近似（scipy.optimizeで正確なプロビット）
from scipy.optimize import minimize
def neg_log_lik_probit(params, Z, S):
pred = Z @ params
ll = np.sum(S * np.log(norm.cdf(pred) + 1e-10)
+ (1 - S) * np.log(1 - norm.cdf(pred) + 1e-10))
return -ll
Z_full = np.column_stack([np.ones(n), X, Z_excl])
gamma_init = np.zeros(3)
result = minimize(neg_log_lik_probit, gamma_init, args=(Z_full, S),
method='BFGS')
gamma_hat = result.x
# 逆ミルズ比
probit_pred = Z_full @ gamma_hat
imr = norm.pdf(probit_pred) / norm.cdf(probit_pred)
# 第2段階: 逆ミルズ比を追加してOLS
X_heckman = np.column_stack([X_sel, imr[S == 1]])
ols_heckman = LinearRegression().fit(X_heckman, Y_sel)
beta_heckman = [ols_heckman.intercept_, ols_heckman.coef_[0]]
delta_hat = ols_heckman.coef_[1]
# (c) 全サンプルOLS（オラクル）
ols_oracle = LinearRegression().fit(X.reshape(-1, 1), Y_star)
beta_oracle = [ols_oracle.intercept_, ols_oracle.coef_[0]]
print(f"\n=== 推定結果 ===")
print(f"真の値: beta_0={beta_0:.4f}, beta_1={beta_1:.4f}")
print(f"オラクル: beta_0={beta_oracle[0]:.4f}, beta_1={beta_oracle[1]:.4f}")
print(f"ナイーブ: beta_0={beta_naive[0]:.4f}, beta_1={beta_naive[1]:.4f}")
print(f"ヘックマン: beta_0={beta_heckman[0]:.4f}, beta_1={beta_heckman[1]:.4f}")
print(f"delta(=rho*sigma): {delta_hat:.4f} (真値: {rho*sigma:.4f})")
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
ax = axes[0]
ax.scatter(X[S==0], Y_star[S==0], alpha=0.1, s=5, color="gray", label="Not selected")
ax.scatter(X[S==1], Y_star[S==1], alpha=0.15, s=5, color="blue", label="Selected")
x_plot = np.linspace(-3, 3, 100)
ax.plot(x_plot, beta_0 + beta_1 * x_plot, "g-", linewidth=2.5, label="True")
ax.plot(x_plot, beta_naive[0] + beta_naive[1] * x_plot, "r--", linewidth=2,
label="Naive OLS")
ax.plot(x_plot, beta_heckman[0] + beta_heckman[1] * x_plot, "b:", linewidth=2,
label="Heckman")
ax.set_xlabel("X", fontsize=12)
ax.set_ylabel("Y", fontsize=12)
ax.set_title("Selection Bias and Heckman Correction", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax = axes[1]
methods = ["True", "Oracle", "Naive", "Heckman"]
b0_vals = [beta_0, beta_oracle[0], beta_naive[0], beta_heckman[0]]
b1_vals = [beta_1, beta_oracle[1], beta_naive[1], beta_heckman[1]]
x_pos = np.arange(len(methods))
w = 0.35
ax.bar(x_pos - w/2, b0_vals, w, label=r"$\beta_0$", color="steelblue", alpha=0.7)
ax.bar(x_pos + w/2, b1_vals, w, label=r"$\beta_1$", color="salmon", alpha=0.7)
ax.set_xticks(x_pos)
ax.set_xticklabels(methods, fontsize=10)
ax.set_ylabel("Coefficient", fontsize=12)
ax.set_title("Coefficient Comparison", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("heckman_correction.png", dpi=150, bbox_inches="tight")
plt.show()
左のグラフでは、灰色の点が選択されなかった個体、青の点が選択された個体です。選択された個体のアウトカムが系統的に高い方に偏っていることがわかります。緑の実線が真の回帰線、赤の破線がナイーブOLS、青の点線がヘックマン補正です。
重要な観察は以下の通りです。
- ナイーブOLSの切片は真の値より大きい: 選択バイアスにより、観測サンプルのアウトカムが系統的に高くなるため、切片が上方にシフトしています。
- ヘックマン補正は真の回帰線に近い: 逆ミルズ比の追加により選択バイアスが補正され、切片と傾きの両方が改善しています。
- $\hat{\delta} = \hat{\rho}\hat{\sigma}$ は正の値: 正の選択バイアス（高いアウトカムを持つ個体ほど選択されやすい）の存在を反映しています。
シミュレーション2: 除外変数の重要性
除外変数がある場合とない場合で、ヘックマン推定の安定性がどう変わるかを確認します。
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression
np.random.seed(0)
n_sims = 500
n = 1000
rho = 0.5
sigma = 1.5
beta_true = [2.0, 1.5]
results_with_excl = {'b0': [], 'b1': []}
results_without_excl = {'b0': [], 'b1': []}
for _ in range(n_sims):
X = np.random.normal(0, 1, n)
Z_excl = np.random.normal(0, 1, n)
errors = np.random.multivariate_normal([0, 0],
[[sigma**2, rho*sigma], [rho*sigma, 1]], n)
epsilon, nu = errors[:, 0], errors[:, 1]
Y_star = beta_true[0] + beta_true[1] * X + epsilon
S_star = -0.3 + 0.4 * X + 0.7 * Z_excl + nu
S = (S_star > 0).astype(int)
if S.sum() < 50 or (1-S).sum() < 50:
continue
X_sel = X[S == 1]
Y_sel = Y_star[S == 1]
# --- 除外変数あり ---
def probit_nll(params, Z, S):
p = norm.cdf(Z @ params)
return -np.sum(S*np.log(p+1e-10) + (1-S)*np.log(1-p+1e-10))
Z_with = np.column_stack([np.ones(n), X, Z_excl])
res = minimize(probit_nll, np.zeros(3), args=(Z_with, S), method='BFGS')
imr_with = norm.pdf(Z_with @ res.x) / norm.cdf(Z_with @ res.x)
reg = LinearRegression().fit(
np.column_stack([X_sel, imr_with[S==1]]), Y_sel)
results_with_excl['b0'].append(reg.intercept_)
results_with_excl['b1'].append(reg.coef_[0])
# --- 除外変数なし ---
Z_without = np.column_stack([np.ones(n), X])
res2 = minimize(probit_nll, np.zeros(2), args=(Z_without, S), method='BFGS')
imr_without = norm.pdf(Z_without @ res2.x) / norm.cdf(Z_without @ res2.x)
reg2 = LinearRegression().fit(
np.column_stack([X_sel, imr_without[S==1]]), Y_sel)
results_without_excl['b0'].append(reg2.intercept_)
results_without_excl['b1'].append(reg2.coef_[0])
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
for idx, (key, true_val, label) in enumerate(
[('b0', beta_true[0], r'$\beta_0$'),
('b1', beta_true[1], r'$\beta_1$')]):
ax = axes[idx]
ax.hist(results_with_excl[key], bins=40, alpha=0.6, density=True,
color='steelblue', label='With exclusion restriction')
ax.hist(results_without_excl[key], bins=40, alpha=0.6, density=True,
color='salmon', label='Without exclusion restriction')
ax.axvline(true_val, color='green', linewidth=2, linestyle='--',
label=f'True = {true_val}')
ax.set_xlabel(f'{label} estimate', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'Distribution of {label} estimates', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('heckman_exclusion.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"=== beta_0の推定 ===")
print(f"除外変数あり: mean={np.mean(results_with_excl['b0']):.4f}, "
f"SD={np.std(results_with_excl['b0']):.4f}")
print(f"除外変数なし: mean={np.mean(results_without_excl['b0']):.4f}, "
f"SD={np.std(results_without_excl['b0']):.4f}")
print(f"\n=== beta_1の推定 ===")
print(f"除外変数あり: mean={np.mean(results_with_excl['b1']):.4f}, "
f"SD={np.std(results_with_excl['b1']):.4f}")
print(f"除外変数なし: mean={np.mean(results_without_excl['b1']):.4f}, "
f"SD={np.std(results_without_excl['b1']):.4f}")
このシミュレーション結果から、除外変数の有無がヘックマン推定の安定性に与える影響が明確にわかります。
- 除外変数ありの場合（青）: 推定値の分布は真の値を中心に集中しており、標準偏差も比較的小さいです。推定が安定しています。
- 除外変数なしの場合（赤）: 推定値の分布が大きく広がっており、標準偏差が数倍に増加しています。個々のシミュレーションで推定値が大きく変動し、不安定です。
- 両方とも平均は真の値に近い: 除外変数がなくても平均的にはバイアスは小さいですが、分散が非常に大きいため、個々の推定は信頼できません。
この結果は、除外変数が「識別の点では不要だが、推定の安定性のために極めて重要」であることを実証しています。
まとめ
本記事では、選択バイアスとヘックマン補正について以下を解説しました。
- 選択バイアスは、サンプルに含まれるかどうかがアウトカムと相関する変数に依存する場合に生じる系統的な偏り
- サンプル選択モデルは結果方程式と選択方程式の2つの潜在変数モデルで定式化される。誤差項の相関 $\rho$ が選択バイアスの源
- 逆ミルズ比 $\lambda(a) = \phi(a)/\Phi(a)$ は選択条件付き期待値を閉じた形で表現する関数
- ヘックマン2段階推定法: 第1段階でプロビットモデルを推定し逆ミルズ比を計算、第2段階で逆ミルズ比を追加変数としてOLS推定
- 除外変数の存在が推定の安定性に決定的に重要。除外変数がなくても技術的には識別可能だが、推定が不安定になる
- FIMLは2段階推定法より効率的だが、数値最適化の収束に注意が必要
- DAGの観点から、選択バイアスはコライダーへの条件付けとして理解できる
ヘックマン補正は、選択バイアスに対する理論的にエレガントかつ実用的な解法です。しかし、二変量正規性の仮定と除外変数の選択という2つの重要な判断が必要であり、これらの妥当性を十分に検討することが不可欠です。
次のステップとして、以下の記事も参考にしてください。
- 操作変数法 — 内生性への別のアプローチ
- 因果推論の感度分析 — 仮定のロバスト性の評価
- 交絡バイアスとその制御法 — バイアスの基本的な制御方法