Stata：处理效应的排序检验-permtest

发布时间：2022-01-25 阅读 2507

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下载 - 推文合集

作者：符舜 (中山大学)
邮箱：fush23@mail2.sysu.edu.cn

编者按：本文主要摘译自下文，特此致谢！
Source：Katsouris C. Treatment effect validation via a permutation test in Stata[J]. Available at SSRN 3830527, 2020. -PDF-

1. 理论背景

在经济学研究中，学者经常希望评估某个项目或政策实施后的效果，比如某地政府推出的禁烟政策对当地香烟消费是否有抑制作用，上述政策效应通常也被称为处理效应 (treatment effect)。然而在实际情况中，我们可能会遇到基线不均衡 (baseline imbalance) 和非随机性样本缺失 (non-random attrition) 的问题，因此需要对结果进行检验和修正。

本文介绍的新命令 permtest 旨在通过排序检验 (permutation test) 来对处理效应进行非参数检验，并通过相应的方法来修正基线不均衡和非随机性样本缺失问题。

1.1 处理效应

处理效应是基于 Rubin 于 1974 年提出的 “反事实框架” 下得出的。以虚拟变量 $D_{i} = {0, 1}$ 表示个体 $i$ 是否参加项目，称 $D_{i}$ 为处理变量 (treatment variable)。记我们感兴趣的结果变量为 $y$ ，对于每个个体 $y$ 有两种状态，即：

y = {\begin{cases} y_{1 i}, d = 1 \\ y_{0 i}, d = 0 \end{cases}

其中， $(y_{1 i} - y_{0 i})$ 为个体 $i$ 参加项目的处理效应，显然现实中我们不可能同时观测到 $y_{0 i}$ 和 $y_{1 i}$ 。由于处理效应 $(y_{1 i} - y_{0 i})$ 为随机变量，我们关心其期望值，分别定义 $A T E$ 、 $A T T$ 、 $A T U$ 为平均处理效应、参与者平均处理效应、非参与者平均处理效应，即

\begin{aligned} A T E & = E {y_{1 i} - y_{0 i}} \\ A T T & = E {y_{1 i} - y_{0 i} | D = 1} \\ A T U & = E {y_{1 i} - y_{0 i} | D = 0} \end{aligned}

对于政策制定者而言，可能更为关注 $A T T$ ，因为它衡量的是对项目参与者的效果。但是由于 $y_{0 i}$ 和 $y_{1 i}$ 不能同时得到，如果我们直接根据观测到的变量进行比较会导致选择偏差：

\begin{aligned} E (y_{1 i} | D = 1) - E (y_{0 i} | D = 0) & = \underset{A T E}{\underset{⏟}{E (y_{1 i} | D = 1) - E (y_{0 i} | D = 1)}} \\ + \underset{选择偏差}{\underset{⏟}{E (y_{0 i} | D = 1) - E (y_{0 i} | D = 0)}} \end{aligned}

为了解决选择偏差，我们可以采用随机分组 (random assignment) 的方法，使得 $D_{i}$ 独立于 $(y_{0 i}, y_{1 i})$ ，这样的实验方式称为随机干预实验。此时我们有 $A T E = A T T = A T U$ ，这就消除了选择偏差的影响。

y_{i} = γ D_{i} + ϵ_{i}

然而，现实情况下随机分组并非在所有情况下都可行 (可能成本很高)。

在多数情况下，我们可以假设个体是否参加项目由一些可测变量决定，如年龄、性别等，称这些变量称为协变量。如果个体对 $D_{i}$ 的选择完全取决于 $x_{i}$ ，则在给定 $x_{i}$ 的情况下，潜在结果 $(y_{0 i}, y_{1 i})$ 将独立于 $D_{i}$ ，这就是 Rosenbaum 和 Rubin (1983) 所引入的 “可忽略性” 假设 ( $Y ⊥ ⊥ D | X$ )：

\begin{aligned} E (y_{0 i} | x_{i}, D i) & = E (y_{0 i} | x_{i}) \\ E (y_{1 i} | x_{i}, D i) & = E (y_{1 i} | x_{i}) \end{aligned}

此时，我们就可以得到处置效应的线性回归模型，其中 $γ$ 的估计系数代表平均处理效应：

y_{i} = γ D_{i} + β X_{i} + ϵ_{i}

即使如此，在现实情况下，我们并不知道有多少个协变量，也不确定方程是否为线性形式，以及可能会遇到基线不均衡和样本缺失问题，导致仍然会出现偏差。

1.2 排序检验

排序检验，又称置换检验，最早是由 R.A.Fisher 于 1953 年提出，目前是在因果推断中常用的一种非参数检验方法，特别是对于一些样本容量小且分布未知的情况。相比于传统参数检验，置换检验并不需要正态性和同方差假定，仅依赖于样本的随机排列。

排序检验的基本思想和 Bootstrap 方法类似，都是通过重复抽样的方法来近似得到统计量的精确分布。不同点在于 Bootstrap 方法是有放回抽样，而排序检验是不放回抽样，对样本进行排列组合。下面简要介绍排序检验的基本思想。

假定 $(X_{1}, . . ., X_{n})$ 是来自同一总体的样本， $T_{n} (X_{1}, . . ., X_{n})$ 是我们检验的统计量。通过对样本进行重复排列组合，我们可得到复制的统计量 $T_{n}^{*} (X_{1}^{*}, . . ., X_{n}^{*})$ ，重复这一过程，总共有 $n!$ 个全排列，可得到 $n!$ 个 $T_{n}^{*}$ 。由于 $(X_{1}, . . ., X_{n})$ 来自同一分布，故每一个 $T_{n}^{*}$ 的取值可能相等，通过下式我们即可得到 $T_{n}$ 的 $p$ 值：

\hat{p} = \hat{P r} (T^{*} > T) = \frac{# (T^{r *} > T)}{n!}

如果我们考虑所有可能的排列组合，此时的置换检验称作 “精确” 检验。随着样本量的增加，获取所有可能排列的计算量会非常大。这种情况下，我们可以从所有可能的排列中进行抽样，获得一个近似的检验。

1.3 基线不均衡问题

基线不均衡问题通常表现为协变量在处理组和对照组之间不均衡，这会造成混淆 (confounding)，违背了 $y_{i}$ 独立于 $D_{i}$ 的假设，导致我们的估计是有偏的。

如下图所示，我们以一个母亲是否吸烟和出生婴儿体重的例子来说明基线不均衡问题。红点表示怀孕期间吸烟的母亲，绿点表示怀孕期间未吸烟的母亲，母亲年龄在是否吸烟这一对照实验之间存在不均衡，因此有可能是吸烟导致出生婴儿体重偏轻，也有可能是母亲年龄较大导致出生婴儿体重偏轻。因此直接将吸烟与不吸烟母亲的出生婴儿体重作比较会导致结果存在选择偏差。

为了对基线不均衡问题进行修正，Christis 和 Katsouris (2020) 采用了 Freedman 和 Lane (1983) 提出的对残差进行置换检验的方法。考虑线性回归方程 $y_{i} = γ D_{i} + β X_{i} + ϵ_{i}$ ，现对处理效应即系数 $γ$ 进行检验 $H_{0} : γ = 0$ ，具体过程如下：

第一步，通过最小二乘法对参数 $γ$ 进行估计，求得其 $t$ 统计量 $t = \frac{\hat{γ}}{s e (\hat{γ})}$ ；
第二步，将 $y_{i}$ 对 $X_{i}$ 进行带截距项回归，得到残差 ${\hat{v}}_{i} = y_{i} - {\hat{c}}_{0} - \hat{β} X_{i}$ ；
第三步，对残差 ${\hat{v}}_{i}$ 进行随机置换得到置换的误差 $v_{i}^{*}$ ，令 $y_{i}^{*} = {\hat{c}}_{0} + \hat{β} X_{i} + v_{i}^{*}$ ；
第四步，基于 $(y_{i}^{*}, X_{i}, D_{i})$ 对原方程进行拟合，得到 ${\hat{γ}}_{*}$ 及其置换后的 $t$ 统计量 $T^{*}$ ；
第五步，重复步骤 2-3 $n p$ 次，得到相应的复制统计量 $T_{1}^{*}, . . ., T_{n p}^{*}$ ，求得 $p$ 值的近似值为 $\hat{p} = \hat{P r} (T^{*} > T) = \frac{# (T^{*} > T)}{n p}$ 。

其背后原理也非常简单，通过对 $v_{i}^{*}$ 的随机排列组合，我们能够实现 $v ⊥ ⊥ D | X$ ，从而得到 $Y ⊥ ⊥ D | X$ 。

1.4 样本缺失问题

在随机干预实验实施过程中可能会出现样本中途退出或未能追踪到结果等情况，即样本缺失问题，这会导致实验估计的处理效应存在偏差，甚至无效。

一方面，可能流失的样本与最初样本有显著差异。比如评估一禁烟政策对香烟消费量的影响时，平均香烟消费量最高的一批受试者的数据缺失，则评估结果可能会低估该政策的影响。

另一方面，也可能流失的样本在处理组和控制组之间有显著差异。比如在处理组中缺失的样本普遍平均香烟消费量很高，而在控制组中缺失的样本普遍平均香烟消费量很低，此时干预组与对照组不再具有可比性，评估结果也将会存在偏差。

为了对样本缺失问题进行修正，Christis 和 Katsouris (2020) 采用了 Wooldridge (2010) 提出的逆概率加权法 (IPW，inverse probability weights)。下面介绍 IPW 的基本思想：

设向量 $(w_{1 i}, . . ., w_{n i}, s_{i})$ 为样本观察变量，其中 $s_{i}$ 为二元变量，我们有以下假设：

当 $s_{i} = 1$ 时， $w_{i}$ 可被观测到，当 $s_{i} = 0$ 时， $w_{i}$ 不能或仅能部分被观测到；
存在随机向量 $z_{i} = (z_{1 i}, . . ., z_{J i})$ 使得 $P (s_{i} = 1 | z_{i}, w_{i}) = P (s_{i} = 1 | z_{i}) \equiv p (z_{i})$ ；
对所有 $z \in ζ \subset R^{J}$ ，均有 $p (z) > 0$ ；
不论 $s_{i}$ 是否为零， $z_{i}$ 总能够被观测到。

对于 $p (z_{i})$ 的计算，我们可以采用 Probit 或 Logit 模型来进行估计，以其倒数 $1 / p (z_{i})$ 为权重进行加权运算，就能够解决样本缺失带来的结果偏差，证明过程如下：

令 $g (w)$ 为任意使得 $E [g (w)]$ 存在的数学函数，则有：

\begin{aligned} E [s_{i} g (w) / p (z_{i})] & = E {E [s_{i} g (w) / p (z_{i}) | w_{i}, z_{i}]} \\ = E {E [s_{i} | w_{i}, z_{i}] g (w) / p (z_{i})} \\ = E {P (s_{i} = 1 | w_{i}, z_{i}) g (w) / p (z_{i})} \\ = E {p (z_{i}) g (w) / p (z_{i})} \\ = E [g (w)] \end{aligned}

2. 命令介绍

* 命令安装
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 文件的格式输出。

3. 蒙特卡洛模拟

在本节，我们将基于 permtest 命令，通过蒙特卡洛模拟来探究在不同的情境下排序检验和 $t$ 检验的表现。本文使用 Stata 的 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--------------

3.1 基础情况

3.1.1 Model 1：Treatment Effect estimation with no covariates

通过模型 1 我们可以比较排序检验和 $t$ 检验在不同残差分布下的有效性，其中 $t$ 分布和威尔布尔分布可以看作偏态分布，柯西分布代表厚尾分布，而对数正态分布视为偏态分布和厚尾分布的结合。

y_{j} = γ D_{j} + ϵ^{(k)}, ϵ^{(k)} \sim F (ϵ^{(k)}) j = {1, . . ., n}

ϵ^{(1)} \sim N (0, 1) ϵ^{(2)} \sim t (n - 1) ϵ^{(3)} \sim W (a = 1, b = 1.5)

ϵ^{(4)} \sim \exp (N (0, 1)) ϵ^{(5)} \sim C (a = 1, b = 0.25)

*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 中可以看到，残差服从正态分布时，排序检验的效果和传统 $t$ 检验基本是一样好的，这有力地证明了排序检验的有效性。在样本容量小且残差偏度越大时，排序检验的优势就更为明显。特别是当样本量为 20，残差服从柯西分布时，虽然相比于正态分布条件，两种检验的拒绝率均降低，但是 $t$ 检验出现了严重的偏差，而排序检验的偏差相对较小。

3.1.2 Model 2：Treatment Effect estimation with covariates

y_{j} = γ D_{j} + 0.2 * X_{1} + 0.05 * X_{2} + ϵ^{(k)}, ϵ^{(k)} \sim F (ϵ^{(k)}) j = {1, . . ., n}

X_{1} \sim N (0, 2) X_{2} \sim r o u n d (U (0, 1))

ϵ^{(1)} \sim N (0, 1) ϵ^{(2)} \sim t (n - 1) ϵ^{(3)} \sim W (a = 1, b = 1.5)

ϵ^{(4)} \sim \exp (N (0, 1)) ϵ^{(5)} \sim C (a = 1, b = 0.25)

* 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 的基础上加入了一个连续变量 $X_{1}$ 和一个二元变量 $X_{2}$ 。根据 Nagelkerke 等 (2000)，在估计中加入协变量能够修正混淆效应，从而识别出处理效应是否真实存在。从表 2 中可以明显的看到，在残差服从柯西分布，样本容量为 20 的条件下，加入协变量使得 $t$ 检验和排序检验的效果相对表 1 有了显著改善。

3.2 存在基线不均衡的情况

这一节我们考虑了两种基线不均衡的情况，一是是协变量 $X_{1}$ 的均值在处理组和对照组之间有差异；二是接受处理的个体数并非占总个数的一半 ( $D \sim b i n o m i a l (1, 0.3)$ )，正如前面所说，这会带来混淆，导致估计结果是有偏的。

我们的预期在原样本数据存在混淆和偏差的情况下，排序检验能够比 $t$ 检验更稳健地识别出处理效应。

3.2.1 Model 3：Unbalanced continuous covariate

y_{j} = γ D_{j} + 0.2 * X_{1} + 0.05 * X_{2} + ϵ^{(k)}, ϵ^{(k)} \sim F (ϵ^{(k)}) j = {1, . . ., n}

X_{1} \sim N (0.25, 2) if d = 1 and X_{1} \sim N (- 0.25, 2) if d = 0

ϵ^{(1)} \sim N (0, 1) ϵ^{(2)} \sim t (n - 1) ϵ^{(3)} \sim W (a = 1, b = 1.5)

ϵ^{(4)} \sim \exp (N (0, 1)) ϵ^{(5)} \sim C (a = 1, b = 0.25)

*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 即为表格中的拒绝率

3.2.2 Model 4：Unbalanced treatment

y_{j} = γ D_{j} + 0.2 * X_{1} + 0.05 * X_{2} + ϵ^{(k)}, ϵ^{(k)} \sim F (ϵ^{(k)}) j = {1, . . ., n}

D_{j} \sim B i n o m i a l (1, 0.3)

ϵ^{(1)} \sim N (0, 1) ϵ^{(2)} \sim t (n - 1) ϵ^{(3)} \sim W (a = 1, b = 1.5)

ϵ^{(4)} \sim e x p (N (0, 1)) ϵ^{(5)} \sim C (a = 1, b = 0.25)

* 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，可以发现排序检验和 $t$ 检验在 $γ$ 较大时的拒绝率都有了明显的下降，这代表了基线不均衡所带来的偏差问题。另外，在某些情况下，排序检验相较于 $t$ 检验也确实能够改善基线不均衡问题，比如在 $n = 60, γ = 0.5$ 的情况下，排序检验在五种分布下的拒绝率均大于等于 $t$ 检验。

3.3 存在样本缺失的情况

这一节我们考虑了样本缺失的情况，样本缺失问题可能会使处理组和对照组不再具有代表性，带来估计结果的偏差并降低我们识别出显著的处理效应的可能性。

样本缺失问题可分为随机性和非随机性的，基于其重复抽样的性质，排序检验可以较好地修正随机性样本缺失问题。但是当缺失问题是由于实验个体的一些特点所决定即非随机性的时候，排序检验不再有效，我们有必要采用 IPW 的方法对其进行修正以减小偏差。 Christis 和 Katsouris (2020) 假定样本缺失问题是基于可观测变量引起的，借鉴 Huber (2012) 的实验设计，考虑了如下模型：

Model 5：Attrition

y_{j} = γ D_{j} + ϵ^{(k)}, ϵ^{(k)} \sim F (ϵ^{(k)}) j = {1, . . ., n}

R_{j} = 1 {0.1 + 0.25 D_{j} + 3.75 * X_{1} + v > 0}

v \sim N (0, 1) X_{1} \sim N (0, 2), y_{j} = . if R_{j} = 0

ϵ^{(1)} \sim N (0, 1) ϵ^{(2)} \sim t (n - 1) ϵ^{(3)} \sim W (a = 1, b = 1.5)

ϵ^{(4)} \sim e x p (N (0, 1)) ϵ^{(5)} \sim C (a = 1, b = 0.25)

*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 方法时， $γ = 0$ 的拒绝率偏高而 $γ \neq 0$ 的拒绝率偏低，这表明样本缺失问题对估计结果造成的偏差较大，需要加以修正。

通过比较排序检验和 $t$ 检验的拒绝率，我们发现排序检验在 IPW 方法下的改进更为明显，特别是当残差服从柯西分布时，排序检验在 IPW 方法下的表现要优于在 $t$ 检验在 IPW 方法下的表现。

4. 总结

理论上，在满足随机分组或者可忽略性的假设下，我们能够通过线性回归模型得到无偏的处理效应，但是现实情况中经常会出现基线不均衡问题和样本缺失问题，使得估计结果出现偏差。

permtest 命令旨在通过运用排序检验和 IPW 的方法来修正基线不均衡和样本缺失问题，更稳健地识别出处理效应的存在。

基于蒙特卡洛模拟，我们发现排序检验相比于 $t$ 检验有如下优点：

在正态分布假定下，排序检验效果和 $t$ 检验至少一样好；当样本容量较小且残差分布的偏度较大时，排序检验要明显优于 $t$ 检验；
排序检验能够在一定程度上修正基线不均衡问题；
在使用 IPW 方法修正样本缺失问题时，排序检验在残差分布的偏度较大时效果更好。

5. 参考资料

Christis, Katsouris, Treatment Effect Validation via a Permutation Test in Stata, ArXiv Preprint ArXiv:2110.12268, 2020. -PDF-
Wooldridge J. M, Econometric analysis of cross section and panel data, MIT press, 2010: 821-827.
David Freedman, David Lane, A Nonstochastic Interpretation of Reported Significance Levels, Journal of Business & Economic Statistics, 1:4, 292-298, 1983. -PDF-
陈强. 高级计量经济学及 Stata 应用 (第二版)[M]. 高等教育出版社, 2014.
STATA中的治疗效果：RA：回归调整、 IPW：逆概率加权、 IPWRA、 AIPW -Link-

6. 相关推文

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

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]