B-Stata模拟：控制变量！控制变量！Good-Controls-Bad-Controls

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
  . ssc install lianxh
详情参见帮助文件 (有惊喜)：
  . help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg
 

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：李烨阳 (浙江大学)
邮箱：li_yeyang95@163.com

编者按：本文主要整理自下文，特此致谢！
Source：Cinelli, C., A. Forney, J. Pearl, 2022, A crash course in good and bad controls, Sociological Methods & Research, Forthcoming. -PDF- -Link- -R- -Stata-

本系列分成理论和实操两部分：

• A-理论部分：控制变量！控制变量！Good-Controls-Bad-Controls
• B-Stata模拟：控制变量！控制变量！Good-Controls-Bad-Controls

目录

1. 导言

在连享会推文「控制变量！控制变量！Good-Controls-Bad-Controls」中，我们已经介绍了如何使用因果图 (有向无环图，DAG) 筛选控制变量的理论知识，以及几种社会科学研究中较为典型的良好控制、糟糕控制和中性控制情形。在本文中，我们将使用 Stata 软件，对前文提到的 18 中模型进行模拟分析，以检验理论是否成立。

2. 良好的控制变量

为了方便读者阅读，本文模拟分析的顺序与「控制变量！控制变量！Good-Controls-Bad-Controls」一致，模型的编号与原文一致。

Stata 社区发布的 dag 程序包提供了对有向无环图进行模拟分析的外部命令。但由于这一程序包主要是面向流行病学研究设计的，它生成的模拟变量为 0-1 变量。而在社会科学研究中，我们常常需要使用到连续变量，同时为与原文的 R 语言模拟分析保持一致，本文并没有直接使用该程序包，感兴趣的读者可以自行尝试。

2.1 解释变量和被解释变量的混淆情形

模型 1-3 是存在同时影响解释变量和被解释变量的混淆的情形，我们已经知道，图中的 $Z$ 都是良好的控制变量。

首先我们对模型 1 进行模拟分析。

第一步，先根据 DAG 图中变量间的关系，定义模型 1 的模拟程序：每次模拟我们随机生成 10000 个样本，$z$ 是一个符合 $N(0,1)$ 标准正态分布的随机变量，$X$ 等于 $Z$ 加上一个符合 $N(0,1)$ 标准正态分布的随机变量，$Y$ 等于 $X$ 加 $Z$ 再加上一个符合 $N(0,1)$ 标准正态分布的随机变量。然后进行两次回归估计，第一次不控制变量 $Z$，第二次控制变量 $Z$，并分别记录每次回归 $X$ 系数估计结果的点估计和 95% 置信区间的上界和下界。

*- step 1 定义模型的数值模拟程序

	* model 1 模拟程序
	capture program drop model_1 // 模型命名为model_1, 若有更换注意修改
	program model_1,rclass       // 若有更换注意修改
	clear
	set obs 10000

	* 此部分为针对不同模型结构需要修改的部分
	***************************
	* model 1: DAG结构
	gen z = rnormal()
	gen x = z + rnormal()
	gen y = x + z + rnormal()
	***************************

	* 不控制 z
	reg y x
	mat A=r(table)
	return scalar mean0=A[1,1]
	return scalar low0=A[5,1]
	return scalar up0=A[6,1]
    
	* 控制 z
	reg y x z 
	mat A=r(table)
	return scalar mean1=A[1,1]
	return scalar low1=A[5,1]
	return scalar up1=A[6,1]
	end

第二步，进行重复模拟分析。我们进行 100 次模拟，为了保证结果可复现，设定了种子值。

*- step 2 模拟模块
	clear
	set seed 11
	simulate mean0=r(mean0) low0=r(low0) up0=r(up0) ///
	mean1=r(mean1) low1=r(low1) up1=r(up1), reps(100): model_1 // 对 model_1 进行 100 次模拟, 若有更换注意修改
	gen id=_n

第三步，我们对模拟结果进行可视化。左图为不控制变量 $Z$ 的情形，右图为控制变量 $Z$ 的情形。横坐标是 100 次模拟的序列，纵坐标是 $X$ 系数的估计值。蓝色附加线为 100 次模拟的 $X$ 系数点估计值的均值，红色附加线为真实的数据模拟设定的 $X$ 系数，本文中所有模型均设置为 1。

*- step 3 模拟结果可视化

	set scheme s1color //设定绘图模板
	sum mean0
	local b_mean `r(mean)'
	twoway (scatter mean0 id, ms(Oh) mc(gs5)) (rcap low0 up0 id, lc(gs5)), ///
	ytitle("E(β{sub:x}_hat)") legend(order(1 "E（β{sub:x}_hat）点估计" 2 "E（β{sub:x}_hat）置信区间")) ///
	yline(`b_mean', lcolor(blue)) 
	graph save model_1a.gph, replace //图名model_1a，若有更换注意修改

	sum mean1
	local b_mean `r(mean)'
	twoway (scatter mean1 id, ms(Oh) mc(gs5)) (rcap low1 up1 id, lc(gs5)), ///
	ytitle("E(β{sub:x}_hat|z)") legend(order(1 "E（β{sub:x}_hat|z）点估计" 2 "E（β{sub:x}_hat|z）置信区间")) ///
	yline(`b_mean', lcolor(blue))  yline(1, lcolor(red)) 
	graph save model_1b.gph, replace  //图名model_1b，若有更换注意修改
	
	graph combine model_1a.gph  model_1b.gph , /// //若有更换注意修改图名
	xsize(10) ysize(5) title("Model 1") //图形合并，图标题名Model 1，若有更换注意修改

结果如下图所示。可以看出，在不控制 $Z$ 的情况下，100 次模拟的估计系数均值约为 1.5，且从图中可以看出，95% 置信区间均没有包括真实系数，估计系数是有偏的。而引入了控制变量 $Z$ 之后，100 次模拟的估计系数均值与真实值非常接近，且绝大部分的单次模拟的 95% 置信区间均包含真实值。因此数值模拟结果支持了模型 1 中的 $Z$ 是良好的控制变量这一观点。

为了节约篇幅，在后续的模型数值模拟介绍中，只给出步骤 1 定义模型的数值模拟程序步骤中 DAG 结构部分的代码，读者只需要将这一部分代码替换到上述例子的 step 1 中同样用一组星号分割开的位置中即可，并更改相应的模型命名、图名和图的标题名，其余部分毋须改动。当然，大家也可以通过点击「good_bad_controls_simulation.do」，下载完整的模拟代码。

模型 2 DAG 结构部分的代码为：

	***************************
	* model 2: DAG 结构
	gen u = rnormal()
	gen z = u + rnormal()
	gen x = z + rnormal()
	gen y = x + u + rnormal()
	***************************

模型 2 模拟结果为如下图所示，同样印证了 $Z$ 是一个好的控制变量。

模型 3 DAG 结构部分的代码为：

	***************************
	* model 3: DAG结构
	gen u = rnormal()
	gen z = u + rnormal()
	gen x = u + rnormal()
	gen y = x + z + rnormal()
	***************************

其结果同样与理论一致。

2.2 解释变量和中介变量的混淆情形

模型 4、5、6 是存在同时影响解释变量和中介变量的混淆的情形，此三种情况同样是良好的控制变量。

下面分别给出模型 4、5、6 DAG 结构部分的代码，模拟分析结果均显示，不控制 $Z$ 时存在混淆偏差，而控制了 $Z$ 之后都可以得到 $X$ 系数无偏的估计。

	***************************
	* model 4: DAG 结构 
	gen z = rnormal()
	gen x = z + rnormal()
	gen m = x + z + rnormal()
	gen y = m + rnormal()
	***************************

	***************************
	* model 5: DAG 结构 
	gen u = rnormal() 
	gen z = u + rnormal()
	gen x = z + rnormal()
	gen m = x + u + rnormal()
	gen y = m + rnormal()
	***************************

	***************************
	* model 6: DAG 结构 
	gen u = rnormal() 
	gen z = u + rnormal()
	gen x = u + rnormal()
	gen m = x + z + rnormal()
	gen y = m + rnormal()
	***************************

模型 4-6 的模拟结果为：

2.3 处理后变量一定是坏的控制变量吗

虽然大部分情况下处理后变量是坏的控制变量，但也有例外。例如模型 15 中，如果以 $W$ 为条件是给定的前提，此时 $Z$ 是一个好的控制变量。

定义模型 15 的数值模拟程序为：

*- step 1 定义模型的数值模拟程序

	* model 15 模拟程序
	capture program drop model_15 
	program model_15,rclass 
	clear
	set obs 10000
	
	***************************
	* model 15: DAG 结构 
	gen x = rnormal() 
	gen u = rnormal() 
	gen z = x + rnormal()
	gen w = z + u + rnormal()
	gen y = x -2*u + rnormal()
	***************************

	* 控制 w, 不控制 z
	reg y x w
	mat A=r(table)
	return scalar mean0=A[1,1]
	return scalar low0=A[5,1]
	return scalar up0=A[6,1]

	* 控制 w 和 z
	reg y x w z 
	mat A=r(table)
	return scalar mean1=A[1,1]
	return scalar low1=A[5,1]
	return scalar up1=A[6,1]
	end

模拟结果如下图所示，左图为只控制 $W$，此时会造成选择性偏误，而当控制 $W$ 的同时也控制 $Z$，由于控制对撞变量而开放的路径又重新截断了，可以获得无偏的估计。

3. 糟糕的控制变量

3.1 处理前变量一定是好的控制变量吗

模型 7 是一个 $M$ 偏误的典型例子，这里的 $Z$ 虽然是处理前变量，但却是糟糕的变量。

模型 7 (a) DAG 结构部分的代码为：

	***************************
	* model 7: DAG结构 
	gen u1 = rnormal() 
	gen u2 = rnormal() 
	gen z = u1 + u2 + rnormal()
	gen x = u1 + rnormal()
	gen y = x -4*u2 + rnormal()
	***************************

模拟结果如下图所示，由于 $Z$ 是一个对撞变量，在不控制 $Z$ 时可以获得无偏估计，控制 $Z$ 反而会引入非因果的相关性。

模型 10 同样是一个控制变量虽然是处理前变量，但仍然是糟糕的控制变量的例子。

模型 10 DAG 结构部分的代码为：

	***************************
	* model 10: DAG结构 
	gen z = rnormal() 
	gen u = rnormal() 
	gen x = -2*z + u + rnormal()
	gen y = x + 2*u + rnormal()
	***************************

模拟结果如下图所示，在不增加控制变量 $Z$ 时，由于有混淆会造成估计偏差，100 次模拟的估计系数均值约为 1.3，但是这里引入控制变量 $Z$ 后，不仅不能截断后门路径，还会放大偏差，100 次模拟的估计系数均值约为 2.0。

3.2 警惕阻断因果路径

模型 11、12 是由于控制中介变量或其子变量造成过度控制偏差的典型范例。

模型 11、12 DAG 结构部分的代码分别为：

	***************************
	* model 11: DAG 结构 
	gen x = rnormal() 
	gen z = x + rnormal() 
	gen y = z + rnormal()
	***************************

    ***************************
	* model 12: DAG 结构 
	gen x = rnormal() 
	gen m = x + rnormal() 
	gen z = m + rnormal() 
	gen y = m + rnormal()
	***************************

模拟结果分别如下图所示，在不控制中介变量或其子变量时，估计是无偏的，但引入他们作为控制变量则会导致间接因果效应路径被截断，从而低估结果，因此中介变量或其子变量是糟糕的控制变量。

3.3 警惕打开对撞变量所在路径

模型 16 和模型 17 中的 $Z$ 均为对撞变量，如果控制对撞变量将引入选择性偏差，因此是糟糕的控制变量。

模型 16、17 DAG 结构部分的代码分别为：

	***************************
	* model 16: DAG 结构 
	gen x = rnormal() 
	gen u = rnormal() 
	gen z = x + u + rnormal() 
	gen y = x + 2*u + rnormal()
	***************************

	***************************
	* model 17: DAG 结构 
	gen x = rnormal() 
	gen y = x + rnormal()
	gen z = x + y + rnormal()
	***************************

模拟结果分别为：

控制模型 18 中的 $Z$ 会导致案例控制偏差。

模型 18 DAG 结构部分的代码为：

	***************************
	* model 18: DAG 结构 
	gen x = rnormal() 
	gen y = x + rnormal()
	gen z = y + rnormal()
	***************************

模拟结果为：

4. 中性的控制变量

4.1 可能提高精度的情形

模型 8 中的 $Z$ 是一个中性的控制变量，但控制它可以提高估计结果的精确度。

模型 8 DAG 结构部分的代码为：

	***************************
	* model 8: DAG 结构 
	gen z = rnormal() 
	gen x = rnormal()
	gen y = x + 2*z + rnormal()
	***************************

从模拟结果可以看出，虽然无论是否控制 $Z$，都可以获得无偏估计，但是控制 $Z$ 可以提高估计精度，即系数估计值的标准误降低，系数置信区间范围缩小。

模型 13 中的 $Z$ 同样是一个提高估计精度的中性的控制变量。

模型 13 DAG 结构部分的代码为：

	***************************
	* model 13: DAG 结构 
	gen z = rnormal() 
	gen x = rnormal()
	gen m = 2*z + rnormal()
	gen y = x + 2*m + rnormal()
	***************************

模拟结果为：

4.2 可能降低精度的情形

模型 9 和 14 中的 $Z$ 均为中性的控制变量。但控制它们会减少 $X$ 的变动性，从而降低估计结果的精确度。

模型 9 和 14 DAG 结构部分的代码分别为：

	***************************
	* model 9: DAG 结构 
	gen z = rnormal() 
	gen x = 2*z + rnormal()
	gen y = x + 2*rnormal()
	***************************

	***************************
	* model 14: DAG 结构 
	gen x = rnormal() 
	gen z = 2*x + rnormal()
	gen y = x + 2*rnormal()
	***************************   

从模拟结果可以看出，虽然无论是否控制 $Z$，都可以获得无偏估计，但是控制 $Z$ 会估计精度，即系数估计值的标准误增加，系数置信区间范围增大。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 控制变量, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：论文写作

  • 控制变量如何选？大牛们的10条建议
• 专题：Stata命令

  • Stata：控制变量组合的筛选-tuples
  • Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？
• 专题：回归分析

  • 控制变量越多越好吗？
  • Stata：控制变量与核心解释变量地位对等吗？
  • 调节效应是否需要考虑对控制变量交乘？
  • 控制变量！控制变量！
  • 不用太关心控制变量，真的！
  • 加入控制变量后结果悲催了！
• 专题：IV-GMM

  • Lasso一下：再多的控制变量和工具变量我也不怕-T217
• 专题：断点回归RDD

  • RDD：断点回归可以加入控制变量吗？
  • Stata：RDD-中可以加入控制变量
• 专题：内生性-因果推断

  • 控制变量！控制变量！Good-Controls-Bad-Controls
  • 敏感性分析B-Stata实操：控制变量内生时的系数敏感性分析-regsensitivity
  • 敏感性分析A-理论基础：控制变量内生时的系数敏感性分析-regsensitivity
• 专题：其它

  • 锚定情境法（一）：有效控制变量自评偏差

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版
• Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 连享会-主页 和 知乎专栏，700+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis

New！ lianxh 和 songbl 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh