当PSM遇上RDD:rddsga命令详解

发布时间:2021-06-18 阅读 795

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

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

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh

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

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

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:张子楠 (浙江财经大学)
邮箱zinanzh@gmail.com


目录


1. 简介

1.1 组间系数差异比较的两个问题及相应对策

在实证分析中,一种常见的操作是通过分组回归来比较组间系数差异,来识别子样本组因果效应的差异。例如不同所有制企业,不同地区、不同种族等。但这种直接分组回归的思路存在两个问题:

问题一: 不同子样本组回归系数虽然都显著,但置信区间存在重叠,从而难以通过系数相对大小来判断谁更经济显著。针对该问题,一种常见的解决思路是不分组回归,而是直接在全样本回归中加入以及分组指示变量(G)与核心解释变量(Z)的交叉项(G*Z)。如果交叉项系数显著,则说明不同子样本组间存在显著差异。

问题二: 不同子样本在协变量上可能存在显著差异。问题一中加入交叉项的思路,假定了协变量组间不存在显著差异。针对该问题,则可以通过继续增加分组指示变量与所有协变量(X)的交叉项来缓解,即增加 G,G*Z,以及 G*X*Z 三个选项。

注:这部分的详细分析可参考连享会推文:Stata: 如何检验分组回归后的组间系数差异?

RDD 回归组间系数比较也存在类似的问题,Geradino etc (2017) 认为如果只是简单增加 G*Z,则有可能出现识别错误。在他们文章附录 C 给出的例子里,不同子样本组的真实断点回归结果都不显著,而如果采用简单增加 G*Z 交叉项,则回归结果会错误报告为子样本组均为正显著。

同时,Geradino etc (2017) 还认为,即使充分增加了交叉项,即同时增加了 G,G*X,以及 G*X*Z 这些交叉项,仍然无法避免协变量组间可能存在系统性差异的问题。

1.2 rddsga 方法的主要原理:倾向得分+逆概率加权+RDD 回归+组间差异显著性检验

针对上述两个问题,Geradino etc (2017) 提出可以事先平衡样本组协变量的组间差异,使得个体处于哪一个子样本组近乎是是随机决定的。其平衡的思路为:

  • 1) 计算每个个体被分到该子样本组的概率,也即倾向得分值;
  • 2) 倾向得分值更高的个体,也即更大概率被分配到该子样本组的个体,在 RDD 回归中被赋予更小的权重,使得加权后子样本组之间的协变量统计特征不存在显著差异。(这就是基于倾向得分的逆概率加权方法的基本思路,以下我们简称 IPSW,Inverse Propensity Score Weighting)。

在完成了 IPSW 后,再通过充分加入指示变量与解释变量和协变量的交叉项,就通过分析交叉项系数的统计显著性来推断比较组间回归系数的差异问题。 Geradino etc (2017) 同时提供了命令来自动实现上述操作,这就是本文接下来将要介绍的 rddsga 命令。

注:关于 IPSW 方法原理的内容可参考赵西亮(2017)。

2. Stata 示例:直接分组 RDD 回归 v.s rddsga 回归

2.1 命令安装和语法介绍

安装最新版本命令:

. ssc install rddsga, replace

完成后,输入 help rddsga,即可查看其帮助文件。 rddsga 的语法格式为:

. rddsga depvar assignvar [indepvars][if] [in] [, options]

其中,

  • depvar 为被解释变量
  • assignvar 为驱动变量
  • indepvars 为协变量

opitions 主要包括四类:

  • 第一类为 subgroup 设置的选项,其中,
    • sgroup() 的括号里填分组指示变量
    • cutoff() 里填驱动变量里的断点,默认为 0
    • bwidt() 用于设定带宽,不同于普通 RDD 回归,这里需要事先指定带宽值。
  • 第二类为平衡的相关选项,其中
    • balance(varlist) 里面填计算倾向得分时方程右边的变量,默认等于 rddsga 回归时的全部协变量
    • probit 指定通过 Probit 估计计算倾向得分,默认为 logit
    • nocomsup 是指定对所有样本进行回归,不写则默认只对处在共同支撑域的样本进行回归。
  • 第三类为 RDD 回归选项,在这里不一一介绍;
  • 第四类为 rddsga 回归结果生成变量的有关选项。其中:
    • dibalance 会输出 IPSW 前和 IPSW 后样本组间差异的比较结果
    • comsup 会生成虚拟变量,赋值为 0 表示不在共同支撑域,为 1 则表示在共同支撑域
    • pscore 生成倾向得分值
    • ipsweight 是生成了基于 pscore 而生成的逆概率加权权重(ipsw)。

注:权重的生成公式见本文 3.1 部分的内容。参考赵西亮(2017)书中的推导,可知经过 IPSW 对抽样个体进行加权,可以获得总体特征的无偏估计。

2.2 模拟数据生成

为推文内容更清晰,生成模拟数据所需的代码放在附录1。读者也可以通过以下链接来直接下载: 数据下载

注:该链接里还包括了一份rddsga作者附带的数据rddsga_synth.dta,本文还将在第 3.5 节利用这份作者附带数据来简单演示rddsga命令的用法。

2.3 检验思路 1:根据分组指示变量 G,直接分样本 RDD 回归

我们先展示了直接分组进行 RDD 回归,不通过 IPSW 来事先平衡子样本组组间协变量差异的回归结果,作为和 rddsga 回归的对比。具体的,我们直接根据分组指示变量 G 来分组 RDD 回归。代码如下:

global dir "d:/RDDStata"   // 模拟数据存放于此,可以自行修改
cd $dir
clear all
use "rddsga_simu_data.dta",clear
global varlist z1 z2 //定义协变量

// 先全样本回归,获取全样本 RDD 回归对应的最优宽带
   /*
ssc install rdrobust, replace // 若未安装,请执行此行
   */
rdrobust y1 xc covs($varlist) 
est store Full
global h = e(h_l)

// 不平衡,直接分组回归
rdrobust y1 xc covs($varlist) , h($h)
est store  full
rdrobust y1 xc covs($varlist)  if G==0 , h($h)
est store  G0
rdrobust y1 xc covs($varlist)  if G==1 , h($h)
est store  G1

// 输出全样本和分样本的回归结果,并报告置信区间
* using regression.rtf
esttab full G0 G1, compress ///
       mtitle("Full Sample" "G=0" "G=1") ///
       nogap noconstant order() ci r2  replace ///
       starlevels(* 0.10 ** 0.05 *** 0.01) 

回归结果如下:

. esttab full G0 G1, compress ///
       mtitle("Full Sample" "G=0" "G=1") ///
       nogap noconstant order() ci r2  replace ///
       starlevels(* 0.10 ** 0.05 *** 0.01) 

----------------------------------------------------------
                      (1)             (2)             (3) 
              Full Sample             G=0             G=1 
----------------------------------------------------------
RD_Estim~e       0.915***        0.987***        0.844***
            [0.672,1.158]   [0.636,1.339]   [0.500,1.188] 
----------------------------------------------------------
N                    4000            2000            2000 
R-sq                                                      
----------------------------------------------------------
95% confidence intervals in brackets
* p<0.10, ** p<0.05, *** p<0.01

可以发现,在这个例子中,直接分组 RDD 回归结果存在问题一,即不同子样本回归虽然回归系数不同且都显著,但由于置信区间存在重叠(G0 组置信区间为 [0.636,1.339],而 G1 组置信区间为 [0.500,1.188]),所以不能认为 G0 组的影响更显著。

同时,直接分样本 RDD 回归可能还存在问题二:即不同组协变量存在显著的差异。为检验是否存在问题二,我们接下来逐个检验协变量 z1 和 z2 的组间特征差异。为此,我们借用 psmatch2 命令提供的辅助命令 pstest 进行检验。代码和检验如下所示,其中 G0 组为表格里的 Control 组,G1 组为表格里的 Treated 组。变量组间差异非常显著。尤其是对于 z1,组间系数偏差为 92.3%,t 检验 p 值为 0。z2 组间系数差异也比较大,偏差为 21.4%,p 值为 0。可见如果直接分组回归,会存在问题二:协变量组间存在显著差异。此外,也可以用ttable2命令,来同样实现检验子样本组间协变量差异的目的。

// 利用 pstest 命令检验不同子样本组协变量是否存在显著差异
pstest z*,raw t(G)

// 也可以用 ttable2 命令来检验
ttable2 z1 z2, by(G)

结果如下:

. pstest z*,raw t(G)

----------------------------------------------------------
         |       Mean             |    t-test    |  V(T)/
Variable | Treated Control  %bias |   t    p>|t| |  V(C)
---------+------------------------+--------------+--------
z1       | .48018   .01847   92.3 | 29.18  0.000 |  1.11*
z2       | 1.6236   .98905   21.4 |  6.77  0.000 |  1.00
----------------------------------------------------------
* if variance ratio outside [0.92; 1.09]

------------------------------------------------------------
Ps R2  LR chi2  p>chi2  MeanBias  MedBias    B     R   %Var 
------------------------------------------------------------
0.145   805.08   0.000    56.8     56.8    94.6*  1.09   50
------------------------------------------------------------
* if B>25%, R outside [0.5; 2]


. ttable2 z1 z2, by(G)
-----------------------------------------------------
Variables    G1(0)  Mean1     G2(1)  Mean2   MeanDiff
-----------------------------------------------------
z1           2000   0.018     2000   0.480  -0.462***
z2           2000   0.989     2000   1.624  -0.635***
-----------------------------------------------------

2.4 检验思路 2:使用 rddsga 命令回归

rddsga 回归检验的代码如下所示:

cd $dir
clear all
use "rddsga_simu_data.dta",clear

global varlist z1 z2  //定义协变量
qui rdrobust y1 xc covs($varlist)  
est store  full 
global h = e(h_l)   //获取全样本回归的带宽

// rddsga回归
rddsga y1 xc $varlist , balance($varlist) ///
       sgroup(G) bwidth($h) dibalance     ///
       comsup(_coms) ipsweight(_ipsw)     ///
       pscore(_ps) reducedform nobootstrap  

运行结果如下表所示,可知组间回归系数差值为 -.0775677,z 值为 -0.87,对应 p 值为 0.382,大于 0.1。可见在本例子中,并不能认为 G0 和 G1 组的因果效应存在显著差异。因而 rddsga 方法可以比较组间系数差异是否统计显著,从而解决了问题一。

|      因变量:y1        |   Coef.   | Std. Err. |   z   |  P>z  | [95% Conf. Interval]
---------+------------------------+--------------+---------------------------+--------
|         G0 组         |  1.05507  | .0853787  | 12.36 | 0.000 |  .8875399 1.222599   
|         G1 组         | .9775018  | .0243976  | 40.07 | 0.000 |  .9296289 1.025375   
| 组间差异:G1 组-G0 组  | -.0775677 | .0887962  | -0.87 | 0.382 |                     

此外,如果在回归命令里加入选项 dibalance,rddsga 会同时显示逆概率加权(IPSW)前后的比较,结果如下表所示。观察可知,对于 z1 而言,IPSW 前后,组间均值差异从 0.493-0.0506=.4424 下降到 0.442-0.453=-0.011,组间标准误也从-.853 下降到.0204,p 值增加到 0.831,可见经过平衡,z1 变量的差异已经大为缓解不显著。z2 也体现了类似的变化。可见,rddsga 回归方法可以缓解问题二,即降低了组间协变量可能存在的系统性差异。

| IPSW 之前 |G0 组均值  | G1 组均值      | std_diff     | p-value     |
---------+------------------------+------------+--------------+-------
|    z1     |   .0506   |     .493      |    -.853     |  1.70e-19   |
|    z2     |   1.12    |      1.6      |    -.157     |    .103     |

| IPSW 之后 | G0 组均值  | G1 组均值      | std_diff     | p-value     |
---------+------------------------+------------+--------------+-------
|    z1     |   .453    |     .442      |    -.0204    |     .831    |
|    z2     |   1.16    |     1.53      |    -.121     |     .205    |

注: 这里的「IPSW 前」,是指已经删除不处于共同支撑域样本,而尚未进行 IPSW 的时候。观察回归结果输出界面可知,处在共同支撑域内的样本只有(122+955)个,相对于删除前,少了(400-122)+(3600-955)个样本。

3. rddsga 回归的补充说明

3.1 回归命令生成变量的解释说明

上面使用 rddsga 语法命令里,我们通过 comsup(_coms) pscore(_ps) ipsweight(_ipsw) 三个选项分别生成了三个变量_coms,_ps 和 _ipsw ,分别表示是否在共同支撑域,倾向得分值和逆概率加权权重。这前两个指标指标和 PSM 方法中psmatch2命令 生成结果含义基本一致。对于第三个变量 ipsw,其计算公式为:

其中,p 为倾向得分值,本文里为回归生成的变量 ps。Gi=1,如果是 G1 组,Gi=0,如果是 G0 组。P(W)为平衡前处于 G1 的比例,本文例子中,值等于 955/(955+122)。公式的具体含义可参考 Gerardino et al. (2017) 或赵西亮 (2017)。

3.2 带宽选择的一致性问题

rddsga 回归中,要事先指定带宽,并且这个带宽要全样本 RDD 回归的带宽相同。在本文里,我们通过 global h = e(h_l) 方式来实现。

3.3 分组指示变量 G 的指定问题

读者可以看到,PSM-DID 和 rddsga 一样都是基于倾向指数来实现因果效应的识别。但对于前者,对于谁是对照组谁是处理组的理解,往往不存在歧义;而后者却并不一定。例如分所有制样本进行回归时,我们既可指定国有企业为对照组,也可指定非国有企业为对照组。在 rddsga 命令里,作者是把 G0 作为对照组,G1 作为控制组。对于上述两种指定,在计算逆概率加权权重(ipsw)时,样本量会出现一些差异,最后回归结果也并不对应。因此在使用时,要注意用对照组与处理组差异的方式解释,强调相对于 G0 组,G1 组的"处理效应"。读者可以在上例中,执行命令 replace G=1-G,来观察这种不同指定方式的区别。

3.4 多子样本组的 rddsga 回归问题

rddsga 命令可很容易扩展到多个子样本组系数比较的场景。例如对于分地区回归时,有东中西三个地区子样本组。就可以均以东部为对照组,然后分别以中部和西部对处理组。具体操作中,可设置两个分组指示变量 G 和 H,其中 G=1 表示中部地区,G=0 表示东部地区;H=1 表示西部地区,H=0 表示东部地区。然后分两次 rddsga 回归就可以。

3.5 使用rddsga命令自带数据来进行演示

本部分使用rddsga命令附带的数据rddsga_synth.dta来进行简单的回归演示。相关代码如下:

clear all
use rddsga_synth.dta  //数据下载地址见文末链接
des  

可知,该数据为一共有 6 个变量,各变量含义如下所示

----------------------------------------------------------------
              storage   display    value
variable name   type    format     label      variable label
----------------------------------------------------------------
Y               float   %9.0g                 Outcome
X               float   %9.0g                 Running variable
D               float   %9.0g                 Dreatment
G               float   %9.0g                 Subgroup
W1              float   %9.0g                 Covariate 1
W2              float   %9.0g                 Covariate 2
----------------------------------------------------------------

下面回归中, W1 为估算倾向得分时加入的协变量, G 为分组指示变量,带宽为 10 。diabl且需要输出组间协变量差异(即 W1 )在IPSW 前后的统计特征比较。reduced表示用OLS估计。

 rddsga Y X, balance(W1) sgroup(G) bwidth(10) dibal reduced

回归结果如下所示,可知在 IPSW 前, W1 在 G0 组与 G1 组之间存在显著差异( p 值几乎为 0) ,而之后则没有显著差异(p值几乎为 0.611) 。同时,G0 组回归系数为 -0.0231459 ,而 G1 组回归系数为 0.3356292。两组组间系数差为 0.3587751 ,但对应 p 值为 0.166 ,并不显著。

. rddsga Y X, balance(W1) sgroup(G) bwidth(10) reduced dibal

* Unweighted
----------------------------------------------------------
             |   mean_G0    mean_G1   std_diff    p-value 
-------------+--------------------------------------------
          W1 |      .722     -.0222       .777   2.53e-97 
----------------------------------------------------------
Obs. in subgroup 0: 1355
Obs. in subgroup 1: 1325
Mean abs(std_diff): .7770883
F-statistic: 476.31398
Global p-value: 0

Inverse Propensity Score Weighting
----------------------------------------------------------
             |   mean_G0    mean_G1   std_diff    p-value 
-------------+--------------------------------------------
          W1 |      .425       .406      .0194       .611 
----------------------------------------------------------
Obs. in subgroup 0: 1353
Obs. in subgroup 1: 1325
Mean abs(std_diff): .01936049
F-statistic: .25918458
Global p-value: .61072298

Bootstrap replications (50)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50


------------------------------------------------------------------------
           Y |     Coef.   Std. Err.    z   P>|z|   [95% Conf. Interval] (P)
-------------+----------------------------------------------------------
Subgroup     |                                                          
           0 | -.0231459   .1608766  -0.14  0.886  -.3386011    .2923094
           1 |  .3356292   .2235642   1.50  0.133  -.1027474    .7740058
-------------+----------------------------------------------------------
Difference   |  .3587751   .2588159   1.39  0.166  -.2380554    .9556056
------------------------------------------------------------------------

4. rddsga 回归的进阶:命令拆解

根据上文的介绍,我们可知 rddsga 实质上就是在 RDD 回归前,通过倾向得分+逆概率加权的方式来缓解子样本组间的协变量差异。因而其分析过程可分解为如下四步:1)计算倾向得分值;2)计算逆概率加权值(ipsw);3)断点回归;4)计算组间系数差异及统计显著性。

接下来,我们不使用 rddsga 命令,而是上述这四个步骤来复现上文中 rddsga 回归的结果。此外,我们同时也复现了rddsga方法下协变量组间差异的前后对比情况(即 2.4 节 第二和第三张表)。

4.1 计算倾向得分(ps),统计非共同支撑域(coms)样本

cls
clear all
cd $dir
use "rddsga_simu_data.dta",clear

global varlist z1 z2  // 定义协变量
qui rdrobust y1 xc covs($varlist)  
global h = e(h_l)   // 获取全样本回归的带宽

// 计算倾向得分
logit G z1 z2   if  (xc>=-$h)&(xc<=$h)
predict double _ps if (xc>=-$h)&(xc<=$h)
ereturn clear

4.2 计算逆概率加权权重值(ipsw

// 计算子样本组间的共同支撑域
qui sum _ps if G == 1 
gen _coms = (_ps>= `r(min)' & _ps <= `r(max)')

// 计算处于共同支撑域内的G0组和G1组的样本数量
  qui count if  (xc>=-$h)&(xc<=$h)&(G==0)&(_coms==1)
  local N_G0 = `r(N)'
  qui count if  (xc>=-$h)&(xc<=$h)&(G==1)&(_coms==1)
  local N_G1 = `r(N)'
  
// 计算逆加权概率
  qui gen _ipsw=`N_G1' /(`N_G1' +`N_G0' ) /_ps*(G==1)+`N_G0' /(`N_G1' +`N_G0' ) /(1-_ps)*(G==0)

4.3 RDD 回归并获取子样本组间系数

// 这里我们直接使用 OLS 回归命令,来实现和 RDD 回归同样的结果
  reg y1 i.G#1.T i.G i.G#(c.z1 c.z2 c.xc c.xc#i.T ) ///
    [iw=_ipsw] if  (xc>=-$h)&(xc<=$h)&(_coms==1),vce(robust)

上述命令中,比较关键的是通过 [iw=_ipsw]选项设定 OLS 回归的权重,通过 (xc>=-$h)&(xc<= $h)将样本限制在带宽以内,以及通过(_coms==1)将样本限制在共同支撑域内。回归结果中中 G#T 分类下的结果是我们关心的变量,也就是即 G0 组和 G1 组 RDD 回归的因果效应,如下表所示。其中(0 1)所在行表示 G=0,T=1,即 G0 组的处理效应;类似地,(1 1)所在行表示 G1 组的处理效应。对比第一个表的回归结果,可知完全相同。

|  y1    |     Coef.    Robust SE     t    p value  [95% Conf.Interval] 
---------+--------------------------------------------------------------
|   G#T  |                                                              
---------+--------------------------------------------------------------
|   0 1  |  1.047867    0.0771577   13.58     0    [0.8964683 1.199265] 
|   0 1  | 0.9715827    0.0771577   13.58     0    [0.8964683 1.199265] 

4.4 组间系数差异的统计显著性检验

接下来,我们构建检验组间系数差异的统计显著性。

qui matrix m = r(table)
qui scalar b1=m[1,2]   // G1组回归系数
qui scalar b0=m[1,1]   // G0组回归系数
qui scalar se1=m[2,2]  // G1组标准误
qui scalar se0=m[2,1]  // G0组标准误
  
scalar b_diff=  (b1-b0)
scalar se_diff =((se1^2+se0^2))^0.5
scalar t_diff = abs(b_diff/se_diff)
scalar p_diff = 2*(1-normal(abs(t_diff))) 

dis "子样本组间系数的差值、t值和p值分别为:"
scalar list b_diff t_diff p_diff

运行结果为:b_diff = -.0775682,z_diff = -0.87355314,p_diff = 0.38236166。对比第 2.4 节第一个表的结果,可知完全相同。

到此,我们就完成了复现 rddsga 回归结果的任务。

4.5 IPSW 前后协变量的组间差异检验

此外,除了输出子样本组之间回归系数差异的统计显著性外,rddsga 另外一个特点提供了基于倾向得分的逆概率加权方法的,平衡子样本组之间协变量差异的效果,即本文 2.4 的第二个和第三个表。下面我们将复现这部分结果。

由于 Stata 的ttest命令无法比较加权样本的变量差异,在这里我们无法直接加入[iw=_ipsw]来进行 IPSW 后的组间 t 检验,所以我们采用reg命令来间接地实现。本部分的代码如下所示:

//协变量组间差异的比较:IPSW前
qui reg z1  if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==0)
qui matrix m = r(table)
qui scalar b0_nw=m[1,1] // G0组样本均值

qui reg z1  if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==1)
qui matrix m = r(table)
qui scalar b1_nw=m[1,1]  // G1组样本均值

qui reg z1 i.G  if (xc>=-$h)&(xc<=$h)&(_coms==1)
qui matrix m = r(table)
qui scalar bdiff_nw=m[1,2] // G1组样本均值-G0组样本均值
qui scalar pval_nw = m[4,2] // p-value 


//协变量组间差异的比较: IPSW后
qui reg z1 [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==0)
qui matrix m = r(table)
qui scalar b0_ipsw=m[1,1] // G0组样本均值

qui reg z1 [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==1)
qui matrix m = r(table)
qui scalar b1_ipsw=m[1,1]  // G1组样本均值

qui reg z1 i.G [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)
qui matrix m = r(table)
qui scalar bdiff_ipsw=m[1,2] // G1组样本均值-G0组样本均值
qui scalar pval_ipsw = m[4,2] // p-value 


//显示平衡前后的统计特征
qui dis "z1变量IPSW前:G0组样本均值、G1组样本均值、均值组间差异、p值"
scalar list b0_nw b1_nw bdiff_nw  pval_nw

qui dis "z1变量IPSW后:G0组样本均值、G1组样本均值、均值组间差异、p值"
scalar list b0_ipsw b1_ipsw bdiff_ipsw  pval_ipsw

5. 参考资料

6. 附:全文 dofile

6.1 附录1:生成模拟数据的代码

* ### 2.2 生成模拟数据
*-Note:
*  下面的例子采用蒙特卡洛模拟生成数据,
*  你只需执行命令便可以自动产生数据文件

global dir "d:/RDDStata"   // 模拟数据存放于此,可以自行修改
capture mkdir $dir
cd $dir

clear 
set obs 4000
set seed 123

gen x = runiform()  //分配变量
gen xc = x-0.5      //分配变量去中心化

gen e = rnormal()/5     // noise
gen z1 = rnormal()*0.5  //控制变量z1
gen z2=1+3*invnormal(uniform())+sin(x*5)/3+e  // 另一个控制变量z2

gen T=0               
replace T=1 if x>0.5   //treatment 

gen g0 = 0 + 3*log(x+1) + sin(x*6)/3
gen g1 = T + 3*log(x+1) + sin(x*6)/3
gen y1 = g1 + 0.5*z1 +0.3*z2+ e   // outcome variable,with    cutoff effect
gen y0 = g0 + 0.5*z1 +0.3*z2+ e   // outcome variable, without cutoff effect

label var y1 "Outcome variable (y)"
label var y0 "Outcome variable (y)"
label var x  "Assignment variable (x)"
label var xc "Centered Assignment variable (x-c)"
label var T  "T=1 for x>0.5, T=0 otherwise"

drop e g* 
save "RDD_simu_data0.dta", replace  

*  对数据RDD_simu_data0进行调整,以更好展示rddsga命令的回归结果
clear all
use RDD_simu_data0, clear

set seed 1234
gen t=rnormal()
sum t,detail
local t_50=r(p50)
gen G=0   
replace G=1 if t>=`t_50'

// 调整z1和z2的赋值
replace z1=z1+0.5*G    
replace z2=z2+0.5*G

drop y0 x t
label var G "subgroup indicator variable"
label var z1 "control variable 1"
label var z2 "control variable 2"

save "rddsga_simu_data.dta", replace 

6.2 附录2:本推文所使用到的其它代码


* ### 2.3 检验思路 1:根据分组指示变量 G,直接分样本 RDD 回归

cd $dir
clear all
use "rddsga_simu_data.dta",clear
global varlist z1 z2 //定义协变量

// 先全样本回归,获取全样本 RDD 回归对应的最优宽带
   /*
ssc install rdrobust, replace // 若未安装,请执行此行
   */
rdrobust y1 xc covs($varlist) 
est store Full
global h = e(h_l)

// 不平衡,直接分组回归
rdrobust y1 xc covs($varlist)  if G==0 , h($h)
est store  G0
rdrobust y1 xc covs($varlist)  if G==1 , h($h)
est store  G1
rdrobust y1 xc covs($varlist)  , h($h)
est store  full

// 输出全样本和分样本的回归结果,报告置信区间而非标准误
esttab full G0 G1, compress ///
       mtitle("Full Sample" "G=0" "G=1") ///
       nogap noconstant order() ci r2  replace ///
       starlevels(* 0.10 ** 0.05 *** 0.01) 
	  
// 利用 pstest 命令检验不同子样本组协变量是否存在显著差异
/*
ssc install pstest, replace //外部命令
*/

pstest z*, raw t(G)

// 也可以用 ttable2 命令来检验
ttable2 z1 z2, by(G)



* ### 2.4 检验思路 2:使用 rddsga 命令回归
cd $dir
clear all
use "rddsga_simu_data.dta",clear

/*
  ssc install rddsga, replace //外部命令,执行此行安装
*/

global varlist z1 z2  //定义协变量
qui rdrobust y1 xc covs($varlist)  
est store  full 
global h = e(h_l)   //获取全样本回归的带宽

// rddsga回归
rddsga y1 xc $varlist , balance($varlist) ///
       sgroup(G) bwidth($h) dibalance     ///
       comsup(_coms) ipsweight(_ipsw)     ///
       pscore(_ps) reducedform nobootstrap 


* ### 4. rddsga 回归的进阶理解:命令拆解
* #### 4.1 计算倾向得分(ps),统计非共同支撑域(coms)样本
cls
clear all
cd $dir
use "rddsga_simu_data.dta",clear

global varlist z1 z2  // 定义协变量
qui rdrobust y1 xc covs($varlist)  
global h = e(h_l)   // 获取全样本回归的带宽

// 计算倾向得分
logit G z1 z2   if  (xc>=-$h)&(xc<=$h)
predict double _ps if (xc>=-$h)&(xc<=$h)
ereturn clear


* #### 4.2 计算逆概率加权权重值(ipsw)
// 计算子样本组间的共同支撑域
qui sum _ps if G == 1 
gen _coms = (_ps>= `r(min)' & _ps <= `r(max)')

// 计算处于共同支撑域内的G0组和G1组的样本数量
  qui count if  (xc>=-$h)&(xc<=$h)&(G==0)&(_coms==1)
  local N_G0 = `r(N)'
  qui count if  (xc>=-$h)&(xc<=$h)&(G==1)&(_coms==1)
  local N_G1 = `r(N)'
  
// 计算逆加权概率
  qui gen _ipsw=`N_G1' /(`N_G1' +`N_G0' ) /_ps*(G==1)+`N_G0' /(`N_G1' +`N_G0' ) /(1-_ps)*(G==0)
  
  
  
* #### 4.3 RDD 回归并获取子样本组间系数
// 这里我们直接使用 ols 回归命令,来实现和 RDD 回归同样的结果
  reg y1 i.G#1.T i.G i.G#(c.z1 c.z2 c.xc c.xc#i.T ) ///
    [iw=_ipsw] if  (xc>=-$h)&(xc<=$h)&(_coms==1),vce(robust)
    
* #### 4.4 组间系数差异的统计显著性检验
qui matrix m = r(table)
qui scalar b1=m[1,2]   // G1组回归系数
qui scalar b0=m[1,1]   // G0组回归系数
qui scalar se1=m[2,2]  // G1组标准误
qui scalar se0=m[2,1]  // G0组标准误
  
scalar b_diff=  (b1-b0)
scalar se_diff =((se1^2+se0^2))^0.5
scalar t_diff = abs(b_diff/se_diff)
scalar p_diff = 2*(1-normal(abs(t_diff))) 

dis "子样本组间系数的差值、t值和p值分别为:"
scalar list b_diff t_diff p_diff


* #### 4.5 IPSW 前后协变量的组间差异检验    
//协变量组间差异的比较:IPSW前
qui reg z1  if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==0)
qui matrix m = r(table)
qui scalar b0_nw=m[1,1] // G0组样本均值

qui reg z1  if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==1)
qui matrix m = r(table)
qui scalar b1_nw=m[1,1]  // G1组样本均值

qui reg z1 i.G  if (xc>=-$h)&(xc<=$h)&(_coms==1)
qui matrix m = r(table)
qui scalar bdiff_nw=m[1,2] // G1组样本均值-G0组样本均值
qui scalar pval_nw = m[4,2] // p-value 


//协变量组间差异的比较: IPSW后
qui reg z1 [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==0)
qui matrix m = r(table)
qui scalar b0_ipsw=m[1,1] // G0组样本均值

qui reg z1 [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)&(G==1)
qui matrix m = r(table)
qui scalar b1_ipsw=m[1,1]  // G1组样本均值

qui reg z1 i.G [iw=_ipsw] if (xc>=-$h)&(xc<=$h)&(_coms==1)
qui matrix m = r(table)
qui scalar bdiff_ipsw=m[1,2] // G1组样本均值-G0组样本均值
qui scalar pval_ipsw = m[4,2] // p-value 


//显示平衡前后的统计特征
qui dis "z1变量IPSW前:G0组样本均值、G1组样本均值、均值组间差异、p值"
scalar list b0_nw b1_nw bdiff_nw  pval_nw

qui dis "z1变量IPSW后:G0组样本均值、G1组样本均值、均值组间差异、p值"
scalar list b0_ipsw b1_ipsw bdiff_ipsw  pval_ipsw

7. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh RDD PSM
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

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

连享会小程序:扫一扫,看推文,看视频……

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

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

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