Stata：处理衡量偏误-变量误差模型的一些建议
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 面板数据因果推断专题
作者：王舒瑶 (吉林大学)
邮箱：378807478@qq.com
[编者按]： 本文介绍的内容来自如下文献，特此致谢！
[Source]：Lockwood J. R., Daniel F. McCaffrey, 2020, Recommendations about estimating errors-in-variables regression in Stata, Stata Journal, 20(1): 116–130. -PDF-, -PDF2-, -Link-
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 引言
2. 背景
3. EIV 回归
3.1 模型假设
3.2 EIV 回归估计量
3.3 标准误估计
3.4 eivreg标准误的偏差
3.5 Stata 命令用法
4. 模拟研究
4.1 模拟研究概述
4.2 Stata 模拟程序和过程
5. 结论
参考资料
1. 引言
内生性有 3 个主要来源：「遗漏变量」、「互为因果」和「衡量偏误」。其中第三种「衡量偏误」受到的关注相对较少，但若存在这一问题，会导致研究者无法一致地估计解释变量的系数。
在以往的推文 第三种内生性：衡量偏误(测量误差)如何检验-dgmtest？中，我们介绍了如何判断是否存在「衡量偏误」。而在推文 xtewreg：面板数据存在衡量偏误-测量偏误时如何估计？中，我们进一步对「衡量偏误」的应对办法 xtewreg 进行探讨。
在本文中，我们对 eivreg 和 sem 这两种常见的变量误差 (errors-in-variables, EIV) 回归方法进行比较，并给出建议。
2. 背景
举个例子来导入：在可以用来估计教育政策对学生成就的影响的传统教育生产函数模型中，当前的研究成果是兴趣的投入和先前的成就之间的函数。 但是，先前的成就是不能被观测到的。 相反，它是通过测试分数来衡量的，该分数通常是从美国各州或学区进行的标准化评估中获得的。 使用易于出错的度量代替相应的潜在变量来拟合模型通常会导致所有模型参数的估计值不一致，而不仅仅是与误差测量的变量相对应的回归系数。
在已知有关衡量误差文献中，我们通常使用两种主要方法来获得一致的估计量：eivreg 和 sem 。 eivreg 使用矩量法调整来解决测量中的误差。 sem 是为回归模型中的因变量、自变量和测量误差指定联合高斯分布，然后使用最大似然估计回归系数。
3. EIV 回归
3.1 模型假设
对 i = 1 到 N，让 ( , , , , ) 与具有有限第四矩的分布是独立且均匀分布的 (IID) 。
, 和 是长度均为p的向量。
随机变量 和 是标量。
, 是可以观测到的数据，而 , , 是不可以观测到的数据。
我们的目标是估计出 ，但问题是“真实”的 是错误地由 估计的。
3.2 EIV 回归估计量
如果测量误差 的方差-协方差矩阵是已知的或可以估计的，则 EIV 回归估计使用矩量法来估计 。在 3.1 的模型假设下，假设其他标准正则条件成立，则 的 EIV 回归估计量是一致的。
我们将注意力集中在 是具有元素的对角矩阵的情况下，因此 的不同分量中的测量误差是相互不相关的。 我们专注于这种情况，因为此假设通常在应用程序中进行，并且 eivreg 要求此限制。 在 不是对角矩阵的情况下，仍然可以很好地定义 EIV 回归估计量，在这种情况下，Stata 用户可以使用 sem 而不是 eivreg ，因为 sem 的语法足够通用，可以允许 是非对角矩阵。
3.3 标准误估计
如上所述，在假设成立的情况下， eivreg 使用的矩方法算法和 sem 使用的最大似然算法在理论上得出相同的 值。原因是，两种算法都产生了一组相同的估计方程组，其解为 。
在实践中，只要存在给定的观测数据和假定的可靠性，并且只要 sem 使用的迭代算法收敛，那么这两个命令将报告相同的 值，直到数值上的差异很小。
但是，这些命令没有使用相同的方法来估计 的方差-协方差矩阵，而 eivreg 所使用的方法将倾向于低估在潜在变量建模中通常做出的假设下估算器的实际抽样方差。原因在于，对于下述公式 (3) ：
的替代表达式 (4)：
而 eivreg 估计出来的 是 , 是上述表达式 (4) 右侧的第二项，第一项被忽略了。
仅在所有协变量及其对应的测量误差均已固定的前提下， eivreg 报告的 才是适当的。 该假设通常与从总体中随机抽取单位不一致，并且在具有测量误差的应用中特别受限制，因为它也以未观察到的测量误差的固定值为条件。
此外，即使在需要保证这些假设的应用中，eivreg 使用的方差估计量也将适合于表征估计回归系数的抽样变异性，但不适用于表征这些估计量的均方误差。因为 通常不等于 。
另外，由 sem 报告的 估计值会产生标准误差估计值，该估计值与协变量的随机抽样、测量误差以及总体分布的结果一致，隐含地解释了上述 的替代表达式 (4) 右侧的两个项。
3.4 eivreg标准误的偏差
通常很难评估表达式 (4) 中两个项的相对大小，但是一些基本结果是显而易见的。 对于固定的 N 和固定的 , 收敛到零且 收敛到 ，即在所有变量的随机抽样下， 的 OLS 估计量的方差。
还要注意， 不依赖于 。因此，对于固定的 N 和
小于 1 的固定可靠性, 由于 趋近 0, 占主导地位。因为 eivreg 忽略了这一项，eivreg 报告的标准误差将趋向于负偏，并且相关的置信区间的覆盖率将小于名义水平，无论样本量大小。
3.5 Stata 命令用法
eivreg 在实际操作时与 reg 类似，如：
. eivreg y x c
sem ( structural equation modeling ) 是一种结构方程模型，它可以同时处理潜变量 ( 难以直接准确测量 ) 和显变量 ( 利用显变量去间接测量显变量 ) 。以单因素测量模型为例：
. sem ( x1 x2 x3 x4 <- X )
其中，x1 x2 x3 x4 为显变量，X 为潜变量。
4. 模拟研究
4.1 模拟研究概述
我们进行了简单的模拟研究，以证明 eivreg 报告的标准误差与 sem 报告的标准误差之间的实际差异。
我们考虑 （即截距和标量预测变量）的情况，并且仅关注预测变量的估计系数。 我们的模拟改变了三个因素：样本量 N ，真实回归的 和可靠性 r。
对于每个模拟条件和蒙特卡洛复制，我们使用模拟数据, 分别基于eivreg 和 sem 来计算 的 EIV 回归估计 b 及其相关的标准误差估计。 对于 eivreg ，我们既跟踪了 b 的报告标准差，又跟踪了使用 250 个自抽样估计的标准误。 我们使用250个自抽样的样本，因为根据 Efron 和 Tibshirani (1993, 52) 的大多数情况，这个数量应该足够了，而且计算时间并不是很长。
在对模拟研究的初步探索中，我们发现了具有 EIV 回归估计量并已通过 eivreg 成功计算但其 sem 并未从其默认起始值收敛到该解的模拟数据集。因此，我们修改了对 sem 的调用，以使用模型参数的 MLE 作为起始值。 以这种方式初始化参数可导致 sem 迅速收敛，并在 sem 和 eivreg 之间估计出回归系数，这表明所有模拟数据集的数值差异都可以忽略不计。
模拟结果与 eivreg 的负偏差分析结果一致。 使用 eivreg 报告的标准误， 的 95% 置信区间的覆盖率不到 95% 。 在 100 个模拟条件下，使用 eivreg 报告的标准误的覆盖概率范围为 0.58 至 0.95 ，平均值为 0.87 。 当 大且 r 小时，无论样本大小 N 如何，覆盖率都更差。当 小而 r 大时，无论 N 大小，覆盖率都接近名义水平。 eivreg 报告的平均标准误差与 b 的估计抽样标准误之比在 0.42 至 1.02 之间，平均值为 0.82 ，这与置信区间的掩盖度一致。
使用 sem 报告的标准误或自举标准误计算的置信区间更接近名义覆盖率。 对于 sem ，覆盖范围从 0.94 到大于 0.99 ，平均值为 0.97。 由 sem 报告的平均标准误差与 b 的估计抽样标准误之比在 0.95 至 1.74 之间，平均值为 1.14 ，与置信区间覆盖范围略有超出名义水平一致。 对于自举标准误，覆盖范围在 0.92 至 0.96 之间，平均值为 0.95 ，平均标准误差与 b 的估计采样标准误之比在 0.95 至 1.05 之间，平均值为 0.99。
4.2 Stata 模拟程序和过程
4.2.1 运行一次的模拟迭代程序
先使用 program drop 将原有程序从内存中删除。如果内存中原本不存在这个程序，那么 program drop 将返回一个错误。在 program drop 之前输入 capture 将捕获此错误并允许代码块继续运行。
然后定义我们此次的程序 simit ，定义程序后，输入 simit ，Stata 将在 program 和 end 之间运行所有命令。选项 rclass 告诉 Stata 我们想要使用 return 返回程序中的值。
接着用 syntax 定义程序的语法，nobs ( ) 表示观察值的计数，integer 表示需要输入整数；rsq ( ) 和 lambda ( ) 是参数，real 表示需要输入实数。
随后就都是对运行一次的模拟迭代程序的语句进行设计。
最后一行是 end ，它告诉 Stata 你已经完成了对程序 simit 的定义。
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
4.2.2 条件和蒙特卡罗复制的循环模拟
蒙特卡罗模拟是指通过大量随机抽样进行模拟，最终得出变量的累计概率分布。设置 seed 是为了让实验可以重复。主体循环是先取 obs_seq ，再取 rsq_seq ，最后取 lambda_seq ，例如： 第一次取 100 0.10 0.50 ，第二次取 100 0.10 0.60 ，…… ，第六次取 100 0.30 0.50 ，以此类推。 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
需要注意的是，4.2.1 和 4.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 。
5. 结论
变量误差（ EIV ）回归是在具有易错协变量的线性模型中进行一致估计的标准方法， Stata 命令 eivreg 和 sem 均可实现这一目的。 但是，在通常做出的假设下, eivreg 报告的标准误会负偏，从而导致置信区间覆盖率低于名义水平。 因此，在 EIV 回归的大多数实际应用中，单独使用 sem 或使用自抽样标准误增强的 eivreg 比单独使用 eivreg 更好。
参考资料
Lockwood J. R., Daniel F. McCaffrey, 2020, Recommendations about estimating errors-in-variables regression in Stata, Stata Journal, 20(1): 116–130. -PDF-, -PDF2-, -Link- Meijer, E., E. Oczkowski, T. Wansbeek, 2021, How measurement error affects inference in linear regression, Empirical Economics, 60 (1): 1-25. -PDF- 第三种内生性：衡量偏误(测量误差)如何检验-dgmtest？ xtewreg：面板数据存在衡量偏误-测量偏误时如何估计？
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。