内生性とは、一般的には \[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i (1) \] というモデルで、 \[ \mathrm{E}(\epsilon_i | X_i) \neq 0 \] という状態をいう（ここでは横断的データの分析を考えているが、考え方はパネル・データでも同じ）。ただし\(E(a | b)\) は条件付きの期待値で、\(b\)が真である（実現している）という条件下での\(a\)の期待値である。 つまり内生性とは、回帰モデルで独立変数の値によって残差がプラスになったり、マイナスになったりしやすくなるということである。これは独立変数と残差の間に相関があるということとは少し違うが、とりあえず同じ意味で理解しておいて構わない。内生性は傾きや切片の推定値をゆがませる原因になる。時間的に変化しない観察されない異質性 (Time Constant Unobservable Heterogeneity: TCUH) も内生性を生む一因である。
しかし、内生性は TCUH 以外の原因によっても生み出される。例えば、時間的に変化する観察されない異質性もあるし、独立変数と従属変数の間に双方向の因果関係がある場合にも内生性が生じる。これらは固定効果モデルを推定するだけではコントロールできない。
双方向の因果関係が内生性を生じさせることを簡単に解説しておこう。仮に上の (1) 式が真のモデルで、同時に以下のような関係が成り立っているとする。 \[ X_i = \gamma_0 + \gamma_1 Y_i + \upsilon_i (2) \] つまり、X と Y の間には双方向の因果関係があるということである。式(2) に式(1) を代入すると、 \[ X_i = \gamma_0 + \gamma_1 (\beta_0 + \beta_1 X_i + \epsilon_i) + \upsilon_i (3) \] となる。\(\epsilon_i\) に比例して \(X_i\) も変化するので、両者に相関が生じる。それゆえ \[ \mathrm{E}(\epsilon_i | X_i) \neq 0 \] である。直感的に考えても、両者に双方個の因果関係がある場合、それを考慮しないと X から Y への効果と Y から X への効果が混ざり合って、推定がゆがむ（両方ともプラスまたは両方ともマイナスなら効果の大きさを過大に推定するだろうし、一方がプラスで他方がマイナスならば打ち消しあって過少に推定する）のは理解できるだろう。
内生性による歪みを取り除く方法の一つに、操作変数法がある。操作変数とは、従属変数に直接効果を持たないが、内生性による歪みを取り除きたい独立変数には強い効果を持つような変数のことである。
今、上の式 (1) のようなモデルを推定したいとする。しかし、 \(X\) と \(\epsilon\) の間には関連があり、内生性があるとする。\(X\) には \(Z\) という変数が影響を与えているが、\(Z\) は \(X\) を媒介して \(Y\) に影響するものの直接的な影響を持っていないとする。それゆえ、 \(Z\) と \(\epsilon\) の相関はゼロ。このような前提が成り立つならば、以下のような計算をすることで \(\beta_1\) の一致推定量が得られる（つまり、\(\beta_1\) の真の値と、\(\beta_1\) のデータからの推定値の期待値が一致する）。
\(X\) を \(Z\) に回帰させて OLS で傾きと切片を推定する。
上で推定した切片と傾きから、\(X\) の予測値 \(\hat X\) を計算する。
\(Y\) を \(\hat X\) に回帰させて、傾きと切片を OLS で推定する。このとき得られる傾きが \(\beta_1\)の一致推定量になる。
このような推定法を二段階最小二乗法 (two stage least square estimation) と呼ぶ。以下は計算例。
# 架空のデータの準備
N <- 1000 # サンプル・サイズ
epsilon <- rnorm(N) # 残差
Z <- 1:N # 操作変数
X <- epsilon + 0.01 * Z + rnorm(N, sd=2) # 残差と相関する独立変数
Y <- 2 + 0.1 * X + epsilon # 従属変数。Xの真の係数は 0.1
d0 <- data.frame(Y, X, Z)
# 変数の相関のチェック
round(cor(d0), 2)## Y X Z
## Y 1.00 0.6 0.31
## X 0.60 1.0 0.80
## Z 0.31 0.8 1.00# 2段階最小二乗法の推定
lm.x <- lm(X ~ Z) # 一段階目のOLS
hatX <- predict(lm.x) # 上のモデルからの予測値を計算
tsls.y <- lm(Y ~ hatX) # 2段階目の OLS
lm.y <- lm(Y ~ X) # 比較のために単純に Y を X に回帰させた結果も
library(texreg)
screenreg(list(lm.y, tsls.y), single.row = TRUE)##
## ===================================================
## Model 1 Model 2
## ---------------------------------------------------
## (Intercept) 1.48 (0.05) *** 1.81 (0.07) ***
## X 0.19 (0.01) ***
## hatX 0.13 (0.01) ***
## ---------------------------------------------------
## R^2 0.36 0.10
## Adj. R^2 0.36 0.10
## Num. obs. 1000 1000
## RMSE 0.96 1.14
## ===================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05しかし、上のような計算では、係数は正確に推定されるものの標準誤差の推定が不正確である。そのため、専用の関数を使うのが一般的である。 R の場合、 AER パッケージの ivreg() 関数や sem パッケージの tsls() が使える。
ivreg(モデル式 | 操作変数, …) # “| 操作変数” 以外は lm() と同じ
tsls(モデル式, ~ 操作変数, …) # “~ 操作変数” 以外は lm() と同じ
library(AER)
tsls.y2 <- ivreg(Y ~ X | Z)
summary(tsls.y2)##
## Call:
## ivreg(formula = Y ~ X | Z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.01255 -0.67496 0.01523 0.68229 2.66230
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.81481 0.06065 29.92 <2e-16 ***
## X 0.12654 0.01060 11.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9943 on 998 degrees of freedom
## Multiple R-Squared: 0.3189, Adjusted R-squared: 0.3182
## Wald test: 142.6 on 1 and 998 DF, p-value: < 2.2e-16library(sem)
tsls.y3 <- tsls(Y ~ X, ~Z)
summary(tsls.y3) ##
## 2SLS Estimates
##
## Model Formula: Y ~ X
##
## Instruments: ~Z
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.01300 -0.67500 0.01523 0.00000 0.68230 2.66200
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.81481145 0.06064753 29.92391 < 2.22e-16 ***
## X 0.12653606 0.01059735 11.94035 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9942793 on 998 degrees of freedom一般的には、内生性のある説明変数 (X と略称)が複数ある場合もあるし、内生性の無い説明変数 (W と略称）も複数ありうる。操作変数 (Z と略称）も複数指定しうる。そのような場合、二段階最小二乗法を用いるためには、少なくとも個々の X に対して1つ以上の Z が必要である。そして、以下の手順で推定を行う。
例えば、内生性のある説明変数が２つ、内生性のない説明変数が２つ、操作変数が３つという場合を架空データで計算してみよう。
# 架空のデータの準備
N <- 1000 # サンプル・サイズ
epsilon <- rnorm(N) # 残差
W1 <- 1:N # 内生性のない独立変数
W2 <- scale(W1)^2 # 〃
Z1 <- rnorm(N) # 操作変数
Z2 <- scale(W1) + rnorm(N) # 〃 操作変数は W と相関していてもよい
Z3 <- -0.5 * scale(W2) + rnorm(N) # 〃
X1 <- epsilon - 2 * Z1 + 0.5 * Z2 + 0.005 * W1 + rnorm(N, sd=2) # 残差と相関する独立変数
X2 <- 2 * epsilon + 1.5 * Z2 + 0.3 * Z3 - 0.5 * W2 + rnorm(N, sd=2) # 〃
Y <- 2 + 0.1 * X1 + 0.003*W1 + 0.5* W2 + epsilon # 従属変数。Xの真の係数は 0.1
d0 <- data.frame(Y, X1, X2, W1,W2, Z1, Z2, Z3)
round(cor(d0), 2)## Y X1 X2 W1 W2 Z1 Z2 Z3
## Y 1.00 0.70 0.62 0.66 0.27 -0.10 0.49 -0.13
## X1 0.70 1.00 0.42 0.57 0.03 -0.53 0.48 -0.02
## X2 0.62 0.42 1.00 0.42 -0.18 0.05 0.58 0.13
## W1 0.66 0.57 0.42 1.00 0.00 -0.01 0.71 -0.01
## W2 0.27 0.03 -0.18 0.00 1.00 -0.03 0.00 -0.44
## Z1 -0.10 -0.53 0.05 -0.01 -0.03 1.00 0.03 0.00
## Z2 0.49 0.48 0.58 0.71 0.00 0.03 1.00 -0.04
## Z3 -0.13 -0.02 0.13 -0.01 -0.44 0.00 -0.04 1.00lm.x1 <- lm(X1 ~ W1 + W2 + Z1 + Z2 + Z3) # 第1ステップ： W と Z に X を回帰させる
lm.x2 <- update(lm.x1, X2 ~.) # 〃
hatX1 <- predict(lm.x1) # 第2ステップ： 第一ステップのモデルから X の予測値を計算
hatX2 <- predict(lm.x2) # 〃
tsls.y <- lm(Y ~ hatX1 + hatX2 + W1 + W2) # 第3ステップ：第2ステップの予測値を X の代わりに使って OLS
lm.y <- lm(Y ~ X1 + X2 + W1 + W2) # 比較のために通常の OLS でも計算
tsls.y2 <- ivreg(Y ~ X1 + X2 + W1 + W2 | W1 + W2 + Z1 + Z2 + Z3)
screenreg(list(lm.y, tsls.y, tsls.y2), digits=3,
custom.coef.names = c(rep(NA, 5), "X1", "X2"))##
## =====================================================
## Model 1 Model 2 Model 3
## -----------------------------------------------------
## (Intercept) 2.553 *** 2.018 *** 2.018 ***
## (0.059) (0.103) (0.092)
## X1 0.158 *** 0.081 *** 0.081 ***
## (0.008) (0.019) (0.017)
## X2 0.191 *** 0.011 0.011
## (0.008) (0.024) (0.022)
## W1 0.002 *** 0.003 *** 0.003 ***
## (0.000) (0.000) (0.000)
## W2 0.615 *** 0.490 *** 0.490 ***
## (0.028) (0.045) (0.040)
## -----------------------------------------------------
## R^2 0.784 0.515 0.613
## Adj. R^2 0.783 0.513 0.612
## Num. obs. 1000 1000 1000
## RMSE 0.763 1.143 1.020
## =====================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05以上のように操作変数があれば内生性の問題はすべて解消するので、固定効果モデルも必要ない。しかし、現実には操作変数を見つけるのは困難である。なぜなら、操作変数候補が従属変数に対して直接的な効果を持つと想定することはしばしば可能だからである。また、データから帰納法的に操作変数を見つけることも容易ではない。Z1, Z2, Z3 は、X1 と X2 を媒介して間接的にしか Y に影響しないように私が作ったものだが、上の Y を予測する単純な OLS モデルに操作変数である Z1, Z2, Z3 を投入すると有意になってしまうのである。
lm.y2 <- update(lm.y, ~.+ Z1 + Z2 + Z3) # 操作変数もモデルに投入してOLS
summary(lm.y2) # 内生性のせいで操作変数も直接効果が有意になってしまう！##
## Call:
## lm(formula = Y ~ X1 + X2 + W1 + W2 + Z1 + Z2 + Z3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.16339 -0.46408 0.02466 0.42668 2.35175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0681212 0.0616329 33.555 < 2e-16 ***
## X1 0.2168245 0.0100213 21.636 < 2e-16 ***
## X2 0.2287587 0.0078151 29.271 < 2e-16 ***
## W1 0.0023486 0.0001152 20.384 < 2e-16 ***
## W2 0.6004861 0.0264839 22.674 < 2e-16 ***
## Z1 0.2487609 0.0295948 8.406 < 2e-16 ***
## Z2 -0.3821291 0.0233537 -16.363 < 2e-16 ***
## Z3 -0.0870341 0.0217490 -4.002 6.76e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6618 on 992 degrees of freedom
## Multiple R-squared: 0.8378, Adjusted R-squared: 0.8366
## F-statistic: 731.7 on 7 and 992 DF, p-value: < 2.2e-16これは X1, X2 が内生性を持っているために生じたバイアスである。それゆえ、操作変数は理論上、間接的にしか影響しないはずだ、という想定を持てるような場合にしか使えないのである。
また、操作変数は、内生性を持つ説明変数と強く相関している必要がある。相関が弱いと二段階最小二乗法の推定結果の誤差が大きくなるのである。これらの条件を満たす変数を横断的データで見つけるのはしばしば困難である。
しかし、パネル・データでは一時点前の X を操作変数として用いるという方法が考えられる。例えば既婚女性のサンプルに関して、性役割意識と働いているかどうか（就業）の間に相関があったとしよう。固定効果モデルで、性役割意識が就業に影響を持っているという結果が出ても、因果の向きは逆ではないのか、という疑いが当然生じる。そこで、操作変数として一時点前の性役割意識を使うことが考えられる。一年前の性役割意識と現在の性役割意識の間にはある程度の相関が期待できる。そして、一年前の性役割意識が直接現在の就業に影響を及ぼすことはないと考えられるならば、二段階最小二乗法が使えるということである。
このようにパネル・データは操作変数を使う上でも有効なデータであると考えられる。
plm() は操作変数を指定することができる。書式は ivreg() と同じである。以下は計算例。
library(plm)
data("EmplUK")
d1 <- pdata.frame(EmplUK)
fixed1 <- plm(wage ~ emp + lag(capital) + lag(output), data=d1) # 操作変数なし
iv1 <- update(fixed1, ~. | # ".-emp" で "lag(capital) + lag(output)"
. -emp + lag(emp)) # と指定した場合と同じ
iv2 <- update(fixed1, ~. | # lag(emp, 2) で2時点前の emp の値を返す
. -emp + lag(emp) + lag(emp, 2), data=d1) # 操作変数を増やしたモデル
screenreg(list(fixed1, iv1, iv2), single.row = TRUE, digits=3)##
## ===========================================================================
## Model 1 Model 2 Model 3
## ---------------------------------------------------------------------------
## emp -0.109 (0.031) *** -0.082 (0.059) -0.120 (0.051) *
## lag(capital) 0.204 (0.058) *** 0.200 (0.058) *** 0.145 (0.066) *
## lag(output) -0.089 (0.008) *** -0.091 (0.009) *** -0.098 (0.008) ***
## ---------------------------------------------------------------------------
## R^2 0.213 0.213 0.271
## Adj. R^2 0.179 0.178 0.219
## Num. obs. 891 891 751
## ===========================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05fixed3 <- plm(emp ~ lag(emp) + lag(emp, 2), data=d1) # 内生変数と操作変数の関係を確認
summary(fixed3)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = emp ~ lag(emp) + lag(emp, 2), data = d1)
##
## Unbalanced Panel: n=140, T=5-7, N=751
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -12.40000 -0.14300 -0.00276 0.13400 25.20000
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## lag(emp) 0.873987 0.042359 20.6330 < 2.2e-16 ***
## lag(emp, 2) -0.214848 0.043234 -4.9694 8.739e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 3369.6
## Residual Sum of Squares: 1979.6
## R-Squared : 0.4125
## Adj. R-Squared : 0.33451
## F-statistic: 213.802 on 2 and 609 DF, p-value: < 2.22e-16