合成控制法R及Stata操作和应用

计量经济学服务中心专辑汇总！计量百科·资源·干货：

Stata  |Python  |Matlab  |Eviews  |R 

Geoda  |ArcGis  |GeodaSpace  |SPSS 

一文读懂   |数据资源   |回归方法  |网络爬虫                                                               

门限回归   |工具变量   |内生性   |空间计量 

因果推断   |合成控制法   |倾向匹配得分   |断点回归   |双重差分   

面板数据  | 动态面板数据

香烟控制法在美国有着悠久的历史，至少可以追溯到1893年，当时的华盛顿成为第一个禁止香烟销售的州。在在接下来的30年里，其他15个州也相继出台了类似的禁烟措施(Dinan和Heckelman 2005)。这些早期的反烟草法主要是出于道德考虑， 健康问题是次要的(Tate 1999)。近100年来，后来，在这些早期的法律被废除很久之后，对吸烟的健康风险的认识有了一个新的广泛传播，美国各州和联邦政府纷纷出台反烟草法以及最终的海外市场。1988年，加州选民发起的99号提案，第一个现代化的大规模控烟项目，这也就是著名的99香烟法案。

本文基于R和Stata为大家讲解如何进行合成控制法相关操作及应用。

1

背景介绍

99法案将加州的卷烟消费税每包提高25美分，这部分指定税收将用于健康和反吸烟教育预算，资助反吸烟媒体，并推动了当地的室内空气清洁条例遍及全州((Siegel 2002)。在最初的实现中, 99号提案每年在学校、社区、县、州一级开展反烟草项目产生超过1亿美元的收入 。每年几乎有2000万美元可用于烟草有关的研究。正如Glantz和Balbach(2000)所说，他说:“这些项目使任何其他州或地区相形见绌 ，联邦政府也曾经对烟草采取过措施。” 99号提案引发了一波当地清洁空气条例的浪潮，在加州。在99号提案之前，加州没有城镇要求餐厅100%禁烟。从1989年到 2000年通过了大约140部这样的法律(Siegel 2002)。到1993年，当地颁布法令禁止在工作场所吸烟，保护了加州近三分之二的工人(Glantz) 和Balbach 2000)。1994年，加利福尼亚州通过了附加法案 禁止在封闭工作场所吸烟的立法。到1996年，超过90%的加州工人获得了医疗保险，通过无烟工作场所政策(西格尔2002)。不吸烟的权利倡导者们看到了一波又一波地方法令通过。

1988年11月美国加州通过了当代美国最大规模的控烟法（anti-tobacco legislation），并于1989年1月开始生效。该法将加州的香烟消费税（cigarette excise tax）提高了每包25美分，将所得收入专项用于控烟的教育与媒体宣传，并引发了一系列关于室内清洁空气的地方立法（local clean indoor-air ordinances），比如在餐馆、封闭工作场所等禁烟。Abadie et al. (2010)根据美国1970-2000年的州际面板数据，采用合成控制法研究美国加州1988年第99号控烟法（Proposition 99）的效果。

合成控制法Stata应用

2

语法格式

合成控制法命令为：synth，安装方法为

ssc install synth

基本语法格式为（下述四行为一行完整命令，出现断行了）：

synth depvar predictorvars , trunit(#) trperiod(#)  ///[counit(numlist) xperiod(numlist) mspeperiod()   /// resultsperiod() nested allopt unitnames(varname) /// figure keep(file) customV(numlist) optsettings ]

3

选项含义介绍

depvar 表示结果变量（outcome variable）

predictorvars为预测变量（predictors）（ the list of predictor variables）（ By default, all predictor variables are averaged over the entire pre-intervention period, which ranges from the earliest time period available in the panel time variable specified in tsset timevar to the period immediately prior to the intervention specified in trperiod. Missing    values are ignored in the averages.）

trunit(#)用于指定处理地区（trunit表示 treated unit）。trunit(#) the unit number of the unit affected by the intervention as given in the panel id variable specified in tsset panelvar（表示需要前期进行面板数据设定）; see tsset. Notice that only a single unit number can be specified. If the intervention of interest affected several units the user may chose to combine these units first and then treat them as a  single unit affected by the intervention.（表示一次只能一个地区受到影响，若是有多个地区受到影响，需要前期把其合成）

trperiod(#)用于指定政策干预开始的时期（trperiod表示 treated period）（the time period when the intervention occurred）。

counit(numlist)用于指定潜在的控制地区（即donor pool，其中counit表示 control units），默认为数据集中的除处理地区以外的所有地区。

xperiod(numlist)用于指定将预测变量（predictors）进行平均的期间，默认为政策干预开始之前的所有时期（theentire pre-intervention period）。

mspeperiod()用于指定最小化均方预测误差（MSPE）的时期，默认为政策干预开始之前的所有时期。

figure表示将处理地区与合成控制的结果变量画时间趋势图，而选择项“resultsperiod()”用于指定此图的时间范围（默认为整个样本期间）。

nested表示使用嵌套的数值方法寻找最优的合成控制（推荐使用此选项），这比默认方法更费时间，但可能更精确。在使用选择项“nested”时，如果再加上选择项“allopt”（即“nested allopt”），则比单独使用“nested”还要费时间，但精确度可能更高。

keep(filename)将估计结果（比如，合成控制的权重、结果变量）存为另一Stata数据集（filename.dta），以便进行后续计算。

replace 表示将keep(filename) 中的data若是存在可以对其进行替换

keep(filename)中的数据主要包括如下内容：

 _time: 一个变量，包含resultsperiod()中指定的所有时间段的相应时间段(来自tsset面板时间变量(timevar))。

 _Y_treated: 对于resultsperiod()中指定的每个时间段，tr()中指定的处理单元的观察结果depvar。

_Y_synthetic: 估计合成控制单元的估计结果

_Co_Number: 包含在co()中指定的每个控制单元的单元号(来自tsset面板单元变量(panelvar))的变量。

_W_weight:包含co()中指定的每个控制单元的估计单位权重的变量。

更多选择项，详见help synth。

4

案例介绍：99法案香烟控制法案政策效应影响

1、加载示例数据:这个面板数据集包含了美国39个州在1970-2000年间的信息(参见Abadie、Diamond和Hainmueller(2010)的详细信息)。

sysuse synth_smoking

2、面板数据设定

xtset state year

查看数据结果截图：

相关数据含义如下：

该研究的结果变量为cigsale（人均香烟消费量，包/年）

预测变量包括retprice（平均香烟零售价格）、lnincome（人均收入对数）、age15to24（15-24岁人口所占总人口比重，年轻人为吸烟主力）、beer（人均啤酒消费量，烟酒不分家）。这些预测变量均为1980-1988年的州平均值，另外还使用1975、1980与1988年的人均香烟消费量作为三个额外的预测变量。另外，面板变量为 state（州），而时间变量为 year（年）。

3、合成控制法操作

synth cigsale retprice lnincome age15to24 beer cigsale(1975) cigsale(1980) cigsale(1988), trunit(3)trperiod(1989) xperiod(1980(1)1988) figure nested keep(smoking_synth)

上述代码全部为：synth cigsale retprice lnincome age15to24 beer cigsale(1975) cigsale(1980) cigsale(1988), trunit(3)trperiod(1989) xperiod(1980(1)1988) figure nested keep(smoking_synth)

相关含义表示为：

其中，“cigsale(1975)cigsale(1980) cigsale(1988)”分别表示人均香烟消费在1975、1980与1988年的取值。

必选项“trunit(3)”表示第3个州（即加州）为处理地区；

必选项“trperiod(1989)”表示控烟法在1989年开始实施。

选择项“xperiod(1980(1)1988)”表示将预测变量在1980-1988年期间进行平均，其中“1980(1)1988”表示始于1980年，以1年为间隔，而止于1988年；其效果等价于“1980 1981 1982 1983 1984 1985 1986 1987 1988”，而前者的写法显然更为简洁。

选择项“keep(smoking_synth)”将估计结果存为Stata数据集smoking_synth.dta（放在Stata的当前工作路径）。估计结果如下。

上表显示，大多数州的权重为0，只有以下五个州的权重为正，即Colorado (0.161)，Connecticut (0.068)，Montana (0.201)，Nevada (0.235)与Utah (0.335)，此结果与Abadie et al. (2010)汇报的结果非常接近（细微差别或由于计算误差）。

上表显示了合成中的每个控件状态的权重，表中报告的权重表明：加州的吸烟趋势在提案通过之前最好是由这Colorado 、Connecticut、Montana、Nevada 与Utah表示。其他状态的w权值为零。

考察加州与合成加州的预测变量是否接近：

在表1中，比较了实际的加利福尼亚和合成的加利福尼亚的预测变量均十分接近，故合成加州可以很好地复制加州的经济特征。

然后比较二者的人均香烟消费量在1989年前后的表现：

上图显示了加州和合成加州在1970-2000年期间的人均烟草销售额。注意，合成加州的人均销售量非常好的拟合追踪了加州整个99号提案前的变化趋势。

高程度的拟合表明，在1989年控烟法之前，合成加州的人均香烟消费与真实加州几乎如影相随，表明合成加州可以很好地作为加州如未控烟的反事实替身。在控烟法实施之后，加州与合成加州的人均香烟消费量即开始分岔，而且此效应越来越大。

我们估计第99号提案对香烟消费的影响分水岭是在99通过这个时间点，该法案一通过两条线开始明显地分开。

而香烟消费在合成加州继续其温和下降趋势来看，真正的加州经历了大幅下滑。这两条线之间的差异表明第99号提案对人均香烟销量的有很大的负面影响。

上图表明99号提案对人均香烟销量有很大影响，这种效应随着时间的推移而增强。我们的结果表明，对于整个1989-2000年期间的卷烟每包的消费量平均减少了近20包，人均下降约25%。

4、加州与合成加州差异对比

use 计量经济学服务中心smoking_synth.dta, cleargen effect= _Y_treated - _Y_syntheticlabel variable _time "year"label variable effect "gap in per-capita cigarette sales (in packs)"line effect _time, xline(1988,lp(dot) lc(black)) yline(0,lp(dash) lc(black)) ///text(-25 1988 "Passage of Proposition 99 {&rarr} ", placement(sw)) /// placement() et a. options for text()xtitle(year) ytitle("gap in per-capita cigarette sales (in packs)") lc(black) xlabel(1970(5)2000) ///title("政策效应：加州与合成加州人均香烟销售量差异") /// subtitle("2017.10.01")  ///note("计量经济学服务中心2017.10")

为了评估我们的结果的稳健性，我们使用了其他其他预测变量，包括各州的失业率，收入不平等，贫困，福利转移，犯罪率，毒品相关的逮捕率，香烟税，人口密度，以及种族等变量，发现结果都是稳健的。

5

安慰剂检验

为了评估我们的估计的显著性，我们提出：我们的结果是否是随机偶然的。如果我们随机选择一个州而不是加州来进行研究? 我们多久可以得到这样一个显著的结果，为了回答这个问题，我们使用安慰剂测试。

类似于 Abadie和Gardeazabal (2003)， Bertrand, Duflo和Mullainathan (2004)，我们使用合成药物进行安慰剂研究的方法，具体应用到香烟控制法案研究。

“安慰剂”（placebo）一词来自医学上的随机实验，比如要检验某种新药的疗效。此时，可将参加实验的人群随机分为两组，其中一组为实验组，服用真药；而另一组为控制组，服用安慰剂（比如，无用的糖丸），并且不让参与者知道自己服用的究竟是真药还是安慰剂，以避免由于主观心理作用而影响实验效果，称为“安慰剂效应”（placebo effect）。

为此，Abadie et al. (2010)进行了一系列的安慰剂检验，依次将donor pool中的每个州作为假想的处理地区（假设也在1988年通过控烟法），而将加州作为控制地区对待，然后使用合成控制法估计其“控烟效应”，也称为“安慰剂效应”。通过这一系列的安慰剂检验，即可得到安慰剂效应的分布，并将加州的处理效应与之对比。

在上图中，黑线表示加州的处理效应（即加州与合成加州的人均香烟消费之差），而灰线表示其他19个控制州的安慰剂效应（即这些州与其相应合成州的人均香烟消费之差）。显然，与其他州的安慰剂效应相比，加州的（负）处理效应显得特别大。假如加州的控烟法并无任何效应，则在这20个州中，碰巧看到加州的处理效应最大的概率仅为1/20 = 0.05，而这正好是常用的显著性水平。

政策实施前均方预测误差的平方根

tempname resmat forvalues i = 1/4 { synth cigsale retprice cigsale(1988) cigsale(1980) cigsale(1975) , trunit(`i') trperiod(1989) xperiod(1980(1)1988) matrix `resmat' = nullmat(`resmat') \ e(RMSPE) local names `"`names' `"`i'"'"' } mat colnames `resmat' = "RMSPE" mat rownames `resmat' = `names'        matlist `resmat' , row("Treated Unit")

合成控制法R应用

1

R语法格式

合成控制法命令为：synth，安装方法为

install.packages('synth')

导入该安装包为

library(Synth)

2

命令描述

在Abadie and Gardeazabal(2003)和Abadie, Diamond, Hainmueller(2010, 2011, 2014)的比较案例研究中实现因果推理的综合控制方法。synth通过比较受干预影响的单位的总体结果演变与合成对照组的相同总体结果演变来估计干预的效果。

synth通过搜索控制单元的加权组合来构建这个综合控制组，选择这些控制单元来根据预测结果的特征来近似受干预影响的单元。合成控制组结果的演变是对在没有干预的情况下对受影响单位观察到的反事实的估计。

synth还可用于进行各种安慰剂和置换试验，从而产生信息推断，而不考虑可用比较单位的数量和可用时间周期的数量。详见Abadie和Gardeazabal(2003)、Abadie、Diamond和Hainmueller(2011,2011,2014)。

synth要求用户提供四个矩阵作为它的主要参数。这些矩阵被命名为X0 X1 Z1 Z0。X1和X0分别包含处理单元和控制单元的预测值。Z1和Z0分别为治疗单元和对照组干预前时间的结果变量。预干预期是指干预前的一段时间，在此期间均方根预测误差(MSPE)应达到最小。MSPE是指在Z1和Z0中规定的所有干预前期间，处理单元和综合控制单元的结果之间的平方偏差。

从(面板)数据集创建矩阵X1、X0、Z1和Z0可能很繁琐。因此，Synth库提供了一个名为dataprep的预备函数，它允许用户轻松创建Synth所需的所有输入。通过首先调用dataprep，用户创建了一个名为data.prep.obj的对象，它包含运行synth所需的所有基本数据元素。

因此，实现综合控制方法的通常命令序列是首先调用dataprep来准备要加载到synth中的数据。然后调用synth构建合成控制组。最后，利用函数合成对结果进行了总结。包括如下函数：synth.tab, path.plot, or gaps.plot.

dataprep的文档中提供了这个序列的示例。强烈推荐这种方法。或者，用户可以提供自己的预处理数据矩阵，并通过X0、X1、Z1和Z0参数将它们加载到synth中。

synth的输出是一个列表对象，其中包含预测变量(solution.V)的权重和控制单元(solution.W)的权重，它们定义了对合成控制单元的贡献。

3

合成控制法语法格式：

代码为

dataprep(foo = NULL, predictors = NULL, predictors.op = "mean", special.predictors = NULL, dependent = NULL, unit.variable = NULL, time.variable = NULL, treatment.identifier = NULL, controls.identifier = NULL, time.predictors.prior = NULL, time.optimize.ssr = NULL, time.plot = time.optimize.ssr,   unit.names.variable = NA)

选项含义为：

4

案例

1、导入数据并查看数据

2、从为synth()提供输入的面板数据创建矩阵

3、运行synth命令以确定为被处理对象创建可能的最佳合成控制单元的权重。

代码为

synth.out <- synth(dataprep.out)

4、获取单元与合成单元预测值

代码为

synth.tables <- synth.tab( dataprep.res = dataprep.out, synth.res = synth.out)print(synth.tables)

结果为：

5、绘图

使用path.plot()和gap .plot()命令，可以获得经过处理的和综合控制单元的结果轨迹的总结图

代码为

path.plot(dataprep.res = dataprep.out,synth.res = synth.out)
gaps.plot(dataprep.res = dataprep.out,synth.res = synth.out)


结果为：

5

加州案例

1、代码为

library(Synth) 
# dataprep: prepare data for synthdataprep.out <- dataprep( foo = smoking, predictors= c("lnincome" ,"age15to24" ,"retprice" ) ,predictors.op = "mean" ,dependent = "cigsale" ,unit.variable = "id" ,time.variable = "year" ,special.predictors = list( list("cigsale", 1975, "mean"), list("cigsale", 1980, "mean"), list("cigsale", 1988, "mean") ), treatment.identifier = 3 ,controls.identifier = c(1,2,4:39) ,time.predictors.prior =c(1980:1988) ,time.optimize.ssr = c(1970:1987) ,unit.names.variable = "name" ,time.plot = 1970:2000 )

结果为：

2、进行合成分析

# run synthsynth.out <- synth(data.prep.obj = dataprep.out)


结果为：

3、获得合成分析表结果并进行打印输出

# Get result tablessynth.tables <- synth.tab( dataprep.res = dataprep.out, synth.res = synth.out)

结果为：

5、合成控制法结果表1

# plot results:# pathpath.plot(synth.res = synth.out, dataprep.res = dataprep.out, Ylab = c("gap in per-capita cigarette sales (in packs)"), Xlab = c("year"), Legend = c("加州","合成加州"), Legend.position=c("topright"), Main = NA, Z.plot = FALSE)

结果为：

6、合成控制法结果图2：加州与合成加州的差距图

## gapsgaps.plot(synth.res = synth.out, dataprep.res = dataprep.out, Ylab = c("gap in per-capita cigarette sales (in packs)"),          Xlab = c("year"), tr.intake = NA, Ylim = NA, Z.plot = FALSE))

结果为：

6

References

Abadie, A., Diamond, A., Hainmueller, J. (2014). Comparative Politics and the Synthetic Control Method. American Journal of Political Science Forthcoming 2014.

Synthetic : An R Package for Synthetic Control Methods in Comparative Case Studies. Journal of Statistical Software 42 (13) 1–17.

Abadie, A., Diamond, A., Hainmueller, J. (2011). Synth: An R Package for Synthetic Control Methods in Comparative Case Studies. Journal of Statistical Software 42 (13) 1–17.

Abadie A, Diamond A, Hainmueller J (2010). Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of California's Tobacco Control Program. Journal of the American Statistical Association 105 (490) 493–505.

Abadie, A. and Gardeazabal, J. (2003) Economic Costs of Conflict: A Case Study of the Basque Country American Economic Review 93 (1) 113–132.