Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈卓然(中山大学)
邮箱:chenzhr25@mail2.sysu.edu.cn
目录
在因果推断中一个标准的问题,是识别和估计某个处理变量
尽管如此,实证研究者们还是想知道这一假设在实证分析中究竟有多么重要,或者说结论对于无混淆性因素假设有多么敏感?
定义如下的随机变量:
实际观测到的结果满足:
我们不能观测到全部的 DGP
在这一假设下,一个处理效应统计量的可识别集合将会是一个闭集,这一闭集依赖于 tesensitivity
命令就是用来计算这些区间,并展现这些处理效应统计量的可识别集合是如何随着敏感性参数
除了对这一系列的
命令安装:
ssc install tesensitivity, replace
我们使用的数据有两份:一份是实验数据,一份是可观测数据。有关数据的详细内容,请参考连享会推文「因果推断:混杂因素敏感性分析理论(上)」。
. net get tesensitivity, replace
. use lalonde1986, clear
. global outcome "re78"
. global treatment "treat"
. global controls "married age black hispanic education re74 re75 re74pos re75pos"
这份数据集中包含我们在理论篇中提到的两份数据。其中,sample1 指示的是实验数据,也就是 NSW 项目中的处理组和控制组;sample2 指示的是 NSW 中的处理组和从 PSID 中构建的控制组;sample3 是在 sample2 的基础上去除项目开始前 3 到 4 年工资高于 5000 美元的工人。
我们首先对于实验数据和非实验数据分别采用逆概率加权的方法估计平均处理效应。
. * 实验数据
. teffects ipw ($outcome) ($treatment $controls) if sample1
Treatment-effects estimation Number of obs = 445
Estimator : inverse-probability weights
Outcome model : weighted mean
Treatment model: logit
------------------------------------------------------------------------------
| Robust
re78 | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
treat |
(1 vs 0) | 1632.930 649.886 2.51 0.012 359.177 2906.683
-------------+----------------------------------------------------------------
POmean |
treat |
0 | 4578.176 343.789 13.32 0.000 3904.362 5251.990
------------------------------------------------------------------------------
. * 非实验数据:
. teffects ipw ($outcome) ($treatment $controls) if sample3
Treatment-effects estimation Number of obs = 390
Estimator : inverse-probability weights
Outcome model : weighted mean
Treatment model: logit
------------------------------------------------------------------------------
| Robust
re78 | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATE |
treat |
(1 vs 0) | 3336.676 769.180 4.34 0.000 1829.112 4844.241
-------------+----------------------------------------------------------------
POmean |
treat |
0 | 2868.201 331.833 8.64 0.000 2217.820 3518.582
------------------------------------------------------------------------------
* POmean 表示潜在产出均值 The POM for treatment level t is the average potential outcome for that treatment level
可以看到在无混淆因素假设下,处理效应在两个样本下都是显著为正。接着,我们使用 tesensitivity
分析这些结果对于假设有多么敏感。我们用到的子命令是 tesensitivity cpi
,cpi
就是条件部分独立。这一命令在给定一系列 c 值的前提下计算处理效应统计量边界,并计算处理效应统计量高于一个给定阈值这一结论的截断点。
首先我们计算 ATE 的边界,默认情况下,这一命令会给出 40 个 c 之下的边界和 ATE
. tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ate
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 445
Outcome model : linear quantile Breakdown = 0.075
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 1,419, 1,419]
0.026 | [ 920, 1,930]
0.051 | [ 433, 2,454]
0.077 | [ -43, 2,986]
0.103 | [ -508, 3,529]
...
0.897 | [-11,594, 21,952]
0.923 | [-11,594, 21,952]
0.949 | [-11,594, 21,952]
0.974 | [-11,594, 21,952]
1.000 | [-11,594, 21,952]
---------------------------------------
我们也可以通过 cpiplot
命令将上述结果直观的呈现出来。
. tesensitivity cpiplot
从上图中可以发现,当 c 相对较大时,这个边界范围变得非常宽,但是 c 应该多大算是比较合理的呢?cscale
命令对此做出了回答。这一命令主要实现的是理论篇中提到的去一分析法,即比较去除某一个协变量之后的倾向得分和包含全部协变量时的倾向得分之间的差异。它提供给我们一个 c 值的参考,也就是什么水平的 c 值是相对较为合理的。
因为上述这种倾向得分的差异是随着被丢掉的协变量不同而不同的,这也就产生了差异的分布,cscale
子命令可以帮助我们计算这一分布的上确界和一些分位数。默认情况下,tesensitivity cscale
计算每一个协变量产生分布的上确界和 50th、75th、90th 分位点,其他分位点可以通过 quantiles
选项得到。
. tesensitivity cscale
Treatment-effects sensitivity analysis
Analysis : leave one out prop. score diff. Number of obs = 445
Treatment model : logistic
--------------------------------------------------
quantile | 0.500 0.750 0.900 | max
---------------+-------------------------+--------
married | 0.006 0.012 0.032 | 0.042
age | 0.015 0.024 0.034 | 0.099
black | 0.007 0.009 0.014 | 0.082
hispanic | 0.007 0.017 0.099 | 0.124
education | 0.012 0.022 0.031 | 0.087
re74 | 0.002 0.011 0.035 | 0.209
re75 | 0.001 0.004 0.008 | 0.053
re74pos | 0.002 0.010 0.018 | 0.034
re75pos | 0.013 0.017 0.062 | 0.082
--------------------------------------------------
我们也可以通过加入 density
选项的方式,绘制某一个协变量造成的倾向得分差异的分布,比如教育。
. tesensitivity cscale education,density
正如在理论篇中所提及的,我们可以将去一分析法中得到的倾向得分偏离作为 c 的一个参考。具体而言,我们不妨取去一分析法中得到的最大偏离作为 c 的参考,当 c 超过这一个最大偏离就说明我们的结论是相当稳健的,反之就要怀疑结论的稳健性了。
. tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ate creference
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 445
Outcome model : linear quantile Breakdown = 0.075
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 1,419, 1,419]
0.026 | [ 920, 1,930]
re74pos 0.034 | [ 751, 2,109]
married 0.042 | [ 616, 2,254]
0.051 | [ 433, 2,454]
re75 0.053 | [ 398, 2,492]
0.077 | [ -43, 2,986]
re75pos 0.082 | [ -128, 3,083]
black 0.082 | [ -134, 3,090]
education 0.087 | [ -234, 3,206]
age 0.099 | [ -448, 3,458]
0.103 | [ -508, 3,529]
hispanic 0.124 | [ -888, 3,996]
0.128 | [ -963, 4,091]
0.154 | [ -1,411, 4,674]
0.179 | [ -1,851, 5,313]
0.205 | [ -2,282, 6,023]
re74 0.209 | [ -2,354, 6,146]
0.231 | [ -2,714, 6,765]
0.256 | [ -3,150, 7,554]
...
0.949 | [-11,594, 21,952]
0.974 | [-11,594, 21,952]
1.000 | [-11,594, 21,952]
---------------------------------------
当然我们也可以通过图形的方式将上表直观地呈现出来。
. tesensitivity cpiplot, creflines
tesensitivity cpi
的结果会被存储在 e()
中。首先我们可以通过 estat summarize
命令,查看 tesensitivity cpi
命令中用到的全部变量的描述性统计。
. estat summarize
Estimation sample tesensitivity Number of obs = 445
-------------------------------------------------------------------
Variable | Mean Std. dev. Min Max
-------------+-----------------------------------------------------
re78 | 5300.764 6631.492 0 60307.93
treat | .4157303 .4934022 0 1
married | .1685393 .3747658 0 1
age | 25.37079 7.100282 17 55
black | .8337079 .3727617 0 1
hispanic | .0876404 .2830895 0 1
education | 10.19551 1.792119 3 16
re74 | 2102.265 5363.582 0 39570.68
re75 | 1377.138 3150.961 0 25142.24
re74pos | .2674157 .4431092 0 1
re75pos | .3505618 .4776829 0 1
-------------------------------------------------------------------
此外注意全部协变量被存储在矩阵 e(covsupp)
中,注意这里包括常数项。我们也可以将估计结果通过 est store
这一命令存储下来,然后通过 est replay
的命令加以调用。
. est store ate1
. est replay ate1
----------------------------------------------------------------------------------
Model ate1
----------------------------------------------------------------------------------
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 445
Outcome model : linear quantile Breakdown = 0.075
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 1,419, 1,419]
0.026 | [ 920, 1,930]
re74pos 0.034 | [ 751, 2,109]
married 0.042 | [ 616, 2,254]
0.051 | [ 433, 2,454]
re75 0.053 | [ 398, 2,492]
0.077 | [ -43, 2,986]
re75pos 0.082 | [ -128, 3,083]
black 0.082 | [ -134, 3,090]
education 0.087 | [ -234, 3,206]
age 0.099 | [ -448, 3,458]
0.103 | [ -508, 3,529]
hispanic 0.124 | [ -888, 3,996]
0.128 | [ -963, 4,091]
0.154 | [ -1,411, 4,674]
0.179 | [ -1,851, 5,313]
0.205 | [ -2,282, 6,023]
re74 0.209 | [ -2,354, 6,146]
0.231 | [ -2,714, 6,765]
0.256 | [ -3,150, 7,554]
...
0.949 | [-11,594, 21,952]
0.974 | [-11,594, 21,952]
1.000 | [-11,594, 21,952]
---------------------------------------
我们也可以将不同的回归综合起来进行比较,譬如将上述实验组的数据和下面非实验组的数据进行比较。
. tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample3, ate
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 390
Outcome model : linear quantile Breakdown = 0.020
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 2,865, 2,865]
0.026 | [ -663, 11,806]
0.051 | [ -2,879, 15,385]
...
0.949 | [-13,374, 27,392]
0.974 | [-13,376, 27,397]
1.000 | [-13,376, 27,399]
---------------------------------------
. est sto ate3
我们可以使用 tesensitivity cpitable
命令来进行简单的比较。
. tesensitivity cpitable ate1 ate3
Treatment effects sensitivity
Analysis : conditional partial independence, multiple results
------------------------------------------------------------
c | ate1 | ate3
-------------------+--------------------+-------------------
0.000 | [ 1,419, 1,419] | [ 2,865, 2,865]
0.026 | [ 920, 1,930] | [ -663, 11,806]
re74pos 0.034 | [ 751, 2,109] | [ ., .]
married 0.042 | [ 616, 2,254] | [ ., .]
0.051 | [ 433, 2,454] | [ -2,879, 15,385]
re75 0.053 | [ 398, 2,492] | [ ., .]
0.077 | [ -43, 2,986] | [ -4,389, 17,794]
re75pos 0.082 | [ -128, 3,083] | [ ., .]
black 0.082 | [ -134, 3,090] | [ ., .]
education 0.087 | [ -234, 3,206] | [ ., .]
age 0.099 | [ -448, 3,458] | [ ., .]
0.103 | [ -508, 3,529] | [ -5,356, 19,300]
hispanic 0.124 | [ -888, 3,996] | [ ., .]
0.128 | [ -963, 4,091] | [ -6,121, 20,318]
0.154 | [ -1,411, 4,674] | [ -6,764, 21,055]
0.179 | [ -1,851, 5,313] | [ -7,296, 21,687]
0.205 | [ -2,282, 6,023] | [ -7,746, 22,256]
re74 0.209 | [ -2,354, 6,146] | [ ., .]
0.231 | [ -2,714, 6,765] | [ -8,232, 22,747]
0.256 | [ -3,150, 7,554] | [ -8,678, 23,234]
...
0.949 | [-11,594, 21,952] | [-13,374, 27,392]
0.974 | [-11,594, 21,952] | [-13,376, 27,397]
1.000 | [-11,594, 21,952] | [-13,376, 27,399]
------------------------------------------------------------
Analysis : breakdown point, multiple
-------------------------------------------
conclusion | ate1 ate3
------------------+------------------------
stat > 0 | 0.075 0.020
-------------------------------------------
当然我们也可以将上表直观的呈现出来。
. tesensitivity cpiplot ate1 ate3
tesensitivity
可以支持的处理效应统计量有:
对于处理组的平均处理效应 ATET,我们下面比较一下两个样本中 ATET 的敏感性。
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, atet
. qui estimates store atet1
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample3, atet
. qui estimates store atet3
. tesensitivity cpiplot atet1 atet3
对于分位数处理效应,我们必须通过 quantiles
这一选项来指定所需要估计的分位数,默认情况下是估计中位数。我们下面将实验组数据中的中位数处理效应和 ATE 进行比较。
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ate
. est sto ate1
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, qte
. qui estimates store qte1_50
. tesensitivity cpitable ate1 qte1_50
Treatment effects sensitivity
Analysis : conditional partial independence, multiple results
------------------------------------------------------------
c | ate1 | qte1_50
-------------------+--------------------+-------------------
0.000 | [ 1,419, 1,419] | [ 988, 988]
0.026 | [ 920, 1,930] | [ 321, 2,009]
0.051 | [ 433, 2,454] | [ -541, 2,857]
...
0.949 | [-11,594, 21,952] | [-60,308, 60,308]
0.974 | [-11,594, 21,952] | [-60,308, 60,308]
1.000 | [-11,594, 21,952] | [-60,308, 60,308]
------------------------------------------------------------
Analysis : breakdown point, multiple
-------------------------------------------
conclusion | ate1 qte1_50
------------------+------------------------
stat > 0 | 0.075 0.038
-------------------------------------------
. tesensitivity cpiplot ate1 qte1_50
对于两个条件统计量 CATE 和 CQTE 需要指定协变量的具体值,默认情况下将使用每一个协变量的均值,通过 median
选项可以使用每一个协变量的中位数。要是想要指定一些协变量的其他值,我们可以有两种方法:一是使用 covariates
选项,二是使用 qcovariates
选项来指定协变量的分位数,两者也可以同时使用 (在下述第三种情形之下)。
这两种方法均接受三种输入 (input):
如果使用第一种或者第二种,全部协变量的值都需要被提供,他们也必须和 tesensitivity
中的变量保持相同的顺序。如果使用第三种,则不需要对于全部的协变量都指定值,那些没有指定的变量将会默认采用均值,或者如果 median
选项被指定,这些没有指定的变量将会采用中位数。
例如,如果我们将要比较已婚人士在中位数收入下的 CATE 和未婚人士在 1974 年 90th 分位数收入下的 CATE,然后将其他未指定的控制变量设定为其中位数。
. matrix cov = (1)
. matrix colnames cov = married
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ///
> cate median cov(cov)
. qui estimates store cate1_50
. matrix cov[1,1] = 0
. matrix qcov = (.9)
. matrix colnames qcov = re74
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ///
> cate median qcov(qcov) cov(cov)
. qui estimates store cate1_90
. tesensitivity cpiplot cate1_50 cate1_90
我们也可以通过 e(covsupp)
的方式来查看协变量的具体设定。
. matrix list e(covsupp)
e(covsupp)[1,10]
married age black hispanic education re74 re75 re74pos re75pos _const
r1 0 24 1 0 10 7914.1309 0 0 0 1
. estimates restore cate1_9
. matrix list e(covsupp)
e(covsupp)[1,10]
married age black hispanic education re74 re75 re74pos re75pos _const
r1 0 24 1 0 10 7914.1309 0 0 0 1
我们也可以采用相同的方式对 50th 条件分位数处理效应展开分析。
. matrix cov = (1)
. matrix colnames cov = married
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ///
> cqte median cov(cov)
. qui estimates store cqte1_50
. matrix cov[1,1] = 0
. matrix qcov = (.9)
. matrix colnames qcov = re74
. qui tesensitivity cpi ($outcome $controls) ($treatment $controls) if sample1, ///
> cqte median qcov(qcov) cov(cov)
. qui estimates store cqte1_90
. tesensitivity cpiplot cqte1_50 cqte1_90
在之前的分析中,我们均假设结果变量是连续函数。对于二元结果变量而言,目前我们的命令只支持 ATE。我们考虑在 1978 年收入是否为正的结果变量。
. gen re78pos = re78 > 0
. tesensitivity cpi (re78pos $controls) ($treatment $controls) if sample1, ate
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 445
Outcome model : logistic Breakdown = 0.086
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78pos
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 0.1121, 0.1121]
0.026 | [ 0.0805, 0.1431]
0.051 | [ 0.0476, 0.1741]
...
0.949 | [-0.4787, 0.5213]
0.974 | [-0.4787, 0.5213]
1.000 | [-0.4787, 0.5213]
---------------------------------------
在进行上述 ate
、atet
、qte
等计算过程中可能会花费较长的时间,为了了解运算的进度,可以使用 verbose
选项。
. tesensitivity cpi (re78pos $controls) ($treatment $controls) if sample1, ate verbose
calculating ate (41)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.........................................
calculating breakdown point...
Treatment effects sensitivity
Analysis : cond. partial independence Number of obs = 445
Outcome model : logistic Breakdown = 0.086
Treatment model : logistic Conclusion = ate > 0
Outcome variable : re78pos
---------------------------------------
c | ate
-------------------+-------------------
0.000 | [ 0.1121, 0.1121]
0.026 | [ 0.0805, 0.1431]
0.051 | [ 0.0476, 0.1741]
...
---------------------------------------
在 tesensitivity cpiplot
绘图过程中,Stata 默认的 twoway_options
的选项都可以被传递到tesensitivity cpiplot
中。我们不妨对之前的 CQTE 比较图中的标题、背景颜色等做一些调整,以使得图形更加美观。
. qui tesensitivity cpiplot cqte1_50 cqte1_90, ///
> xtitle(c-dependence) ytitle(conditional quantile treatment effect) ///
> graphregion(color(%8) margin(vsmall)) ///
> title(CQTE at 90th and 50th percentiles of income)
我们也可以通过 boundpatterns
和 boundcolors
选项来对边界的形状和颜色做出调整。此外也可以将 graph twoway line
命令中的选项传递到 tesensitivity cpiplot
的选项 boundoptions
中。
. qui tesensitivity cpiplot cqte1_50 cqte1_90, ///
> xtitle(c-dependence) ytitle(conditional quantile treatment effect) ///
> graphregion(color(%8) margin(vsmall)) ///
> title(CQTE at 20th and 50th percentiles of income) ///
> boundcolors(navy ltblue) boundpatterns(solid) boundoptions(lwidth(vthick))
最后我们也可以对于中间那条截断点横线做出调整。例如,如果不想要这条线,就采用 nobreakdown
选项,也可以通过 breakdownoptions
的选项来控制其颜色形状等等,Stata 中 yline
的选项也都可以丢进这一选项中。
. qui tesensitivity cpiplot cqte1_50 cqte1_90, ///
> xtitle(c-dependence) ytitle(conditional quantile treatment effect) ///
> graphregion(color(%8) margin(vsmall)) ///
> title(CQTE at 20th and 50th percentiles of income) ///
> boundcolors(navy ltblue) boundpatterns(solid) boundoptions(lwidth(vthick)) ///
> breakdownoptions(lcolor(orange))
当然也可以去掉图例。
. qui tesensitivity cpiplot cqte1_50 cqte1_90, ///
> xtitle(c-dependence) ytitle(conditional quantile treatment effect) ///
> graphregion(color(%8) margin(vsmall)) ///
> title(CQTE at 20th and 50th percentiles of income) ///
> boundcolors(navy ltblue) boundpatterns(solid) boundoptions(lwidth(vthick)) ///
> breakdownoptions(lcolor(orange)) nolegend
也可以调整图例的形式。
. qui tesensitivity cpiplot cqte1_50 cqte1_90, ///
> xtitle(c-dependence) ytitle(conditional quantile treatment effect) ///
> graphregion(color(%8) margin(vsmall)) ///
> title(CQTE at 20th and 50th percentiles of income) ///
> boundcolors(navy ltblue) boundpatterns(solid) boundoptions(lwidth(vthick)) ///
> breakdownoptions(lcolor(orange)) legoptions(region(fcolor(%8)))
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