31. 重回帰分析（操作変数法）

このページは、浅野正彦先生の「計量分析（政治）」補助教材をRマークダウンで学習した記録です。
　　　31. 重回帰分析 13（操作変数法：IV 基礎）
乱数モデルを用いた操作変数法の解説は秀逸です。この場を借りて御礼申し上げます。

0.1 x2iと ziの生成

まず N = 1000 で交絡因子 x2と操作変数 zを生成する.
交絡因子 x2と操作変数zは無関係なので、rnorm()関数を使って、それぞれ個別に生成する.

set.seed(1)
n1 <- 1000　　　　# N = 1000 と設定
x2 <- rnorm(n1)　 # N = 1000 の x2 を生成する   
z <- runif(n1)   # N = 1000 の z を生成する 
cor(x2,z) #無関係であるが無相関までは言っていない

## [1] 0.03504146

操作変数 Zは X1に対して直接影響を及ぼす
操作変数 Zは Yに対して直接影響を及ぼさない
操作変数 Zは Yに対して直接影響を及ぼさず、X1に対してのみ直接影響を及ぼすため、緑色の点線が単独で結果変数とオーバーラップするところはない

0.2 操作変数に関する仮定のまとめ

操作変数の外生性 (instrument exogeneity)

cov[Z,U]=0
Ziに向けて X2iからも Yiからも矢印が出ていない　= Ziと X2iは関係がない

操作変数の関連性 (instrument relevance)

cov[Z,X1]≠0
Ziから X1iに矢印が向いている　= Ziと X1iは関係がある

除外制約 (exclusion restriction)

操作変数 Ziと誤差項 Uiと無相関である -Ziから X1iに矢印が出ている
Ziから Yiには矢印が出ていない　→ Ziは Yiに間接的に影響を与えている　→ Xiの値が決まれば、 Ziがあってもなくても、Xiから Yiへの影響は変わらない
操作変数推定量 β1を求める式　\(\beta_1 = \frac{cov[Z, Y]}{cov[Z, X_1]}\)

0.3 誤差項 e1を生成する

e1 <- rnorm(n1)
mean(e1)

## [1] -0.01021977

0.4 x1iの生成

誤差項 e1を加味して、x2iと ziから影響を受けている共変量 x1iを生成する

x1 <- 1 + 2 * z - 3 * x2 + e1

0.5 yiの生成

#誤差項 e2を生成する
e2 <- rnorm(n1)
#観察されない共変量 x2iをモデルに含めない時の誤差項 u1を作成する
u1 <- x2 + e2
#x2iと x1iから影響を受けている共変量 yiを生成する
y1 <- 3 + 1.5 * x1 + 1.5 * x2 + e2

par(mfrow = c(2,2))
plot(x1, y1)
plot(x2, y1)
plot(x2, z)

操作変数 Ziと観察されない共変量 x2iをモデルに含めない時の誤差項 u1の共分散がゼロであることを確認

cor(z, u1)

## [1] 0.01051698

cor(x2, u1)

## [1] 0.7009528

共変量 x1と観察されない共変量 x2iをモデルに含めない時の誤差項 u1の共分散がゼロでないことの確認.

cor(x1, u1)

## [1] -0.6537806

0.6 シミュレーション結果

Yi	結果変数	観測可
x1i	処置変数（共変量）	観察可
x2i	共変量（＝交絡因子）	観察不可
zi	操作変数	観察可

0.6.1 重回帰モデル（真の値）

m_1 <- lm(y1 ~ x1 + x2)
tidy(m_1)

term	estimate	std.error	statistic	p.value
(Intercept)	3.07	0.0631	48.7	1.47e-265
x1	1.49	0.0273	54.8	8.16e-303
x2	1.47	0.0868	17	7.8e-57

par(mfcol=c(1,2))
plot(m_1,which=c(1,2))

0.6.2 単回帰モデル

m_2 <- lm(y1 ~ x1)
tidy(m_2)

term	estimate	std.error	statistic	p.value
(Intercept)	3.91	0.0439	89.1	0
x1	1.06	0.0114	93.4	0

par(mfcol=c(1,2))
plot(m_2,which=c(1,2))

0.6.3 操作変数モデル

m_3 <- cov(z, y1)/cov(z, x1)
tidy(m_3)

## Warning:  'tidy.numeric' は廃止予定です 
##  help("Deprecated") を見て下さい

x
1.57

0.6.4 シミュレーションの結果まとめ

・真の推定値は 1.49 ・単回帰モデルの推定値 (1.06) は真の推定値 (1.49) を過小評価している・操作変数モデルの推定値 (1.57) は真の推定値 (1.49) に近い値を推定している →　未観測の共分散 x2 の交絡を除去している

0.6.5 2 段階最小二乗法

step_1 <- lm(x1 ~ z)     # 第 1 段階の最小二乗法
x1hat <- predict(step_1) # x1 の予測値 x1hat を抽出  

step_2 <- lm(y1 ~ x1hat) # 第 2 段階の最小二乗法
tidy(step_2)　　　　　　 # 結果を表示

term	estimate	std.error	statistic	p.value
(Intercept)	2.9	0.477	6.08	1.68e-09
x1hat	1.57	0.231	6.78	2.03e-11

参考
操作変数法の解析手順【図解】頑健な妻