Rの関数:iv_robust {estimatr}

Rの関数から iv_robust {estimatr} を確認します。

本ポストはこちらの続きです。

Rの関数:horvitz_thompson {estimatr}
Rの関数から horvitz_thompson {estimatr} を確認します。本ポストはこちらの続きです。関数 horvitz_thompson とはhorvitz_thompson は、Horvitz-Thompson 推定量（IP...
www.saecanet.com
2025.12.14

iv_robust とは

iv_robust は、操作変数法（Instrumental Variable Method: IV法） を用いて、回帰モデルのパラメータを推定する関数です。具体的には、2段階最小二乗法（2SLS）を実行します。

この関数の特徴は、関数名が示す通り、「頑健性（Robustness）」 にあります。

通常、操作変数法を行うための関数（ivreg など）では、標準誤差の計算において「誤差項の等分散性」や「観測値の独立性」を仮定することが一般的です。

しかし、現実のデータでは不均一分散やクラスター相関が存在することが多く、標準的な検定は信頼できません。

iv_robust は、内生性（Endogeneity）に対処するための2SLS推定量と、不均一分散やクラスター相関に対処するための頑健標準誤差（Robust Standard Errors）を同時に、かつシームレスに提供します。

iv_robust の引数

library(estimatr)
args(iv_robust)
function (formula, data, weights, subset, clusters, fixed_effects, 
    se_type = NULL, ci = TRUE, alpha = 0.05, diagnostics = FALSE, 
    return_vcov = TRUE, try_cholesky = FALSE) 
NULL
  1. formula

    • 解析モデルを指定する式です。AER::ivreg などと同様の拡張公式構文を使用します。
    • 形式: outcome ~ controls + treatment | controls + instrument

      • ~ の左辺：被説明変数（Y）。
      • ~ の右辺（| の左側）：内生変数（X）と外生変数（W）。
      • | の右側：操作変数（Z）と外生変数（W）。
    • ポイント: 外生変数は | の両側に記述する必要があります。記述された変数のうち、左側にしか存在しない変数が「内生変数」、右側にしか存在しない変数が「操作変数」として自動的に認識されます。
  2. data

    • 解析対象のデータフレームです。
  3. weights

    • 重み付き最小二乗法を行う場合の重みベクトルです。
  4. subset

    • 分析に使用するデータのサブセット条件を指定します。
  5. clusters

    • クラスター変数を指定します。
    • データがグループ構造を持つ場合（例：学校ごとの生徒データ）、クラスターロバスト標準誤差を計算するために使用します。
  6. fixed_effects

    • 固定効果（Fixed Effects）として扱う変数を指定します。
    • グループごとの切片の違いを制御します。
  7. se_type

    • 標準誤差の計算方法を指定する文字列です。
    • デフォルト: NULL（クラスター指定なしなら "HC2"、ありなら "CR2" が自動選択されます）。
    • 選択肢: "HC1", "HC2", "HC3", "stata", "CR0", "CR2", "classical" など。
  8. ci, alpha

    • 信頼区間の計算有無（TRUE/FALSE）と有意水準（デフォルト 0.05）です。
  9. diagnostics

    • 操作変数法の診断テスト結果を出力するかどうかの論理値です。
    • デフォルト: FALSE
    • TRUE にすると、弱操作変数検定（第1段階のF値）、Wu-Hausman検定（内生性の検定）、過剰識別検定（Sargan検定など）の結果が含まれます。
  10. return_vcov

    • 分散共分散行列を返すかどうかの論理値です。

シミュレーションコード

以下に、iv_robust の機能を確認するためのサンプルデータを用いたシミュレーションコードを示します。

シナリオ: 職業訓練プログラムが賃金に与える効果

ある職業訓練プログラム（Training）が、後の賃金（Wage）に与える因果効果を知りたいとします。

  • 課題（内生性）:

    • 「個人の能力（Ability）」や「やる気」といった観察できない要因（交絡因子）が存在します。能力が高い人は、訓練を受けやすく、かつ将来の賃金も高い傾向があります。このため、単に訓練を受けた人と受けていない人の賃金を比較（OLS）すると、訓練の効果が過大評価（上方バイアス）されてしまいます。
  • 解決策（操作変数）:

    • 「訓練会場への距離（Distance）」を操作変数として利用します。
    • 関連性: 距離が近いほど訓練に参加しやすい。
    • 除外制約: 距離そのものは将来の賃金や個人の能力には直接影響しないと仮定します。

この設定において、不均一分散も導入し、iv_robust がバイアスを補正しつつ適切な標準誤差を提供することを確認します。

なお、有意水準は5%とします。

シミュレーションデータの生成

  • 内生性の源泉: Ability が Training と Wage の両方に正の影響を与える。
  • 操作変数: Distance_Near は Training に影響するが Wage には直接影響しない。
# パッケージの読み込み
library(dplyr)
library(texreg) # 結果を見やすく表示するため

# 乱数シードの固定
seed <- 20251216
set.seed(seed)

# サンプルサイズ
N <- 1000

# 1. 未観測の交絡因子 (U: Ability)
# 能力が高いほど値が大きい（平均0、分散1）
U_ability <- rnorm(N, mean = 0, sd = 1)

# 2. 操作変数 (Z: Distance)
# 訓練会場への距離（近い=1, 遠い=0 のランダム割り当てを想定）
# 実際には連続変数でも良いが、ここでは分かりやすく2値とする
Z_distance_near <- rbinom(N, size = 1, prob = 0.5)

# 3. 内生変数 (X: Training)
# 訓練に参加するかどうか (0 or 1)
# 参加確率は以下に依存:
#  - 距離が近いほど参加しやすい (Zの効果: 正)
#  - 能力が高いほど参加しやすい (Uの効果: 正) <- これが内生性の原因
prop_training <- pnorm(-1 + 1.5 * Z_distance_near + 0.8 * U_ability)
X_training <- rbinom(N, size = 1, prob = prop_training)

# 4. 結果変数 (Y: Wage)
# 賃金は以下に依存:
#  - 訓練効果 (真の効果: +1000ドル)
#  - 能力 (Uの効果: 正) <- これが交絡
#  - ノイズ (不均一分散: 訓練参加者は賃金のばらつきが大きいとする)
true_causal_effect <- 1000
noise_sd <- ifelse(X_training == 1, 800, 400) # 不均一分散
epsilon <- rnorm(N, mean = 0, sd = noise_sd)

Y_wage <- 2000 + (true_causal_effect * X_training) + (500 * U_ability) + epsilon

# データフレームの作成
df_sim <- data.frame(
  Wage = Y_wage,
  Training = X_training,
  Distance_Near = Z_distance_near,
  Ability = U_ability # 実際には分析者は観測できない
)

cat("--- データ概要 ---\n")
cat(sprintf("総観測数: %d\n", N))
cat(sprintf("真の訓練効果: +%d\n", true_causal_effect))
--- データ概要 ---
総観測数: 1000
真の訓練効果: +1000

分析1: 通常の最小二乗法 (OLS)

# 内生性を無視して推定
# robust標準誤差は使用 (HC2)
ols_model <- lm_robust(Wage ~ Training, data = df_sim, se_type = "HC2")

cat("--- OLSの結果（Abilityの影響によるバイアスが含まれる） ---\n")
print(summary(ols_model))
--- OLSの結果（Abilityの影響によるバイアスが含まれる） ---

Call:
lm_robust(formula = Wage ~ Training, data = df_sim, se_type = "HC2")

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value   Pr(>|t|) CI Lower CI Upper  DF
(Intercept)     1814      26.30   68.96  0.000e+00     1762     1866 998
Training        1425      52.31   27.24 1.237e-122     1323     1528 998

Multiple R-squared:  0.4485 ,   Adjusted R-squared:  0.448 
F-statistic: 742.2 on 1 and 998 DF,  p-value: < 2.2e-16

分析2: 操作変数法 (iv_robust)

# iv_robust を用いて 2SLS 推定を行う
# モデル式: Y ~ X | Z (内生変数 Training を 操作変数 Distance_Near で予測)
# diagnostics = TRUE で診断テストも実施

iv_model <- iv_robust(
  formula = Wage ~ Training | Distance_Near,
  data = df_sim,
  se_type = "HC2", # 不均一分散頑健標準誤差
  diagnostics = TRUE
)

cat("--- IV推定の結果（バイアスが補正される） ---\n")
print(summary(iv_model))
--- IV推定の結果（バイアスが補正される） ---

Call:
iv_robust(formula = Wage ~ Training | Distance_Near, data = df_sim, 
    se_type = "HC2", diagnostics = TRUE)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value   Pr(>|t|) CI Lower CI Upper  DF
(Intercept)   2041.3       57.4  35.562 1.383e-179   1928.6     2154 998
Training       914.3      128.2   7.133  1.885e-12    662.8     1166 998

Multiple R-squared:  0.3909 ,   Adjusted R-squared:  0.3903 
F-statistic: 50.88 on 1 and 998 DF,  p-value: 1.885e-12

Diagnostics:
                 numdf dendf  value  p.value    
Weak instruments     1   998 200.76  < 2e-16 ***
Wu-Hausman           1   997  23.36 1.55e-06 ***
Overidentifying      0    NA     NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

結果の比較

models_list <- list("OLS (Biased)" = ols_model, "IV (Robust)" = iv_model)
screenreg(models_list, include.ci = FALSE, digits = 2)

======================================
             OLS (Biased)  IV (Robust)
--------------------------------------
(Intercept)  1813.94 ***   2041.27 ***
              (26.30)       (57.40)   
Training     1425.21 ***    914.35 ***
              (52.31)      (128.19)   
--------------------------------------
R^2             0.45          0.39    
Adj. R^2        0.45          0.39    
Num. obs.    1000          1000       
RMSE          786.12        826.17    
======================================
*** p < 0.001; ** p < 0.01; * p < 0.05
OLSモデル

まず、lm_robust（OLS）の結果に注目します。

  • Trainingの推定値: 1425 （標準誤差: 52.31）
  • 真の値: 1000

推定値は真の値よりも 約425 大きく算出されました。

シミュレーションの設定上、「能力（Ability）」が高い人ほど訓練を受けやすく、かつ賃金も高い傾向にあります。

OLSはこの「能力による賃金の高さ」と「訓練による賃金上昇」を区別できません。

それゆえ、能力によるプラス効果が訓練効果に誤って加算され、結果として過大評価（上方バイアス）が生じています。

また、標準誤差が、下記のIVモデルと比較して、小さく（52.31）、結果、p値が設定した有意水準を下回っているため、ミスリーディングな結果を招いています。

IVモデル

次に、iv_robust（操作変数法）の結果を確認します。

  • Trainingの推定値: 914.3 （標準誤差: 128.2）
  • 真の値: 1000

こちらの推定値は、OLSの結果から低下し、真の値である 1000 に近い水準まで補正されました。

95%信頼区間 [662.8, 1166] は真の値 1000 を含んでいます。

操作変数である「会場への距離」を利用することで、能力（Ability）とは無関係な「訓練受講の変動」のみを抽出し、純粋な因果効果を取り出すことに成功しています。

一方で、標準誤差は 128.2 となり、OLSの約2.5倍に拡大しました。

操作変数を用いることは、データに含まれる情報の一部（内生的な変動）を捨て、クリーンな情報（外生的な変動）のみを使って推定を行うことを意味します。

有効な情報量が減るため、推定の「不確実性（標準誤差）」は必然的に大きくなります。

診断テスト（Diagnostics）の評価
  • Weak instruments（弱操作変数検定）:

    • F値: 200.76
    • 解釈: 一般的な目安である「10」を大きく上回っています。操作変数（距離）が内生変数（訓練）を効果的にに説明していることを示しており、分析の前提条件が満たされていることが確認できます。
  • Wu-Hausman（内生性検定）:

    • P値: 1.55e-06 （設定した有意水準より小さい）
    • 解釈: 「OLSとIVの推定値に統計的に有意な差はない（＝内生性はない）」という帰無仮説が棄却されました。

以上です。