Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:符舜 (中山大学)
邮箱:fush23@mail2.sysu.edu.cn
编者按:本文主要摘译自下文,特此致谢!
Source:Katsouris C. Treatment effect validation via a permutation test in Stata[J]. Available at SSRN 3830527, 2020. -PDF-
目录
在经济学研究中,学者经常希望评估某个项目或政策实施后的效果,比如某地政府推出的禁烟政策对当地香烟消费是否有抑制作用,上述政策效应通常也被称为处理效应 (treatment effect)。然而在实际情况中,我们可能会遇到基线不均衡 (baseline imbalance) 和非随机性样本缺失 (non-random attrition) 的问题,因此需要对结果进行检验和修正。
本文介绍的新命令 permtest
旨在通过排序检验 (permutation test) 来对处理效应进行非参数检验,并通过相应的方法来修正基线不均衡和非随机性样本缺失问题。
处理效应是基于 Rubin 于 1974 年提出的 “反事实框架” 下得出的。以虚拟变量
其中,
对于政策制定者而言,可能更为关注
为了解决选择偏差,我们可以采用随机分组 (random assignment) 的方法,使得
然而,现实情况下随机分组并非在所有情况下都可行 (可能成本很高)。
在多数情况下,我们可以假设个体是否参加项目由一些可测变量决定,如年龄、性别等,称这些变量称为协变量。如果个体对
此时,我们就可以得到处置效应的线性回归模型,其中
即使如此,在现实情况下,我们并不知道有多少个协变量,也不确定方程是否为线性形式,以及可能会遇到基线不均衡和样本缺失问题,导致仍然会出现偏差。
排序检验,又称置换检验,最早是由 R.A.Fisher 于 1953 年提出,目前是在因果推断中常用的一种非参数检验方法,特别是对于一些样本容量小且分布未知的情况。相比于传统参数检验,置换检验并不需要正态性和同方差假定,仅依赖于样本的随机排列。
排序检验的基本思想和 Bootstrap 方法类似,都是通过重复抽样的方法来近似得到统计量的精确分布。不同点在于 Bootstrap 方法是有放回抽样,而排序检验是不放回抽样,对样本进行排列组合。下面简要介绍排序检验的基本思想。
假定
如果我们考虑所有可能的排列组合,此时的置换检验称作 “精确” 检验。随着样本量的增加,获取所有可能排列的计算量会非常大。这种情况下,我们可以从所有可能的排列中进行抽样,获得一个近似的检验。
基线不均衡问题通常表现为协变量在处理组和对照组之间不均衡,这会造成混淆 (confounding),违背了
如下图所示,我们以一个母亲是否吸烟和出生婴儿体重的例子来说明基线不均衡问题。红点表示怀孕期间吸烟的母亲,绿点表示怀孕期间未吸烟的母亲,母亲年龄在是否吸烟这一对照实验之间存在不均衡,因此有可能是吸烟导致出生婴儿体重偏轻,也有可能是母亲年龄较大导致出生婴儿体重偏轻。因此直接将吸烟与不吸烟母亲的出生婴儿体重作比较会导致结果存在选择偏差。
为了对基线不均衡问题进行修正,Christis 和 Katsouris (2020) 采用了 Freedman 和 Lane (1983) 提出的对残差进行置换检验的方法。考虑线性回归方程
其背后原理也非常简单,通过对
在随机干预实验实施过程中可能会出现样本中途退出或未能追踪到结果等情况,即样本缺失问题,这会导致实验估计的处理效应存在偏差,甚至无效。
一方面,可能流失的样本与最初样本有显著差异。比如评估一禁烟政策对香烟消费量的影响时,平均香烟消费量最高的一批受试者的数据缺失,则评估结果可能会低估该政策的影响。
另一方面,也可能流失的样本在处理组和控制组之间有显著差异。比如在处理组中缺失的样本普遍平均香烟消费量很高,而在控制组中缺失的样本普遍平均香烟消费量很低,此时干预组与对照组不再具有可比性,评估结果也将会存在偏差。
为了对样本缺失问题进行修正,Christis 和 Katsouris (2020) 采用了 Wooldridge (2010) 提出的逆概率加权法 (IPW,inverse probability weights)。下面介绍 IPW 的基本思想:
设向量
对于
令
* 命令安装
lxhget permtest.pkg, install replace
* 命令语法
permtest [{outcomes}] {if}, treat(varlist) np(#) [options]
其中,
outcomes
:实验的结果变量;treat
:实验的处理变量,必须为为二元变量,用于区分处理组和对照组;np
:进行排序检验的次数。options
选项如下:
blockvars(varlist)
:分组进行排列组合,括号内变量为分组依据变量;naive
:跳过 blockvars
选项,直接对全样本进行操作;lcvars(varlist)
:添加协变量;ipwcovars1-5(varlist)
:IPW 方法中用于预测缺失率的变量,可以是 lcvars
中的变量,最多能添加 5 个变量;link(string)
:IPW 方法选择 Probit 或 Logit 模型来预测概率,但目前仅支持 Logit 模型;reverse
:对结果变量取相反数,计算相应的处理效应;effsize
:报告效应量;savemat(path)
:将矩阵结果以 csv 文件的格式输出。
在本节,我们将基于 permtest
命令,通过蒙特卡洛模拟来探究在不同的情境下排序检验和 simulate
命令来进行模拟。由于模拟次数过少,模拟结果与原文存在较大出入,本文旨在展示论文的思路和结论。
先对蒙特卡洛模拟的基本参数进行一个简要说明:
*=================蒙特卡洛模拟===============
// replications (rp):50 times 注:原论文进行了1000次模拟
// number of permutations (np):50 times 注:原论文进行了1000次排序检验
// treatment indicator variable: d = (rnormal(0,1) <0)
// number of observations: obs = {10,15,20}
// treatment effect: gamma = {0,0.1,0.2 ... ,0.5}
// Distribution of errors: {normal、Student、Weibull、Log-normal、Cauchy}
// H0:gamma = 0; H1:gamma ≠ 0
// 以 P值 = 0.05为标准计算拒绝率(小于记为1,大于记为0,求其平均值)
// (与原文有所不同,原文是检验 H0:gamma=0 和 H0:gamma > 0 的empirical size and power,
// 本文统一计算在H0: gamma=0 下的拒绝率)
// significance level: α = 0.05
// Empirical size and power: reject rate = #(p < α)/rp (输出H0的拒绝率)
// gamma = 0 时拒绝率越低证明效果越好,gamma ≠ 0 时拒绝率越高证明效果越好
*=============================================
本节总共讨论五种模型,在开始模拟前,本人自编了命令 permutations
用于数据生成 (DGP) 和排序检验处理,下文会详细说明每个模型的具体设定。
********************************************************************************
// permutaions - 自编进行DGP和排序检验过程的程序
********************************************************************************
// Options:
// obs() - 观察值个数
// distri() - 残差的分布,取值为1、2、3、4、5,
// 分别为正态分布、t分布、威尔布尔分布、对数正态分布、柯西分布
// gamm() - 处理变量的系数
// cov() - 是否存在协变量,0为不存在,1为存在
// unbal() - 是否存在变量基线不均衡问题及问题类型,
// 0为不存在,1为协变量基线不均衡,2为处理变量基线不均衡
// attri() - 是否存在样本流失问题,0为不存在,1为存在
********************************************************************************
cap program drop permutations
*-------------------------------------------------------------------------------
program define permutations, rclass
version 16
syntax ///
[, ///
obs(integer 20) ///
distri(integer 1) ///
gamm(real 0.1) ///
cov(integer 0) ///
unbal(integer 0) ///
attri(integer 0) ///
]
drop _all
set obs `obs'
gen d = (rnormal(0,1) <0) //生成treatment variable
//对应不同的分布生成残差
if `distri'==1 {
gen e_1 = rnormal(0,1)
}
if `distri' ==2{
gen e_2 = rt(`obs'-1)
}
if `distri' ==3{
gen e_3 = rweibull(1,1.5)
}
if `distri' ==4{
gen e_4 = exp(rnormal(0,1))
}
if `distri' ==5{
gen e_5 = rcauchy(0,1)
}
*1存在协变量情况
if `cov' != 0{
//生成协变量
gen x1 = rnormal(0,2)
gen x2 = round(runiform()) //生成随机0-1变量,可理解成性别
*1-1存在协变量不均衡
if `unbal' == 1{
replace x1 = rnormal(0.25,2) if d==1
replace x1 = rnormal(-0.25,2) if d==0
}
*1-2存在处理不均衡
if `unbal' ==2{
replace d = rbinomial(1,0.3)
}
//对应不同的gamma系数生成对应的outcome variable
gen y = `gamm'*d + 0.2*x1 + 0.05*x2 + e_`distri'
permtest y, treat(d) lcvars(x1 x2) np(50)
local pval_reg1s (r(pval_reg1s)[1,1] < 0.05)
local pval_perm (r(pval_perm)[1,1] < 0.05)
return scalar para_reject =`pval_reg1s' //t检验拒绝情况
return scalar nonpara_reject =`pval_perm' //置换检验拒绝情况
}
*2不存在协变量情况
if `cov' == 0 & `attri' == 0{
gen y = `gamm'*d + e_`distri'
permtest y, treat(d) np(50)
local pval_reg1s (r(pval_reg1s)[1,1] < 0.05)
local pval_perm (r(pval_perm)[1,1] < 0.05)
return scalar para_reject =`pval_reg1s' //t检验拒绝情况
return scalar nonpara_reject =`pval_perm' //置换检验拒绝情况
}
*3存在样本缺失情况
if `cov' == 0 & `attri' !=0{
gen v = rnormal(0,1)
gen x1 = rnormal(0,2)
gen y = `gamm'*d + e_`distri'
gen R = (0.1+0.25*d+3.75*x1+v > 0)
replace y=. if R==0
permtest y, treat(d) np(50) ipwcovars1(x1) //ipw方法采用logit模型
local pval_reg1s_ipw (r(pval_reg1s_ipw)[1,1] < 0.05)
local pval_perm_ipw (r(pval_perm_ipw)[1,1] < 0.05)
local pval_reg1s (r(pval_reg1s)[1,1] < 0.05)
local pval_perm (r(pval_perm)[1,1] < 0.05)
return scalar para_reject =`pval_reg1s' //没有逆概率加权的t检验拒绝率
return scalar nonpara_reject =`pval_perm' //没有逆概率加权的置换检验拒绝率
return scalar para_reject_ipw = `pval_reg1s_ipw' //存在逆概率加权的t检验拒绝率
return scalar nonpara_reject_ipw = `pval_perm_ipw' //没有逆概率加权的置换检验拒绝率
}
end
*----------Note: 执行后续命令前,请先选中上述程序,按快捷键 Ctrl+R--------------
通过模型 1 我们可以比较排序检验和
*1、baseline experimental design
* 1-1、Model 1 (Treatment Effect estimation with no covariates)
forvalues i = 20(20)60{
forvalues j = 1/5{
forvalues k = 0(0.1)0.5{
local k_int = ceil(`k'*10)
simulate model1_para_reject_`i'_`j'_`k_int'=r(para_reject) ///
model1_nonpara_reject_`i'_`j'_`k_int'=r(nonpara_reject), ///
reps(50) seed(12345): permutations, obs(`i') distri(`j') gamm(`k')
gen id=_n
save "model1_`i'_`j'_`k_int'.dta" , replace
clear
}
}
}
mergemany 1:1 all, match(id) all // 将所有的结果 merge
logout, save("./result1") excel replace: tabstat model1*, ///
stat(mean) format(%7.3f) varwidth(30) columns(statistic) // 平均值 mean 即为我们求的拒绝率
从表 1 中可以看到,残差服从正态分布时,排序检验的效果和传统
* 1-2、Model 2 (Treatment Effect estimation with covariates)
clear
forvalues i = 20(20)60{
forvalues j = 1/5{
forvalues k = 0(0.1)0.5{
local k_int = ceil(`k'*10)
simulate model2_para_reject_`i'_`j'_`k_int'=r(para_reject) ///
model2_nonpara_reject_`i'_`j'_`k_int'=r(nonpara_reject), ///
reps(50) seed(12345): permutations, obs(`i') distri(`j') gamm(`k') cov(1)
gen id=_n
save "model2_`i'_`j'_`k_int'.dta" , replace
clear
}
}
}
mergemany 1:1 all, match(id) all // 将所有的结果 merge
logout, save("./result") excel replace: tabstat model2*, ///
stat(mean) format(%7.3f) varwidth(30) columns(statistic) // 平均值 mean 即为表格中的拒绝率
模型 2 在模型 1 的基础上加入了一个连续变量
这一节我们考虑了两种基线不均衡的情况,一是是协变量
我们的预期在原样本数据存在混淆和偏差的情况下,排序检验能够比
*2、experimental design with unbalanced covariates
* 2-1、Model 3 (Treatment Effect estimation with Unbalanced continuous covariate)
clear
forvalues i = 20(20)60{
forvalues j = 1/5{
forvalues k = 0(0.1)0.5{
local k_int = ceil(`k'*10)
simulate model3_para_reject_`i'_`j'_`k_int'=r(para_reject) ///
model3_nonpara_reject_`i'_`j'_`k_int'=r(nonpara_reject), ///
reps(50) seed(12345): permutations, obs(`i') distri(`j') gamm(`k') cov(1) unbal(1)
gen id=_n
save "model3_`i'_`j'_`k_int'.dta", replace
clear
}
}
}
mergemany 1:1 all, match(id) all // 将所有的结果 merge
logout, save("./result") excel replace: tabstat model3*, ///
stat(mean) format(%7.3f) varwidth(30) columns(statistic) // 平均值 mean 即为表格中的拒绝率
* 2-2、Model 4 (Treatment Effect estimation with Unbalanced treatment)
clear
forvalues i = 20(20)60{
forvalues j = 1/5{
forvalues k = 0(0.1)0.5{
local k_int = ceil(`k'*10)
simulate model4_para_reject_`i'_`j'_`k_int'=r(para_reject) ///
model4_nonpara_reject_`i'_`j'_`k_int'=r(nonpara_reject), ///
reps(50) seed(12345): permutations, obs(`i') distri(`j') gamm(`k') cov(1) unbal(2)
gen id=_n
save "model4_`i'_`j'_`k_int'.dta", replace
clear
}
}
}
mergemany 1:1 all, match(id) all // 将所有的结果 merge
logout, save("./result") excel replace: tabstat model4*, ///
stat(mean) format(%7.3f) varwidth(30) columns(statistic) //平均值 mean 即为表格中的拒绝率
根据表 3 和表 4 的结果,再对比表 2,可以发现排序检验和
这一节我们考虑了样本缺失的情况,样本缺失问题可能会使处理组和对照组不再具有代表性,带来估计结果的偏差并降低我们识别出显著的处理效应的可能性。
样本缺失问题可分为随机性和非随机性的,基于其重复抽样的性质,排序检验可以较好地修正随机性样本缺失问题。但是当缺失问题是由于实验个体的一些特点所决定即非随机性的时候,排序检验不再有效,我们有必要采用 IPW 的方法对其进行修正以减小偏差。 Christis 和 Katsouris (2020) 假定样本缺失问题是基于可观测变量引起的,借鉴 Huber (2012) 的实验设计,考虑了如下模型:
Model 5:Attrition
*3、experimental design with attrition
clear
forvalues i = 20(20)60{
forvalues j = 1/5{
forvalues k = 0(0.1)0.5{
local k_int = ceil(`k'*10)
simulate model5_para_reject_`i'_`j'_`k_int'=r(para_reject) ///
model5_nonpara_reject_`i'_`j'_`k_int'=r(nonpara_reject) ///
model5_para_reject_ipw_`i'_`j'_`k_int'=r(para_reject_ipw) ///
model5_nonpara_reject_ipw_`i'_`j'_`k_int'=r(nonpara_reject_ipw), ///
reps(50) seed(12345): permutations, obs(`i') distri(`j') gamm(`k') attri(1)
gen id=_n
save "model5_`i'_`j'_`k_int'.dta", replace
clear
}
}
}
mergemany 1:1 all, match(id) all // 将所有的结果 merge
logout, save("./result") excel replace: tabstat model5*, ///
stat(mean) format(%7.3f) varwidth(30) columns(statistic) // 平均值 mean 即为表格中的拒绝率
从表 5 和表 6 的结果可以看出,通过比较使用和不使用 IPW 方法的拒绝率,我们发现不使用 IPW 方法时,
通过比较排序检验和
理论上,在满足随机分组或者可忽略性的假设下,我们能够通过线性回归模型得到无偏的处理效应,但是现实情况中经常会出现基线不均衡问题和样本缺失问题,使得估计结果出现偏差。
permtest
命令旨在通过运用排序检验和 IPW 的方法来修正基线不均衡和样本缺失问题,更稳健地识别出处理效应的存在。
基于蒙特卡洛模拟,我们发现排序检验相比于
Note:产生如下推文列表的 Stata 命令为:
lianxh 模拟, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh