pacman::p_load(tidyverse,
modelsummary,
skimr,
ivreg,
estimatr)Rで学ぶ計量経済学と機械学習 6
計量経済学２：操作変数法
First version: 2019-05-24
This version: 2026-11-03
1 本講義の目的
- 操作変数法の実施方法を学ぶ。
2 パッケージの読み込み
今回の分析で使うパッケージを読み込む。
3 データの用意
前回と同じく、架空のデータを生成して分析を行う。
前回は個人の所得\(Y\)と学歴\(X\)・能力\(A\)との関係についての架空データを次のように生成した（ただし\(\varepsilon\)は誤差項）。
\[ Y = 200 + 10A + 500X + \varepsilon \]
- サンプルサイズは10,000人
- 能力は0から100まで均等に分布する
- 能力が1上がると所得は10万円上昇する
- 能力が80 以上の約 2,000人中から約 1,000人がランダムに選ばれて大卒となる。
- 大卒だと所得が500万円上昇する
- 切片は200万円（能力=0、非大卒の平均給与）
今回は、前回のデータに新たな変数を加える。「能力に関係なく、30％の確率でランダムに選ばれた人が学費を免除され、入学を許可される」という制度があるとし、さらに、「学費免除に選ばれた人は、50%の確率で大学を卒業できる」（ただし能力が80以上の条件で選ばれて大卒となった人はこの法則の対象外）とする。
データ生成過程は、以下のように設定できる。（難しく感じる人は、コピペでかまわない）
前回と同じく、二項分布の乱数を発生させるrbinom()とif_else()を用いて、mutate()の中で変数を1つずつ作っていく。前回のコードに、学費免除ダミーexemptionの作成と、それによる大卒ダミーuniversityの修正（2行）を加えればよい。
# 乱数の種を固定（同じ種からは常に同じ乱数が生成される。10101を他の数値に変えると異なる乱数となる。）
set.seed(10101)
# n:サンプルサイズ
n <- 10000
# rbinom(n, 1, p): 確率pで1、確率1-pで0となる乱数をn個生成
# if_else(条件, 条件を満たすときの値, 満たさないときの値)
df <- tibble(ID = 1:n,
ability = runif(n, min = 0, max = 100)) |> # 能力は0から100まで均等に分布
mutate(
# 能力が80以上の人は、50%の確率で大卒になる（前回と同じ設定）
university = if_else(ability >= 80, rbinom(n, 1, 0.5), 0),
# 能力に関係なく、30%の確率で学費免除になる
exemption = rbinom(n, 1, 0.3),
# まだ大卒でない学費免除者は、50%の確率で大卒になる
university = if_else(university == 0 & exemption == 1, rbinom(n, 1, 0.5), university),
# DGPに基づいて所得を生成（誤差項は平均=0、標準偏差=50）
income = 200 + 10*ability + 500*university + rnorm(n, mean = 0, sd = 50)
)
# 最初の6行
head(df)# A tibble: 6 × 5
ID ability university exemption income
<int> <dbl> <dbl> <int> <dbl>
1 1 19.0 0 0 434.
2 2 91.1 0 1 1089.
3 3 22.8 0 0 447.
4 4 82.5 1 0 1484.
5 5 91.6 0 0 1124.
6 6 50.5 0 1 684.生成したデータの記述統計を、第3回で学んだskimr::skim()で表示し、どのくらいが学費免除および大卒になったかチェックしてみる（0か1をとるダミー変数exemptionとuniversityの平均値meanは、それぞれ学費免除になった人と大卒の人の割合を表す）。
skimr::skim(df)| Name | df |
| Number of rows | 10000 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| numeric | 5 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| ID | 0 | 1 | 5000.50 | 2886.90 | 1.0 | 2500.75 | 5000.50 | 7500.25 | 10000.00 | ▇▇▇▇▇ |
| ability | 0 | 1 | 50.03 | 29.02 | 0.0 | 25.02 | 50.12 | 75.26 | 100.00 | ▇▇▇▇▇ |
| university | 0 | 1 | 0.24 | 0.42 | 0.0 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| exemption | 0 | 1 | 0.30 | 0.46 | 0.0 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▃ |
| income | 0 | 1 | 818.29 | 406.51 | 110.7 | 494.46 | 770.64 | 1041.76 | 1830.47 | ▆▇▇▂▃ |
こうして生成したデータの所得と能力・学歴の関係をプロットすると次のようになる。前回と異なり、abilityが80を下回っても、学費免除によって大学を卒業し、所得が上昇している人々がいることが分かる。
# 大卒か否かのラベルをデータフレームに加える
df <- df |> mutate(edu_label =
case_when(university == 1 ~ "Grad.",
university == 0 ~ "Not grad."))
# 散布図を描く with 大卒ラベル
ggplot(df, aes(x = ability, y = income, color = edu_label)) +
geom_point(alpha = 0.5) +
labs(title = "Ability and Income")
以下では、「学歴（大卒になること）が所得をどれだけ上昇させるのか」を推定することを目的として分析していく。
4 操作変数法の概要
効果を測りたい処置変数（学歴\(X\)）を通じてのみアウトカム変数（所得\(Y\))に影響を与える変数のことを操作変数（instrumental variable）という。そのような変数が存在することは多くないが、今回の場合は学費免除ダミー\(Z\)を操作変数として用いることができる。
前回、回帰分析などで学歴の効果を推定しようとした際には、交絡要因である能力がデータとして得られないという想定の下では 学歴の効果が過大に推定される問題が起こることを述べた。操作変数法（instrumental variable method、IV method）は操作変数から得られる情報を活用することで、そのような状況においてもよりよい推定を行う方法である。
4.1 操作変数法の求め方
ベーシックな操作変数推定値は、以下の２種類のやり方によって求めることができる（両者の推定値は一致する）。
- 二段階最小二乗法(TSLS)：学歴\(X\)の変動のうち、学費免除\(Z\)の変動による変動分\(\hat{X}\)を推計し（第一段階）、\(\hat{X}\)→所得\(Y\)の因果効果を推定する（第二段階）。
- 誘導型推定値を第一段階推定値で割る：学費免除→学歴の因果効果（第一段階）、学費免除→所得の因果効果（誘導型）はそれぞれ回帰分析で推定できるので、この二つの推定値を利用して、以下のように推定する
\[ 学歴X\to 所得Y \text{の因果効果} = \frac{学費免除Z \to 所得Y\text{の因果効果}}{学費免除Z \to 学歴X\text{の因果効果}} \]
実際には、第一段階の推定結果を開示しつつも、操作変数推定値およびその標準誤差は二段階最小二乗法(TSLS)によって求めることがほとんどである。
なお二段階最小二乗法を用いる場合は、手動で二段階の推定を行うのではなく、下記で提示するパッケージなどを使うこと。自ら\(\hat{X}\)を推定し、それを\(Y\)に回帰して推定値と標準誤差を得ると、誤った標準誤差を算出することになる。
5 Rで実践
5.1 欠落変数を無視した回帰分析
まず、欠落変数（能力\(A\))を無視して、\(Y\)を\(X\)に回帰してみる。単純な単回帰に加えて、学費免除\(Z\)をコントロール変数に加えた重回帰分析も行ってみる。
回帰分析結果は以下のようになる。
reg_biased1 <- lm(formula = income ~ university,
data = df)
reg_biased2 <- lm(formula = income ~ university + exemption,
data = df)
modelsummary(list(reg_biased1, reg_biased2),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 656.331*** | 675.631*** |
| (3.240) | (3.426) | |
| university | 686.845*** | 745.569*** |
| (6.673) | (7.567) | |
| exemption | -110.899*** | |
| (7.016) | ||
| Num.Obs. | 10000 | 10000 |
| R2 | 0.514 | 0.526 |
| R2 Adj. | 0.514 | 0.526 |
| F | 10594.263 | 5553.864 |
「大卒になることによる所得上昇の効果」は、約687（万円）および約746（万円）と過大推定される。
5.2 誘導型推定値を第一段階推定値で割る
二段階最小二乗法を実施する前に、操作変数法をより深く理解するために、第一段階推定と誘導型推定を単純にOLSで実行し、両者の推定値より操作変数推定値を求めてみよう。
それぞれ、以下の回帰式で推定できる。
- 第一段階（first stage）：\(X_i =\pi_0 + \pi_1 Z_i + v_i\)
- 誘導型（reduced form）：\(Y_i = \gamma_0 + \gamma_1 Z_i + w_i\)
求めたい因果効果は
\[ 学歴X\to 所得Y \text{の因果効果} = \frac{\gamma_1}{\pi_1} = \frac{\text{誘導型の}Z_i \text{の係数}}{\text{第一段階の}Z_i\text{の係数}} \]
となる。
第一段階の推定結果は以下のようになる。
first_stage <- lm(formula = university ~ exemption,
data = df)
modelsummary(first_stage,
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| (Intercept) | 0.100*** |
| (0.004) | |
| exemption | 0.455*** |
| (0.008) | |
| Num.Obs. | 10000 |
| R2 | 0.241 |
| R2 Adj. | 0.241 |
| F | 3176.517 |
学費免除ダミーの係数推定値は、約0.46である。本モデルは線形確率モデルであるため、この結果は、学費免除になると大卒になる確率が約46%ポイント上昇することを示している。
なお、データ生成過程では「学費免除に選ばれた人は、50%の確率で大学を卒業できる」と設定したが、なぜここでは50%より小さい値が推定されるのだろうか。
それは、学費免除者（10000×0.3で約3,000人）の中には、能力が80以上で学費免除とは関係なく大卒となる人（3000×0.2×0.5で約300人）が含まれているからである。学費免除によって大卒になるのは、残りの約2,700人のうちの50%の約1,350人なので、学費免除の「因果効果」は、1350/3000=0.45となる（学費免除者の大卒率は(300+1350)/3000=0.55、学費免除でない人の大卒率は0.2×0.5=0.1であり、その差が0.45である）。第一段階推定値はこの値をほぼ正確に推定している。
一方、誘導型の推定結果は以下のようになる。
reduced_form <- lm(formula = income ~ exemption,
data = df)
modelsummary(reduced_form,
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| (Intercept) | 749.965*** |
| (4.692) | |
| exemption | 228.585*** |
| (8.581) | |
| Num.Obs. | 10000 |
| R2 | 0.066 |
| R2 Adj. | 0.066 |
| F | 709.559 |
学費免除ダミーの係数推定値は、約229であり、これは、「学費免除にあたると（大卒になることを通じて）平均所得が約229万円上昇する」ことを意味している。つまり誘導型推定値は、「学費免除の所得に対する平均処置効果」である。
最後に、両者の係数を割ると、操作変数推定値を得ることができる。
reduced_form$coefficients[2] / first_stage$coefficients[2]exemption
502.0146 「学歴（大卒）が所得に与える効果」である500万円を、ほぼ推定できている。ただし、このような方法では、標準誤差を得ることはできない。
5.3 二段階最小二乗法(TSLS)を実行する
{ivreg}パッケージに操作変数法による回帰を行うivreg()という関数がある。これを実行すれば操作変数法で推定し、標準誤差を出すことができる。
（注）ivreg()は、もともと{AER}パッケージに含まれていた関数である。現在は同じ開発者たちによって独立した{ivreg}パッケージとして開発が続けられており、本講義ではこちらを使う。書き方や推定結果はAER::ivreg()と同じであり、{ivreg}版では回帰診断（外れ値の影響の確認など）の機能が追加されている。古い教科書やウェブ上の資料ではAER::ivreg()がよく使われているが、そのまま読み替えてよい。なお、両方のパッケージを読み込むと関数名が重なるので、ivreg::ivreg()のようにパッケージ名を明示して使うとよい。
iv_reg <- ivreg::ivreg(
formula = income ~ university | exemption, # 被説明変数 ~ 説明変数 | 操作変数 と指定
data = df
)
modelsummary(iv_reg,
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| (1) | |
|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |
| (Intercept) | 699.914*** |
| (4.438) | |
| university | 502.015*** |
| (14.102) | |
| Num.Obs. | 10000 |
| R2 | 0.477 |
| R2 Adj. | 0.477 |
TSLS推定値は、約502であり、先に求めた値とも一致する。
5.4 弱操作変数の確認
summary()でdiagnostics = TRUEと指定すると、推定結果の下の「Diagnostic tests」に、操作変数法に関する検定の結果が表示される。
iv_summary <- summary(iv_reg, diagnostics = TRUE)
iv_summary
Call:
ivreg::ivreg(formula = income ~ university | exemption, data = df)
Residuals:
Min 1Q Median 3Q Max
-595.4118 -252.5841 0.6074 253.5689 628.5388
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 699.914 4.438 157.7 <2e-16 ***
university 502.015 14.102 35.6 <2e-16 ***
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 9998 3176.5 <2e-16 ***
Wu-Hausman 1 9997 249.8 <2e-16 ***
Sargan 0 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 293.9 on 9998 degrees of freedom
Multiple R-Squared: 0.4772, Adjusted R-squared: 0.4772
Wald test: 1267 on 1 and 9998 DF, p-value: < 2.2e-16 - Weak instruments：第一段階において、操作変数の係数が0であるという帰無仮説の\(F\)検定（操作変数が1つの場合は、第一段階の回帰の表にある\(F\)値と同じ値になる）。操作変数と処置変数の相関が弱い（弱操作変数）と、操作変数推定値に大きな偏りが生じたり、標準誤差が正しく計算できなくなったりする。目安として、\(F\)値が10を上回るかどうかがよく使われる（近年は、より厳しい基準も提案されている）。今回の\(F\)値は約3,177であり、十分大きい。
- Wu-Hausman：OLS推定値と操作変数推定値に差があるかの検定。帰無仮説は「処置変数は外生である（OLSでも偏りがない）」であり、棄却されれば、処置変数が内生的である（OLSには偏りがある）と考えられる。
- Sargan：操作変数の数が処置変数の数より多い場合に、操作変数の妥当性を調べる検定（過剰識別制約の検定）。今回は操作変数が1つなので計算されない（
NAとなる）。
5.5 ロバスト標準誤差
ただし、TSLSにおいても、一般的には、不均一分散に対して頑健な標準誤差（ロバスト標準誤差）を使うことが望ましい。そこで、{estimatr}パッケージのiv_robust()を使えば、ロバスト標準誤差を用いてTSLSを実施できる。実際にデータ分析する際には、この関数を活用することが望ましい。
iv_reg_robust <- estimatr::iv_robust(
income ~ university | exemption, # 被説明変数 ~ 説明変数 | 操作変数
data = df,
se_type = "HC1"
)
# summary
summary(iv_reg_robust)
Call:
estimatr::iv_robust(formula = income ~ university | exemption,
data = df, se_type = "HC1")
Standard error type: HC1
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 699.9 4.413 158.61 0.000e+00 691.3 708.6 9998
university 502.0 14.243 35.25 1.241e-256 474.1 529.9 9998
Multiple R-squared: 0.4772 , Adjusted R-squared: 0.4772
F-statistic: 1242 on 1 and 9998 DF, p-value: < 2.2e-16# 結果表の出力
modelsummary(list("ivreg（通常の標準誤差）" = iv_reg,
"iv_robust（ロバスト標準誤差）" = iv_reg_robust),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| ivreg（通常の標準誤差） | iv_robust（ロバスト標準誤差） | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 699.914*** | 699.914*** |
| (4.438) | (4.413) | |
| university | 502.015*** | 502.015*** |
| (14.102) | (14.243) | |
| Num.Obs. | 10000 | 10000 |
| R2 | 0.477 | 0.477 |
| R2 Adj. | 0.477 | 0.477 |
（参考）modelsummary()で標準誤差を書き換える
ivreg()で推定した結果を使う場合でも、表を作る段階でmodelsummary()のvcov引数を指定すれば、ロバスト標準誤差に書き換えることができる（計算には内部で{sandwich}パッケージが使われる）。vcov = "HC1"とすると、iv_robust()のse_type = "HC1"と同じロバスト標準誤差になる。
# 同じiv_regを2回並べ、1列目は通常の標準誤差（"iid"）、2列目はロバスト標準誤差（"HC1"）で表示する
modelsummary(list("ivreg（通常の標準誤差）" = iv_reg,
"ivreg（ロバスト標準誤差）" = iv_reg),
vcov = list("iid", "HC1"),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE")| ivreg（通常の標準誤差） | ivreg（ロバスト標準誤差） | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 699.914*** | 699.914*** |
| (4.438) | (4.413) | |
| university | 502.015*** | 502.015*** |
| (14.102) | (14.243) | |
| Num.Obs. | 10000 | 10000 |
| R2 | 0.477 | 0.477 |
| R2 Adj. | 0.477 | 0.477 |
| Std.Errors | IID | HC1 |
表の最後の行（Std.Errors）に、各列の標準誤差の種類が表示されている。2列目の標準誤差は、上のiv_robust()の結果と一致する。
結果の表をhtmlファイルとして保存するには、modelsummary()のoutput引数にファイル名を指定する。
# 結果表の出力：html
modelsummary(list("ivreg（通常の標準誤差）" = iv_reg,
"iv_robust（ロバスト標準誤差）" = iv_reg_robust),
stars = TRUE,
gof_omit = "AIC|BIC|Log.Lik|RMSE",
output = "06ivreg_table.html")6 参考文献
- 矢内勇生「計量経済学応用」（2020年度、高知工科大学）
- 第1回・第4回・第5回でも紹介した講義資料（日本語、Rのコードつき）。「Topic 9 操作変数法」が本講義の内容と対応している。
- Nick Huntington-Klein (2021) The Effect: An Introduction to Research Design and Causality（Chapman & Hall/CRC、ウェブ版は無料、英語）
- 「Chapter 19 - Instrumental Variables」が本講義の内容と対応している。
- Scott Cunningham Causal Inference: The Remix（ウェブ版は無料、英語）
- “Causal Inference: The Mixtape”（2021年）の第2版。「Instrumental Variables」の章が本講義の内容と対応しており、弱操作変数の問題などもより詳しく解説されている。Rのコードつき。
7 （おまけ）Pythonで同様の作業を行う場合
本講義の演習はRで行うが、参考として、ここまで学んだ内容に対応するPythonのコードを簡単に紹介する。
乱数の生成にはNumPy、通常の回帰分析にはstatsmodels、操作変数法（二段階最小二乗法）にはlinearmodelsを使う。なお、PythonとRでは乱数の生成方法が異なるので、同じ設定でデータを生成しても数値は少し異なる。
# ターミナル上で実行（Rのinstall.packages()に相当）
# pip install numpy pandas statsmodels linearmodels matplotlib seabornimport numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS, compare
import matplotlib.pyplot as plt
import seaborn as sns7.1 データの用意
rng = np.random.default_rng(10101) # 乱数の種を固定（Rの set.seed(10101) に相当）
n = 10000
# 能力は0から100まで均等に分布（Rの runif() に相当）
df = pd.DataFrame({"ID": np.arange(1, n + 1),
"ability": rng.uniform(0, 100, n)})
# 能力が80以上の人は、50%の確率で大卒になる（np.where()はRの if_else() に相当）
df["university"] = np.where(df["ability"] >= 80, rng.binomial(1, 0.5, n), 0)
# 能力に関係なく、30%の確率で学費免除になる（Rの rbinom(n, 1, 0.3) に相当）
df["exemption"] = rng.binomial(1, 0.3, n)
# まだ大卒でない学費免除者は、50%の確率で大卒になる
df["university"] = np.where((df["university"] == 0) & (df["exemption"] == 1),
rng.binomial(1, 0.5, n), df["university"])
# DGPに基づいて所得を生成（誤差項は平均0、標準偏差50の正規分布）
df["income"] = 200 + 10 * df["ability"] + 500 * df["university"] + rng.normal(0, 50, n)
# 記述統計（university・exemptionの平均値meanは、大卒・学費免除の人の割合）
df.describe()# 散布図
df["edu_label"] = np.where(df["university"] == 1, "Grad.", "Not grad.")
sns.scatterplot(data=df, x="ability", y="income", hue="edu_label", alpha=0.5)
plt.show()7.2 欠落変数を無視した回帰分析
reg_biased1 = smf.ols("income ~ university", data=df).fit()
reg_biased2 = smf.ols("income ~ university + exemption", data=df).fit()
print(reg_biased1.summary())
print(reg_biased2.summary())大卒の係数は約690・約749となり、R版と同様に過大推定される。
7.3 誘導型推定値を第一段階推定値で割る
first_stage = smf.ols("university ~ exemption", data=df).fit() # 第一段階
reduced_form = smf.ols("income ~ exemption", data=df).fit() # 誘導型
# 操作変数推定値 = 誘導型の係数 / 第一段階の係数
reduced_form.params["exemption"] / first_stage.params["exemption"]第一段階の係数は約0.44、誘導型の係数は約216で、操作変数推定値は約487となる（R版とは乱数が異なるため、数値も少し異なる）。
7.4 二段階最小二乗法(TSLS)と弱操作変数の確認
linearmodelsのIV2SLS.from_formula()では、[処置変数 ~ 操作変数]のように角括弧の中に第一段階を書く。また、切片を入れるために1 +を明示する必要がある（Rのincome ~ university | exemptionに相当）。
iv_reg = IV2SLS.from_formula("income ~ 1 + [university ~ exemption]",
data=df).fit(cov_type="unadjusted") # 通常の標準誤差
print(iv_reg.summary)
# 弱操作変数の確認（第一段階の結果。「Partial F-statistic」が第一段階のF値）
print(iv_reg.first_stage)
# Wu-Hausman検定
print(iv_reg.wu_hausman())操作変数推定値は約487（標準誤差は約15）、第一段階の\(F\)値は約2,956となり、十分大きい。
7.5 ロバスト標準誤差
cov_type="robust"とするとロバスト標準誤差になる。さらにdebiased=Trueとすると、Rのse_type = "HC1"と同じように自由度の調整が行われる。compare()で、通常の標準誤差の結果と並べて表示できる。
iv_reg_robust = IV2SLS.from_formula("income ~ 1 + [university ~ exemption]",
data=df).fit(cov_type="robust", debiased=True)
# 結果の比較表（括弧内は標準誤差）
print(compare({"Unadjusted": iv_reg, "Robust": iv_reg_robust},
precision="std_errors", stars=True))