两阶段最小二乘法与R
1、ivreg包介绍
R语言计量包ivreg用以解决线性回归模型的内生性问题。
描述：工具变量估计的线性模型通过两阶段最小二乘(2SLS) 回归或通过稳健回归M估计(2SM)或MM估计(2SMM)。主要的ivreg()模型拟合函数旨在提供一个工作流程，尽可能类似于标准的lm()回归。大量的方法是被用来拟合ivreg模型对象，除了其他标准模型工具，还包括广泛的功能，计算和图形回归诊断。
作者：Author John Fox, Christian Kleiber, Achim Zeileis
2 、操作流程
在利用工具变量法估计线性回归模型时，往往选择Stata操作，现在介绍R的操作方法：首先，我们先安装工具变量回归安装包ivreg，并加载相关其他计量包；
setwd("D:/Allcode/Rstudy/model/IV_estimate") # 先设置路径
install.packages("ivreg") # 安装ivreg
install.packages("haven") # 用于stata数据导入，默认存在，可以不安装
install.packages("lmtest") # 用于线性回归检验
install.packages("sandwich") # 提供相关异方差稳健标准误
#加载以上所有包
library("haven")
library("ivreg")
library("lmtest")
library("sandwich")
接下来准备数据集，我选用的是陈强老师主页（陈强教授的计量经济学及Stata主页 (econometrics-stata.com)）上的数据集grilic，它是stata的dta格式，因此需要转换导入
grilic <- read_dta("grilic.dta")
names(grilic) # 查看数据框的变量名
# [1] "rns" "rns80" "mrt" "mrt80" "smsa" "smsa80"
# [7] "med" "iq" "kww" "year" "age" "age80"
# [13] "s" "s80" "expr" "expr80" "tenure" "tenure80"
# [19] "lw" "lw80"
该数据集中包括以下变量:lw(工资对数),s (受教育年限) , age(年龄) , expr(工龄) , tenure(在现单位的工作年数),q(智商), med(母亲的受教育年限),kww(在"knowledge of the World ofWork"测试中的成绩),mt(婚姻虚拟变量,已婚=1),rns (美国南方虚拟变量,住在南方=1),smsa(大城市虚拟变量,住在大城市=1),year (有数据的最早年份, 1966-1973年中的某一年)。我们选择lw为被解释变量，其余变量为解释变量。先利用OLS回归作为基准模型
fit_ols1 <- lm(lw ~ s+expr+tenure+rns+smsa,data = grilic) # 没有加入智商iq变量
# Call:
# lm(formula = lw ~ s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
# Min 1Q Median 3Q Max
# -1.11684 -0.22626 -0.01511 0.23103 1.23738
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 4.103675 0.085097 48.223 < 2e-16 ***
# s 0.102643 0.005849 17.549 < 2e-16 ***
# expr 0.038119 0.006327 6.025 2.65e-09 ***
# tenure 0.035615 0.007742 4.600 4.96e-06 ***
# rns -0.084080 0.028797 -2.920 0.00361 **
# smsa 0.139667 0.028082 4.974 8.15e-07 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3464 on 752 degrees of freedom
# Multiple R-squared: 0.3521, Adjusted R-squared: 0.3478
# F-statistic: 81.75 on 5 and 752 DF, p-value: < 2.2e-16
#由于没有加入iq，存在遗漏变量问题，因此加入iq
fit_ols2 <- lm(lw ~ iq+s+expr+tenure+rns+smsa,data = grilic)
# Call:
# lm(formula = lw ~ iq + s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
# Min 1Q Median 3Q Max
# -1.16056 -0.21786 -0.00622 0.22771 1.20580
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 3.895172 0.109110 35.699 < 2e-16 ***
# iq 0.003279 0.001083 3.028 0.00255 **
# s 0.092787 0.006666 13.920 < 2e-16 ***
# expr 0.039344 0.006306 6.239 7.33e-10 ***
# tenure 0.034209 0.007715 4.434 1.06e-05 ***
# rns -0.074532 0.028815 -2.587 0.00988 **
# smsa 0.136737 0.027948 4.893 1.22e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3445 on 751 degrees of freedom
# Multiple R-squared: 0.36, Adjusted R-squared: 0.3548
# F-statistic: 70.39 on 6 and 751 DF, p-value: < 2.2e-16
# 以上回归都是基于同方差假设条件下的结果，我们将系数转换为异方差稳健标准误；这里以fit_ols2为例
coeftest(fit_ols2, vcov = vcovHC, type = "HC1") # 异方差稳健标准误
# t test of coefficients:
#
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 3.8951718 0.1159286 33.5997 < 2.2e-16 ***
# iq 0.0032792 0.0011321 2.8965 0.003883 **
# s 0.0927874 0.0069763 13.3004 < 2.2e-16 ***
# expr 0.0393443 0.0066603 5.9072 5.272e-09 ***
# tenure 0.0342090 0.0078957 4.3326 1.674e-05 ***
# rns -0.0745325 0.0299772 -2.4863 0.013124 *
# smsa 0.1367369 0.0277712 4.9237 1.045e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
使用工具变量法回归，内生解释变量为iq，工具变量选择med、kww、mrt、age;其余控制变量自身视为自身的工具变量；代码如下：
fit_iv <- ivreg(lw ~ iq+s+expr+tenure+rns+smsa |
s+ med+ kww+mrt+age+expr+tenure+rns+smsa ,data = grilic)
# 这里用"|"分隔内生解释变量与工具变量
# 提取稳健标准误
coeftest(fit_iv, vcov = vcovHC, type = "HC0") # 异方差稳健标准误
# t test of coefficients:
#
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 4.8378747 0.3799432 12.7332 < 2.2e-16 ***
# iq -0.0115468 0.0056376 -2.0482 0.040887 *
# s 0.1373477 0.0174989 7.8489 1.446e-14 ***
# expr 0.0338041 0.0074844 4.5166 7.295e-06 ***
# tenure 0.0405640 0.0095848 4.2321 2.602e-05 ***
# rns -0.1176984 0.0359582 -3.2732 0.001112 **
# smsa 0.1499830 0.0322276 4.6539 3.850e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 工具变量法回归还要进行诊断
summary(fit_iv,test = TRUE) # 诊断
# Call:
# ivreg(formula = lw ~ iq + s + expr + tenure + rns + smsa | s +
# med + kww + mrt + age + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
# Min 1Q Median 3Q Max
# -1.3825405 -0.2437078 0.0009735 0.2514625 1.4609417
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 4.837875 0.346424 13.965 < 2e-16 ***
# iq -0.011547 0.005241 -2.203 0.027889 *
# s 0.137348 0.017042 8.059 3.02e-15 ***
# expr 0.033804 0.007302 4.630 4.32e-06 ***
# tenure 0.040564 0.008896 4.560 5.98e-06 ***
# rns -0.117698 0.035468 -3.318 0.000949 ***
# smsa 0.149983 0.031572 4.751 2.43e-06 ***
#
# Diagnostic tests:
# df1 df2 statistic p-value
# Weak instruments 4 748 10.54 2.61e-08 *** # 弱工具变量检验（通过）
# Wu-Hausman 1 750 10.70 0.00112 ** # 内生性检验（通过）
# Sargan 3 NA 61.14 3.36e-13 *** # 过度识别检验（未通过）
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3851 on 751 degrees of freedom
# Multiple R-Squared: 0.2002, Adjusted R-squared: 0.1938
# Wald test: 55.92 on 6 and 751 DF, p-value: < 2.2e-16
由于工具变量个数大于内生解释变量个数，且工具变量过度识别检验未通过，因此需要调整工具变量；这里怀疑age与tenure可能存在过度识别，剔除后进行ivreg回归
fit_iv2 <- ivreg(lw ~ iq+s+expr+tenure+rns+smsa |
med+ kww + s + expr + tenure + rns + smsa ,data = grilic)
summary(fit_iv2,test = TRUE)
# Call:
# ivreg(formula = lw ~ iq + s + expr + tenure + rns + smsa | med +
# kww + s + expr + tenure + rns + smsa, data = grilic)
#
# Residuals:
# Min 1Q Median 3Q Max
# -1.3025533 -0.2405658 0.0005969 0.2349962 1.2621665
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 3.218043 0.384814 8.363 2.97e-16 ***
# iq 0.013928 0.005884 2.367 0.018186 *
# s 0.060780 0.018735 3.244 0.001230 **
# expr 0.043324 0.007038 6.156 1.22e-09 ***
# tenure 0.029644 0.008561 3.463 0.000565 ***
# rns -0.043527 0.034922 -1.246 0.213000
# smsa 0.127222 0.030137 4.221 2.72e-05 ***
#
# Diagnostic tests:
# df1 df2 statistic p-value
# Weak instruments 2 750 14.906 4.49e-07 *** （通过）
# Wu-Hausman 1 750 3.858 0.0499 * （通过）
# Sargan 1 NA 0.130 0.7185 （通过）
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.3661 on 751 degrees of freedom
# Multiple R-Squared: 0.2775, Adjusted R-squared: 0.2718
# Wald test: 61.94 on 6 and 751 DF, p-value: < 2.2e-16
参考文献
陈强.高级计量经济学[M].高等教育出版社
https://cran.r-project.org/web/packages/gmm/gmm.pdf
http://www.econometrics-stata.com/col.jsp?id=101
1255
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
