Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:王舒瑶 (吉林大学)
邮箱:378807478@qq.com
编者按:本文主要整理自以下文献,特此致谢!
Source:Lockwood J R, McCaffrey D F. Recommendations about estimating errors-in-variables regression in Stata[J]. The Stata Journal, 2020, 20(1): 116-130. -PDF-
目录
内生性有 3 个主要来源:遗漏变量、互为因果和衡量偏误。其中第三种衡量偏误受到的关注相对较少,但若存在这一问题,会导致研究者无法一致地估计解释变量的系数。
在连享会之前的推文「第三种内生性:衡量偏误(测量误差)如何检验-dgmtest?」中,我们介绍了如何判断是否存在衡量偏误。在推文「xtewreg:面板数据存在衡量偏误-测量偏误时如何估计?」中,我们进一步对衡量偏误的估计方法进行了探讨。
在已知有关衡量误差的文献中,我们主要使用两种方法来获得一致的估计量:eivreg
和 sem
。其中,eivreg
使用矩量法调整来解决测量中的误差。sem
通过为因变量、自变量和测量误差指定联合高斯分布,使用最大似然估计来解决测量中的误差。在本文中,我们将对 eivreg
和 sem
这两种常见的变量误差 (errors-in-variables,EIV) 回归方法进行比较,并给出建议。
令
我们的目标是估计出
如果测量误差
我们将注意力集中在 eivreg
要求此限制。
在 sem
而不是 eivreg
。因为 sem
的语法足够通用,可以允许
如上所述,在假设成立的情况下,eivreg
使用的矩方法算法和 sem
使用的最大似然算法在理论上得出相同的
在实践中,只要存在给定的观测数据和假定的可靠性,并且只要 sem
使用的迭代算法收敛,那么这两个命令将报告相同的
但是,这些命令没有使用相同的方法来估计 eivreg
所使用的方法倾向于低估在潜在变量建模中通常做出的假设下估算器的实际抽样方差。原因在于,对于以下表达式:
eivreg
估计出来的
仅在所有协变量及其对应的测量误差均已固定的前提下,eivreg
报告的
此外,即使在需要保证这些假设的应用中,eivreg
使用的方差估计量也将适合于表征估计回归系数的抽样变异性,但不适用于表征这些估计量的均方误差。因为
另外,由 sem
报告的
通常很难评估替代表达式中两个项的相对大小,但是一些基本结果是显而易见的。对于固定的 N 和固定的
还要注意,eivreg
忽略了这一项,eivreg
报告的标准误差将趋向于负偏,并且相关的置信区间的覆盖率将小于名义水平,无论样本量大小。
eivreg
在实际操作时与 reg
类似,如:
eivreg y x c
sem
(structural equation modeling) 是一种结构方程模型,它可以同时处理潜变量 (难以直接准确测量) 和显变量 (利用显变量去间接测量显变量)。以单因素测量模型为例:
sem (x1 x2 x3 x4 <- X) // x1 x2 x3 x4 为显变量, X 为潜变量
我们进行了简单的模拟研究,以证明 eivreg
报告的标准误差与 sem
报告的标准误差之间存在的实际差异。具体地,我们考虑
对于每个模拟条件和蒙特卡洛复制,我们使用模拟数据,分别基于 eivreg
和 sem
来计算 eivreg
,我们既跟踪了
在对模拟研究的初步探索中,我们发现了具有 EIV 回归估计量并已通过 eivreg
成功计算,但其 sem
并未从其默认起始值收敛到该解的模拟数据集。因此,我们修改了对 sem
的调用,以使用模型参数的 MLE 作为起始值。以这种方式初始化参数可导致 sem
迅速收敛,并在 sem
和 eivreg
之间估计出回归系数,这表明所有模拟数据集的数值差异都可以忽略不计。
模拟结果与 eivreg
的负偏差分析结果一致。使用 eivreg
报告的标准误,eivreg
报告的标准误覆盖概率范围为 0.58 至 0.95,平均值为 0.87。当 eivreg
报告的平均标准误差与
使用 sem
报告的标准误或自举标准误计算的置信区间更接近名义覆盖率。对于 sem
,覆盖范围从 0.94 到 0.99,平均值为 0.97。由 sem
报告的平均标准误差与
先使用 program drop
将原有程序从内存中删除。如果内存中原本不存在这个程序,那么 program drop
将返回一个错误。在 program drop
之前输入 capture
将捕获此错误并允许代码块继续运行。
然后定义程序 simit
,定义程序后,输入 simit
,Stata 将在 program
和 end
之间运行所有命令。选项 rclass
告诉 Stata 我们想要使用 return
返回程序中的值。
接着用 syntax
定义程序的语法,nobs()
表示观察值的计数,integer
表示需要输入整数,rsq()
和 lambda()
是参数,real
表示需要输入实数。
随后就都是对运行一次的模拟迭代程序的语句进行设计。
capture program drop simit
/* *********************************************** */
/* program for running one iteration of simulation */
/* *********************************************** */
program simit, rclass
syntax, nobs(integer) rsq(real) lambda(real)
preserve
set obs `nobs'
/* generate data */
generate double xstar = rnormal(0.0, 1.0)
generate double err = rnormal(0.0, sqrt( (1.0 - `rsq') / `rsq' ))
generate double u = rnormal(0.0, sqrt( (1.0 - `lambda') / `lambda'))
generate double xobs = xstar + u
generate y = 0.0 + 1.0*xstar + err
/* proceed if EIV is possible given observed data and reliability */
quietly correlate y xobs
if (`lambda' > r(rho)^2) {
return scalar eiv_ok = 1
/* run -eivreg- with reported standard errors ("eiv") */
eivreg y xobs, reliab(xobs `lambda')
local b_init = _b[xobs]
scalar cieiv_l = _b[xobs]-1.96*_se[xobs]
scalar cieiv_u = _b[xobs]+1.96*_se[xobs]
return scalar b_eiv = _b[xobs]
return scalar se_eiv = _se[xobs]
return scalar cover_eiv = cond(cieiv_l<1 & cieiv_u>1,1,0)
/* run -eivreg- with bootstrapped standard errors ("beiv") */
bootstrap, reps(250): eivreg y xobs, reliab(xobs `lambda')
scalar cibeiv_l = _b[xobs]-1.96*_se[xobs]
scalar cibeiv_u = _b[xobs]+1.96*_se[xobs]
return scalar b_beiv = _b[xobs]
return scalar se_beiv = _se[xobs]
return scalar cover_beiv = cond(cibeiv_l<1 & cibeiv_u>1,1,0)
/* run -sem-, initializing regression coefficient at -eivreg- solution, */
/* and initializing var(X) and var(Y|X) to their MLEs */
quietly summarize xobs
local vX_init = r(Var) * ((`nobs' - 1)/`nobs') * `lambda'
quietly summarize y
local veY_init = (r(Var) * ((`nobs' - 1)/`nobs')) - ///
(`b_init'*`b_init'*`vX_init')
capture sem (xobs <- X) (y <- (X, init(`b_init'))), ///
var((X, init(`vX_init'))) var((e.y, init(`veY_init'))) ///
reliab(xobs `lambda')
matrix B = e(b)
matrix vB = e(V)
scalar b_sem = B[1,3]
scalar se_sem = sqrt(vB[3,3])
scalar cisem_l = b_sem -1.96*se_sem
scalar cisem_u = b_sem +1.96*se_sem
return scalar b_sem = b_sem
return scalar se_sem = se_sem
return scalar cover_sem = cond(cisem_l<1 & cisem_u>1,1,0)
return scalar sem_status = _rc
}
else {
return scalar eiv_ok = 0
}
restore
end
蒙特卡罗模拟是指通过大量随机抽样进行模拟,最终得出变量的累计概率分布。设置 seed
是为了让实验可以重复。主体循环是先取 obs_seq,再取 rsq_seq,最后取 lambda_seq。例如,第一次取 100 0.10 0.50,第二次取 100 0.10 0.60,以此类推。nsim 表示循环重复 1800 遍。
关于蒙特卡罗模拟的一些基本介绍,可以参看连享会推文「stata中计算公式命令_Stata:学点蒙特卡洛模拟分析」。
/* ************************************************************ */
/* loop simulation over conditions and Monte Carlo replications */
/* ************************************************************ */
set more off
set seed 1417
set linesize 140
local nobs_seq 100 500 1000 5000
local rsq_seq 0.10 0.30 0.50 0.70 0.90
local lambda_seq 0.50 0.60 0.70 0.80 0.90
local nsim 1800
local filename results_all
tempname simulation
postfile `simulation' numok numsemconv nobs rsq lambda ///
b_eiv_mean b_eiv_sd se_eiv_mean cover_eiv ///
b_beiv_mean b_beiv_sd se_beiv_mean cover_beiv ///
b_sem_mean b_sem_sd se_sem_mean cover_sem ///
using `filename', replace
display "$S_TIME $S_DATE"
foreach nobs of local nobs_seq {
foreach rsq of local rsq_seq {
foreach lambda of local lambda_seq {
simulate eiv_ok=r(eiv_ok) sem_status=r(sem_status) ///
b_eiv=r(b_eiv) se_eiv=r(se_eiv) cover_eiv=r(cover_eiv) ///
b_beiv=r(b_beiv) se_beiv=r(se_beiv) cover_beiv=r(cover_beiv) ///
b_sem=r(b_sem) se_sem=r(se_sem) cover_sem=r(cover_sem), ///
reps(`nsim'): simit, nobs(`nobs') rsq(`rsq') lambda(`lambda')
quietly summarize eiv_ok
scalar numok = r(sum)
generate tmp = (sem_status==0)
quietly summarize tmp
scalar numsemconv = r(sum)
drop tmp
/* compute summary statistics to save, */
/* keeping cases where EIV possible and -sem- converged. */
/* also check that estimated coefficients are the same */
keep if ((sem_status==0) & (eiv_ok==1))
generate d = b_sem - b_eiv
summarize d
drop d
foreach var of varlist b_eiv se_eiv cover_eiv b_beiv se_beiv cover_beiv ///
b_sem se_sem cover_sem {
quietly summarize `var'
scalar `var'_mean=r(mean)
}
foreach var of varlist b_eiv b_beiv b_sem {
quietly summarize `var'
scalar `var'_sd=r(sd)
}
post `simulation' (numok) (numsemconv) (`nobs') (`rsq') (`lambda') ///
(b_eiv_mean) (b_eiv_sd) (se_eiv_mean) (cover_eiv_mean) ///
(b_beiv_mean) (b_beiv_sd) (se_beiv_mean) (cover_beiv_mean) ///
(b_sem_mean) (b_sem_sd) (se_sem_mean) (cover_sem_mean)
clear
}
}
}
postclose `simulation'
display "$S_TIME $S_DATE"
use "`results_all'", clear
sort nobs rsq lambda
list
outsheet using "results_all.csv", comma nolabel replace
需要注意的是,3.2.1 和 3.2.2 是一个整体,如果想要跑出最终的模拟结果,两部分命令缺一不可。其中前面部分的目的在于设计 simit
的 program,后面部分的目的在于运用这个 program 来完成蒙特卡洛模拟。
按照目前的参数设置 (4 个 nobs_seq,5 个 rsq_seq,5 个 lambda_seq),该程序一共跑 100 组 (4 * 5 * 5),每组模拟 1800 次,最后生成 100 条数据,所以 Stata 跑起来还是很耗费时间的。由于 100 个结果展示出来篇幅过于大,笔者在这里只展示 10 个模拟结果 (nobs_seq 取值 100,rsq_seq 取值 0.10、0.30,lambda_seq 取值 0.50、0.60、0.70、0.80、0.90) ,结果如下表所示:
nobs | rsq | lambda | b_eiv_mean | b_eiv_sd | se_eiv_mean | cover_eiv | b_beiv_mean | b_beiv_sd | se_beiv_mean | cover_beiv | b_sem_mean | b_sem_sd | se_sem_mean | cover_sem |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
100 | 0.1 | 0.5 | 0.9688179 | 0.506548 | 0.4473057 | 0.8888889 | 0.9688179 | 0.506548 | 0.4717829 | 0.8888889 | 0.9688179 | 0.5065481 | 0.4821762 | 0.9444444 |
100 | 0.1 | 0.6 | 0.9425998 | 0.4650904 | 0.3935203 | 0.8888889 | 0.9425998 | 0.4650904 | 0.3995831 | 0.9444444 | 0.9425998 | 0.4650904 | 0.4109256 | 0.9444444 |
100 | 0.1 | 0.7 | 1.089905 | 0.3589229 | 0.3425317 | 0.8888889 | 1.089905 | 0.3589229 | 0.3497599 | 0.8888889 | 1.089905 | 0.3589229 | 0.3543588 | 0.9444444 |
100 | 0.1 | 0.8 | 0.976297 | 0.3039494 | 0.347137 | 1 | 0.976297 | 0.3039494 | 0.3498529 | 0.9444444 | 0.976297 | 0.3039494 | 0.3493564 | 1 |
100 | 0.1 | 0.9 | 0.9823672 | 0.2514369 | 0.3138203 | 1 | 0.9823672 | 0.2514369 | 0.3139614 | 1 | 0.9823671 | 0.2514369 | 0.3129175 | 1 |
100 | 0.3 | 0.5 | 1.056897 | 0.2012456 | 0.2108894 | 1 | 1.056897 | 0.2012456 | 0.2190944 | 0.9444444 | 1.056896 | 0.2012456 | 0.2800636 | 1 |
100 | 0.3 | 0.6 | 0.9567788 | 0.1621872 | 0.1953815 | 0.9444444 | 0.9567788 | 0.1621872 | 0.2073208 | 0.9444444 | 0.9567788 | 0.1621872 | 0.228746 | 1 |
100 | 0.3 | 0.7 | 0.9748029 | 0.1727375 | 0.1817665 | 1 | 0.9748029 | 0.1727375 | 0.1892546 | 1 | 0.9748029 | 0.1727375 | 0.2006758 | 1 |
100 | 0.3 | 0.8 | 0.9440436 | 0.1934716 | 0.1777468 | 0.8888889 | 0.9440436 | 0.1934716 | 0.1886845 | 0.9444444 | 0.9440436 | 0.1934716 | 0.1856588 | 0.8888889 |
100 | 0.3 | 0.9 | 1.05179 | 0.1580933 | 0.1671066 | 0.8888889 | 1.05179 | 0.1580933 | 0.1640482 | 0.8888889 | 1.05179 | 0.1580933 | 0.1700332 | 0.8888889 |
表格的前三列显示的是 nobs_seq、rsq_seq 和 lambda_seq 的取值情况,随后依次展示 eivreg
,自抽样标准误增强的 eivreg
和 sem
的系数估计结果,模拟结果显示这三者对系数大小的估计情况差不多,但是 eivreg
的标准误会偏低,置信区间的覆盖率更小。自抽样标准误增强的 eivreg
可以在一定程度上纠正 eivreg
的偏差,但是效果依旧不如 sem
。
变量误差 (EIV) 回归是在具有易错协变量的线性模型中进行一致估计的标准方法, Stata 命令 eivreg
和 sem
均可实现这一目的。但是,在通常做出的假设下,eivreg
报告的标准误会负偏,从而导致置信区间覆盖率低于名义水平。 因此,在 EIV 回归的大多数实际应用中,单独使用 sem
或使用自抽样标准误增强的 eivreg
比单独使用 eivreg
更好。
Note:产生如下推文列表的 Stata 命令为:
lianxh 内生性 衡量偏误, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh