【应用计量系列149】Sant' Anna的DID规定动作（4-2）：处理配置机制～一场星际会议的讨论

从今天开始，我们搬运Scott Cunningham的“Pedro’s diff in diff checklist”【有兴趣的人可以点击这里关注他的blog】。

大家也许还记得，Sant' Anna有一个十步DID：

第一步：画出处理模式图

第二步：呈现每个处理类别的处理个体数量

第三步：画出每个处理类别的平均结果时序图

第四步：【应用计量系列148】Sant' Anna的DID规定动作（4-1）：写好制度背景---处理配置机制

Cunningham看了一部科幻电视剧《The Expanse》，它讲述火星、地球和小行星带上的人们之间的太空冲突。

在第127届太阳系联盟SSF峰会上，各星球代表一起讨论太阳系经济和政治问题。一群社会科学家介绍了Equilift现金补助计划对收入的影响。TIT（泰坦理工学院）的一位青椒走到台前，挥了挥手，一个带有slides的全息影像图就漂浮在房间中央，供所有人观看。

“尊敬的代表们，大家好，”她开始说道。“我是莱拉·索利斯博士，一位未来学家，也是TIT奥利·阿什费尔特社会商业学教授。今天，我将介绍我们对 EquiLift 计划对五个殖民地影响的研究结果。我的演讲题目是‘评估EquiLift ：使用DID和基于结果的处理配置分析收入变化’。”

她讲话时，来自火星、木卫二、土卫六、土卫二和木卫三的代表们全神贯注地看着她。“众所周知，我们的殖民地遭受着严重的不平等。我们中没有一个人能免受这个问题的影响。社会共同体学告诉我们，当不平等现象过于严重时，内乱和两极分化就可能随之而来。我们有责任确保这种情况不会发生。目前，我们已经很接近了，因为我们中的一些人非常富有，而另一些人则非常贫穷。”

她停顿了一下，让大家感受到她话语的严肃性。“所以，三年前，你们所有人都批准了一项针对极度贫困人口的雄心勃勃的计划，名为EquiLift ，该计划将向上一年收入低于 25% 的人们提供食品券、住房券和职业培训。”

全息图显示了处理组和对照组的基线收入分布，揭示了明显的差异。“我们在本文中的目标是估计EquiLift对项目参与者参加两年后的收入影响。我们的主要估计将使用DID，包括事件研究。”

这张幻灯片刚放出来，来自火星的代表就打断了索利斯博士，他举手提问：“索利斯 博士，我最近读了一篇学术文章，是地球大战之前的两位经济学家写的，The Promise and Pitfalls of Differences-in-Differences: Reflections on 16 and Pregnant and Other Applications,’ an academic article by Ariella Kahn-Lang and Kevin Lang。

火星人继续说，“在这篇文章中，Kahn-Lang 和 Lang 写道：

“任何 DiD 论文都应该解决为什么处理组和对照组的原始水平不同，以及为什么这不会影响平行趋势”（Kahn-Lang 和 Lang，2019）。

火星人指着索利斯博士的slides说到，“显然，这两组人在基线时有很大的不同。如果我没记错的话，这不是一个随机实验，这无疑解释了你们两组之间的显著差异。处理组非常贫穷，在参加该计划的前一年，他们的年薪是 125 scripts，而对照组的收入几乎是 660 scripts。这似乎是一个针对一群明显与其他人截然不同的人的计划。”

另一位戴着眼镜的中年妇女举起了手。观众安静下来，以便她能发言。索利斯博士点点头，“请讲，罗伊Roy总理。”

“众所周知，计算社会学是计算机科学、社会学、经济学和政治学的后代。它们的交叉合并充分利用了每个学科的优势，包括经济学家对根据潜在结果进行处理选择选择时使用DID的怀疑。如果个人根据他们对项目的预期收益进行分类，那么我根本看不出这项经验研究如何识别这种因果关系。”

索利斯博士回应道，“各位代表，你们说得完全正确。正如我们在计量经济学课程中所学到的知识，在使用DID前，必须要了解处理配置机制。事实上，DID并不与所有形式的处理配置机制相吻合，这也许就是为什么 Kahn-Lang 和 Lang (2019) 指出事件研究本身并不能解决平行趋势是否成立的问题。我们通过几种方式进行了讨论。

首先，我们注意到，罗伊 (Roy) 总理所提到的该计划是否会对参与者有利的分类在这里并不适用。将个人分配到该计划中并不是基于人的欲望（预期到未来是否有利）。诚然，参与者的需求才是关键，但这种需求是用一个分数/标准来决定——他们的基线收入，我们称之为 Y(0)，基“未处理的潜在结果”。参见【应用计量系列148】Sant' Anna的DID规定动作（4-1）：写好制度背景---处理配置机制。我们注意到的第一件事是，该计划使用的处理配置机制确实是基于未处理的潜在结果y0。

该计划采用的选拔机制是评分。即基线的个人按其基线收入进行排名，任何收入低于其所在殖民地收入分布第 25 分位的人都会被纳入 EquiLift 计划。你们可以在下图看到 SSF 中的不平等程度以及资格阈值。

“但我们做的第二件事是检验基线工人的构成，看看我们的处理组和对照组是否存在不同的群体构成。我们使用 OLS 回归的证伪方法：

其中，i表示工人，c表示地区，表示个体固定效应。上述ols回归没有协变量，如果带有协变量，就是条件平行趋势。结果如下表所示：

事实上，虽然殖民地内部存在相当大的异质性，但我们的样本本身在种族、年龄和性别方面非常相似，我们面临的不平等问题也是如此。

如你们所见，只有一个other回归结果在 10% 或更低的水平上具有统计显著性。其余大多数都是不显著异于0。即使是年龄，其系数为 -0.86，也是零，因为在我们的数据中样本平均值为 18.40。虽然这不能证明平行趋势，但这表明处理与除基线贫困和富有（收入）之外的任何特定特征都没有系统相关性。

最后，用事件研究。虽然事件研究本身并不能证明平行趋势假设，但几个世纪以来，它一直是社会科学家评估平行趋势假设的主要方式——通过检验事件研究中的处理前系数。尽管 Kahn-Lang 和 Lang (2019) 认为平行趋势既不是“证明”平行趋势成立的必要条件，也不是充分条件，但它是我们决定使用 25% 以上作为对照组的基础。换句话说，这张图片让我们有信心使用较富裕的群体作为我们的对照组——因为这张图显示，虽然对照组更富有，但我们样本中较贫穷的人的收入变化平均来说也有相似的趋势，至少在我们考虑的前几年是这样。”

当事件研究图出现在他们面前时，房间又恢复了安静。木卫三的代表伸出手，拿出全息图的副本，仔细查看。她用手指将其展开来，以便仔细检查每个点。

“我注意到你自我介绍说，你是泰坦理工学院奥利·阿什费尔特社会协同学教授索利斯博士。我对奥利·阿什费尔特（Orley Ashenfelter ）的唯一了解就是‘阿什费尔特的沉降（Ashenfelter's Dip）’这句话。你能告诉我们关于这个问题以及它在这里出现的情况吗？”代表问道。

“Orley Ashenfelter 是 20 世纪末和 21 世纪初地球上的经济学教授。在他早期研究项目，他发现参加扶贫项目的人在参加之前的收入通常会相对于对照组下降。这种现象被称为“Ashenfelter's Dip”，通常表示与 Y(0) 未来轨迹相关的不可观测因素的组成不再与处理组可比。当人们参加这样的项目时，就会出现这个问题，通常是因为他们生活中与我们正在研究的结果相关的不可观测因素开始在他们周围崩溃。当他们的 Y(0) 时间路径发生这种结构性动态转变时，很难证明使用自愿非参与者作为对照是合理的。

但正如您所说，我们在这里看不到这一点。我们的处理前趋势大多是平缓的。我认为我们没有看到这一点，因为这种情况与奥利·阿什费尔特关注的情况不同。我们样本中的工人并没有选择他们的收入来参与。该计划使用了前一年的收入，在此之前他们没有时间调整收入，即使他们想这样做，而且选择的阈值也没有公布。回想一下Sant‘ Anna的DID第四步：“谁决定处理？他们知道什么？”处理并不是由工人决定的；而是由设定参与资格的政府决定的。即使第 26 分位数的人想参加，他们也不能加入。因此，我们不仅认为平行趋势假设在这种情况下是合理的，而且该计划的性质也机械地保证了无预期假设也可能是合理的。”

另一些参会代表继续问道，“您是否考虑过使用RD设计而不是DID？”。

索利斯博士回答说：“我们考虑过，我们也进行了那项分析。然而，考察该计划对所有参与者的总体影响，而不仅仅是那些处于临界点的参与者LATE，前者更具有政政策重要性。DID使我们能够得到处理胡整个处理人群的效应，而不仅仅是那些处于临界点周围的人群。这种方法为我们关注的人群提供了处理效应的无偏估计，捕捉了所有参与者之间的差异。”

索利斯博士可以看到，代表们很满意，并准备好展示主要结果，尽管事件研究全息图已经提供了大量信息。“现在，我们展示我们的发现，”她说，像发牌一样将一张张写有主要结果的全息卡翻给每位代表。“该计划的ATT是将平均收增加了约 449 scripts。然而，这种影响是两年间效应的平均值。在他们接受援助的那一年，收入相对于反事实增加了约 300 scripts。在 2292 年，收入比反事实增加了 600 scripts。因此，我们发现 ATT 总体为 449 scripts，这只是这两个年度效应的平均值（900/2 = 450）。”

来自恩塞拉多斯的某人举手。“这很有趣。但我仍然感到困惑。如果平行趋势是一个与两组均未接受处理群体的收入反事实趋势相关的假设，并且该项目根据基线收入选择人们进入处理组，导致处理组和对照组的基线收入出现巨大差距，那么为什么这不违反平行趋势？我想我会认为，对基于Y(0) 的任何处理配置都是一个问题。”

索利斯博士花了一点时间思考问题，然后回答。“我同意。我想在我阅读有关选择和平行趋势的文章之前，这些观点对我来说也是不容易理解的。我要说的第一点是，我认为我们永远不会知道平行趋势是否成立。正如 Kahn-Lang 和 Lang (2019) 指出的那样，即使是事件研究也无法告诉我们这一点。处理前趋势只能对现实说话——它们无法告诉我们反事实。因此，在评估具有观测研究设计的项目时需要谦逊——甚至可以说更加谦逊，因为在随机化中，我们对平行趋势的信心远不如我们对随机化将所有事物均匀分布在所有组中的信心。”

“我认为我们可以说，在我们的案例中，选择进入处理是机械的和确定性的。在我们的案例中，没有人主动选择进入这个项目。在很多方面，我们可以轻松地使用 RDD 来研究这个问题，因为基线时的收入是一个驱动变量。这是解释这一原则的一种方式，即只要断点周围两组的平均 Y(0) 决定因素的组成相似，基于基线 Y(0) 的选择处理就不会使得破坏平行趋势。这就是我们检验基线特征与处理相关性的原因。如果我们发现选择进入处理系统地选择了 Y(0) 趋势较高或较低的群体，那么我们需要找到解决方案，但在我们的案例中，这并没有发生。但它可能会发生。”

* 上述所有结果的stata code
* Step4a.do - Heterogeneous trends in Y0 with worker characteristics
clear all
set seed 2

* Set directory
cd "/Users/scunning/Library/CloudStorage/Dropbox-MixtapeConsulting/scott cunningham/0.1 Mixtape Consulting/Substack/Code"

* First create the colonies
quietly set obs 5
gen colony = _n

* Generate 250 workers in each colony
expand 250
bysort colony: gen unit_fe = runiform(1, 1000)
label variable unit_fe "Unique worker fixed effect per colony"
egen id = group(colony unit_fe)

* Generate group variable for different trends
gen group = mod(_n, 40) + 1

* Generate worker characteristics
* Age: skewed distribution with average age 35
gen age = round(rbeta(2, 5) * 65)

* Sex: balanced ratio
gen male = runiform() > 0.5

* Race: specified proportions
gen race = .
replace race = 1 if runiform() < 0.30
replace race = 2 if runiform() >= 0.30 & runiform() < 0.60
replace race = 3 if runiform() >= 0.60 & runiform() < 0.85
replace race = 4 if runiform() >= 0.85 & runiform() < 0.95
replace race = 5 if runiform() >= 0.95

label define race_lbl 1 "Black" 2 "White" 3 "Asian" 4 "Hispanic" 5 "Other"
label values race race_lbl

gen black = race==1
gen white = race==2
gen asian = race==3
gen hispanic = race==4
gen other = race==5

* Generate potential outcomes for the first year (1987)
gen y0_1987 = unit_fe + rnormal(0, 10)

* Generate potential outcomes for subsequent years with more stochastic elements
gen trend = 5
gen y0_1988 = y0_1987 + trend + rnormal(0, 10)
gen y0_1989 = y0_1988 + trend + rnormal(0, 10)
gen y0_1990 = y0_1989 + trend + rnormal(0, 10)
gen y0_1991 = y0_1990 + trend + rnormal(0, 10)
gen y0_1992 = y0_1991 + trend + rnormal(0, 10)

* Reshape data to long format
reshape long y0_, i(id) j(year)
rename y0_ y0

* Determine treatment status in 1990 based on Y0
su y0 if year == 1990, detail
gen treat = 0
replace treat = 1 if y0 <= `r(p25)' & year == 1990

* Ensure treatment status remains consistent across years
bysort id: egen max_treat = max(treat)
bysort id: replace treat = max_treat

* Post-treatment variable
gen post = 0
replace post = 1 if year >= 1991

* Apply group-specific trends in 1991 and 1992
replace y0 = y0 + 10 if year >= 1991 & group == 2
replace y0 = y0 + 20 if year >= 1991 & group == 3
replace y0 = y0 + 30 if year >= 1991 & group == 4

* Generate y1 by adding treatment effect for treated units
gen y1 = y0
replace y1 = y0 + 300 if year == 1991 & treat == 1
replace y1 = y0 + 600 if year == 1992 & treat == 1

* Treatment effect
gen delta = y1 - y0
label var delta "Treatment effect for unit i (unobservable in the real world)"

sum delta if post == 1, meanonly
gen ate = `r(mean)'
sum delta if treat == 1 & post == 1, meanonly
gen att = `r(mean)'

* Generate observed outcome based on treatment assignment
gen earnings = y0
qui replace earnings = y1 if post == 1 & treat == 1

* Show difference in earnings at baseline between treatment and control
bysort treat: sum earnings if year == 1990 // Check baseline differences between groups

* Illustrate parallel trends assumption
su y0 if treat == 1 & post == 0
gen ey0_10 = `r(mean)'
su y0 if treat == 1 & post == 1
gen ey0_11 = `r(mean)'
su y0 if treat == 0 & post == 0
gen ey0_00 = `r(mean)'
su y0 if treat == 0 & post == 1
gen ey0_01 = `r(mean)'

gen parallel_trends = (ey0_11 - ey0_10) - (ey0_01 - ey0_00)

reg y0 post##treat, robust

su parallel_trends

* Diff-in-diff
su earnings if treat == 1 & post == 0
gen ey_10 = `r(mean)'
su earnings if treat == 1 & post == 1
gen ey_11 = `r(mean)'
su earnings if treat == 0 & post == 0
gen ey_00 = `r(mean)'
su earnings if treat == 0 & post == 1
gen ey_01 = `r(mean)'

gen did = (ey_11 - ey_10) - (ey_01 - ey_00)

reg earnings post##treat, robust

su did att

* Run regressions and store results
eststo clear
foreach var in age white black asian hispanic other male {
    eststo: reg `var' i.treat i.colony, cluster(colony)
    
    * Store sample mean
    sum `var'
    estadd scalar mean = r(mean)
}

* Export results to LaTeX
esttab using "./reg_results.tex", ///
    keep(1.treat) ///
    coeflabels(1.treat "Treatment") ///
    cells("b(star fmt(2))" "se(par fmt(2))") ///
    star(* 0.10 ** 0.05 *** 0.01) ///
    label ///
    stats(N mean, fmt(%9.0fc %9.2fc) labels("Observations" "Sample Mean")) ///
    booktabs ///
    replace ///
    nonumbers ///
    noobs ///
    collabels(none) ///
    note("* p<0.10, ** p<0.05, *** p<0.01")
	
	
* Run regressions
eststo clear
eststo: reg earnings post##treat, robust cluster(colony)
estadd local att_label "Estimated ATT from 1991 to 1992"

eststo: reg earnings treat##ib1990.year, cluster(colony)
estadd local att_1991 "ATT in 1991"
estadd local att_1992 "ATT in 1992"

* Calculate pre-treatment mean
sum earnings if year < 1991
local pre_mean = r(mean)

* Export results to LaTeX
esttab using "main_results.tex", ///
    cells("b(star fmt(2))" "se(par fmt(2))") ///
    keep(1.post#1.treat 1.treat#1991.year 1.treat#1992.year) ///
    order(1.post#1.treat 1.treat#1991.year 1.treat#1992.year) ///
    coeflabels(1.post#1.treat "Estimated ATT from 1991 to 1992" ///
               1.treat#1991.year "ATT in 1991" ///
               1.treat#1992.year "ATT in 1992") ///
    stats(N pre_mean, fmt(%9.0fc %9.2fc) labels("Observations" "Pre-treatment mean")) ///
    star(* 0.10 ** 0.05 *** 0.01) ///
    label booktabs ///
    replace ///
    title("Main Results: Estimated total effects and effects by year of EquiLift on Trainees Subsequent Earnings") ///
    addnotes("* p<0.10, ** p<0.05, *** p<0.01") ///
    collabels(none) ///
    scalar(pre_mean)
	
* Event study regression
reg earnings treat##ib1990.year, cluster(colony)

* Generate the coefficient plot
coefplot, keep(1.treat#*) omitted baselevels cirecast(rcap) ///
    rename(1.treat#([0-9]+).year = \1, regex) at(_coef) ///
    yline(0, lp(solid)) xline(1990.5, lpattern(dash)) ///
    xlab(1987(1)1992) ///
    title("Estimated Effect of EquiLift on Earnings") ///
    subtitle("Event Study Analysis")
	
graph export ./equilift_es.png, as(png) replace

local planet_names "Mars Europa Titan Enceladus Ganymede"

* Initialize plot index
local i = 1

* Create individual kernel density plots for each colony with a reference line
foreach planet of local planet_names {
    * Calculate the 25th percentile for earnings for this specific colony at baseline
    sum earnings if year == 1990 & colony == `i', detail
    local p25 = r(p25)
    
    * Get the maximum density value
    kdensity earnings if year == 1990 & colony == `i', nograph generate(x d)
    sum d
    local max_density = r(max)
    drop x d
    
    twoway (kdensity earnings if year == 1990 & colony == `i', lcolor(blue) lwidth(medium) ///
            ylabel(none) ytitle("")) ///
           (pci 0 `p25' `max_density' `p25', lcolor(red) lpattern(dash)), ///
           title("`planet' Colony")  ///
           xtitle("Earnings") ///
           legend(off) ///
           name(plot_`i', replace)
    local ++i
}

* Combine the individual plots into a single graph
graph combine plot_1 plot_2 plot_3 plot_4 plot_5, ///
    title("Distribution of Worker Earnings at Baseline") subtitle("Enrollment cutoff at 25th percentile")
	
* Export the combined graph as a PNG file
graph export "./baseline_earnings_combined.png", as(png) replace


* Keep only the year 1990 data
keep if year == 1990

* Calculate mean and standard deviation for earnings in 1990 by treatment status
egen mean_earnings = mean(earnings), by(treat)
egen sd_earnings = sd(earnings), by(treat)

* Generate variables for whiskers
gen lower_whisker = mean_earnings - sd_earnings
gen upper_whisker = mean_earnings + sd_earnings

* Calculate mean earnings for labels
summarize earnings if treat == 0
local mean_control = r(mean)
summarize earnings if treat == 1
local mean_treatment = r(mean)

* Create variables for label positions
gen label_x = treat
replace label_x = label_x + 0.2 if treat == 0
replace label_x = label_x - 0.2 if treat == 1

* Generate labels for the means
gen mean_label = ""
replace mean_label = string(mean_earnings, "%9.0f") + " Scripts"

* Adjust y-coordinates for label positioning
gen label_y = upper_whisker + 50

* Create the bar chart with different colors and add whiskers for standard deviation
twoway (bar mean_earnings treat, barwidth(0.75) color(gs12 gs14) lcolor(none)) ///
       (rcap lower_whisker upper_whisker treat, lcolor(black)) ///
       (scatter mean_earnings treat, mcolor(black) msymbol(diamond) msize(medium)) ///
       (scatter label_y label_x, mlabel(mean_label) mlabpos(0) mlabcolor(black) mlabsize(large) msymbol(none)), ///
       ytitle("Baseline Earnings") ///
       title("Mean Earnings at Baseline") subtitle("Treatment and Control") ///
       xtitle("Treatment Status") ///
       xlabel(0 "Control" 1 "Treatment") ///
       legend(off)

* Save the graph as a PNG file
graph export ./baseline_earnings.png, as(png) replace

* Restore the original dataset
restore