Rの関数:ivreg {AER}

1. ivreg関数の目的と概要

通常の回帰分析（lm()関数など）では、説明変数が誤差項と相関していない（外生性を満たす）という仮定が置かれています。しかし、以下のような場合、この仮定が崩れ、推定結果にバイアス（偏り）が生じます。これを内生性 (endogeneity) の問題と呼びます。

• 欠落変数バイアス (Omitted Variable Bias): モデルに含まれていない変数が、被説明変数と説明変数の両方に影響を与えている。
  • 例：「教育年数」が「賃金」に与える影響を分析したいが、個人の「能力」という観測できない変数が、教育年数（能力が高い人は高学歴になりやすい）と賃金（能力が高い人は高賃金）の両方に影響している。
• 同時決定性 (Simultaneity): 被説明変数と説明変数が相互に影響を与え合っている。
  • 例：「価格」が「需要」に与える影響を見たいが、需要が高いと価格も上がるという逆の因果関係も存在する。
• 測定誤差 (Measurement Error): 説明変数に測定誤差が含まれている。

このような内生性の問題がある場合、通常の最小二乗法（OLS）は正しい因果効果を推定できません。そこで、操作変数法 (IV) を用います。

操作変数（Instrumental Variable, IV）とは、以下の2つの条件を満たす変数のことです。

1. 関連性 (Relevance): 操作変数は、問題となっている内生変数と相関がある。
2. 除外制約 (Exclusion Restriction): 操作変数は、被説明変数に対しては、内生変数を通じてのみ影響を与え、直接的な影響を持たない（誤差項と相関しない）。

ivreg関数は、この操作変数を用いて、内生性の問題を解消し、バイアスのない推定を行うためのツールです。

2. ivreg関数の使い方：引数の解説

ivregの最も重要な引数は formula, data です。特にformulaの書き方が独特です。

ivreg(formula, data, ...)

3. 具体的な実行例

シナリオ設定: 生徒の観測できない「やる気」が、「追加補習時間」（説明変数）と「テストの点数」（被説明変数）の両方に強い正の影響を与えている状況を考えます。この「やる気」が原因で、「追加補習時間」は内生変数となり、OLSで分析すると補習の効果を過大評価してしまいます。塾から遠い生徒ほど補習に参加しにくいという事実を利用し、「通塾距離」を操作変数とします。

# パッケージの読み込み
library(AER)
library(lmtest)

# 再現性のための乱数シード固定
seed <- 20250828
set.seed(888)

# データ数
n <- 1000

# 1. 観測できない「やる気」を生成 (内生性の根源)
motivation <- rnorm(n, mean = 5, sd = 2)

# 2. 操作変数「通塾距離」と外生変数「基礎学力」を生成
distance <- runif(n, 1, 10)
baseline_ability <- rnorm(n, mean = 50, sd = 10)

# 3. 内生変数「追加補習時間」を生成
# 「やる気」が高いほど補習時間が増える
tutoring_hours <- 1 + 1.0 * motivation - 0.5 * distance + rnorm(n, 0, 1)

# 4. 被説明変数「テストの点数」を生成
# 誤差項に含まれる「やる気」の影響は「5」とする
# 補習時間の真の効果は「8」とする
error <- 5 * motivation + rnorm(n, 0, 5)
test_score <- 20 + 8 * tutoring_hours + 1 * baseline_ability + error

# データフレームにまとめる
sim_data <- data.frame(test_score, tutoring_hours, baseline_ability, distance)

ols_model <- lm(test_score ~ tutoring_hours + baseline_ability, data = sim_data)

cat("--- OLSモデルの結果（バイアスあり）---\n")
print(coeftest(ols_model, vcov = vcovHC, type = "HC1"))

--- OLSモデルの結果（バイアスあり）---

t test of coefficients:

                  Estimate Std. Error t value  Pr(>|t|)    
(Intercept)      33.533581   1.358714   24.68 < 2.2e-16 ***
tutoring_hours   11.033085   0.103325  106.78 < 2.2e-16 ***
baseline_ability  1.033718   0.026229   39.41 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

iv_model <- ivreg(test_score ~ tutoring_hours + baseline_ability | distance + baseline_ability,
  data = sim_data
)

cat("\n--- IVモデルの結果（バイアス修正後）---\n")
print(summary(iv_model, vcov = vcovHC, type = "HC1", diagnostics = TRUE))

--- IVモデルの結果（バイアス修正後）---

Call:
ivreg(formula = test_score ~ tutoring_hours + baseline_ability | 
    distance + baseline_ability, data = sim_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-32.8832  -7.4412   0.3245   7.9295  36.2679 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       44.0052     2.1233   20.73   <2e-16 ***
tutoring_hours     7.9408     0.2994   26.53   <2e-16 ***
baseline_ability   1.0290     0.0369   27.88   <2e-16 ***

Diagnostic tests:
                 df1 df2 statistic p-value    
Weak instruments   1 997     296.0  <2e-16 ***
Wu-Hausman         1 996     369.1  <2e-16 ***
Sargan             0  NA        NA      NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.23 on 997 degrees of freedom
Multiple R-Squared: 0.8618, Adjusted R-squared: 0.8615 
Wald test: 712.5 on 2 and 997 DF,  p-value: < 2.2e-16 

4. 結果の解釈と比較

• OLSの結果: tutoring_hoursの係数は 11.033085 と推定されました。これは、真の効果である「8」よりも、欠落変数バイアスの影響により、過大な評価となっています。

• IV (ivreg)の結果: 一方、ivregによるtutoring_hoursの係数は 7.9408 と推定されました。ivregはバイアスをほぼ除去し、真の値である「8」を捉えています。

• 診断テスト:

  • Weak instruments: F値が296.0であるため、操作変数の「関連性」は無いと考えられます。
  • Wu-Hausman: p値が<2.2e-16であることから、「内生性は存在しない」という帰無仮説を有意に棄却できます。

5. Diagnostic testsの解説