最近几年随着微信公众号和知乎等迅速成为各类信息传播的主要平台，国内科研工作者慢慢的很少”逛”论坛、发帖子、写博客了，有的转战公众号开设自己的账号，有的开始团队运营，这些平台给科学研究提供了不少便利性，扩大了知识获取面，降低了知识检索的时间成本。但是其弊端也毕竟明显，一方面是导致很多人对前沿科学咨讯的获取基本来源于公众号，很少阅读原始文献，另一方面是公众号或者知乎平台分享的信息和知识点没有经过同行评议，其实不乏有错误之处，若不加思索的参考或者照搬很容易导致研究方法误用，特别是对于尚处于学术生涯初期的研究生。
关于介绍因果推断方法的书籍非常多了，包括Angrist & JoshuaDavid著的《Mostly Harmless Econometrics: An Empiricist Companion》，Wooldridge的《计量经济学导论》等等。本书并不打算再次搬运以上这些著作的内容，而是主要记录一个学习者在学习和使用这些方法过程中的心路历程，分享学习和使用这些方法的一些经验和体会。
本章节将对四种常见的因果推断方法进行介绍：
本书的主要参考资料如下：
- Angrist, & JoshuaDavid. (2010). Mostly Harmless Econometrics: An Empiricist Companion. Princeton University Press.
- 伍德里奇, & 费剑平校. (2010). 计量经济学导论: 第4版. 中国人民大学出版社.
- Scott Cunningham. (2020). Causal Inference: The Mixtape. Yale University Press.
- https://lost-stats.github.io/Model_Estimation/Research_Design/event_study.html
双重差分法（Difference-in-Difference，DID）通常又可称作倍差法或差异中的差异法，是目前社会科学领域用于开展政策效果评估和识别两个因素因果关联的最常用方法之一。DID最早由Orley Ashenfelter和David Card提出，两位学者于1985年在学术期刊《Review of Economics and Statistics》上发表了运用DID的方法估计美国综合就业和培训法案对收入的影响 Ashenfelter and Card (1985) 。DID本质是一种准实验设计（quasi-experimental design），利用干预组和对照组的纵向数据来获得适当的反事实（counterfactual）来估计因果效应，具体而言就是比较干预组人群和对照组人群之间的结果随时间的变化，其通常用于估计特定的干预或治疗（如某项法律的通过、某种政策的实施或大规模项目实行）的效果。
DID的基本原理很简单，但是如果据此就认为采用DID识别框架的方法都很简单就错了，随着Joshua D. Angrist、Guido W. Imbens、Gary King、Susan Athey等经济学和统计学家对因果识别方法体系的贡献，时至今日，DID已经衍生出了非常多的类型来适用不同的研究设计和数据结构，可以说已经是一大方法家族了吧。
DID的统计学基础是哈佛大学著名统计学家 (Rubin 1980) 提出的潜在结果模型（Potential Outcomes Model）或称为反事实框架 (Counterfactual Framework)，其核心思想在于利用“反事实”假设构造一个可比的对照组来消除样本选择偏误（Selection Bias）。根据潜在结果模型可知，要想获得某种干预的因果效应就必须比较同一个研究对象在接受干预和不接受干预时的结果差异，此时这种差异才是接受干预相对于不接受干预的效果。然而，对于同一研究对象而言，通常无法既观察其接受干预的结果，又观察其不接受干预的结果。对于接受干预的研究对象而言，不接受干预时的状态是一种 “反事实” 状态，对于不接受干预的研究对象而言，接受干预时的状态也是一种 “反事实” 状态。正是由于在同一研究对象中同一时点接受干预与非干预的状态存在互斥性，因此无法通过直接观察“反事实”状态下的效应值。
在医学研究中，为了解决这个问题通常采用随机对照试验（RCT）的方法，通过随机化的方法将一组具有高度同质性的研究对象随机分配去接受干预（Trentment Group）或不接受干预（Control Group）。由于研究对象具有高度同质性且通过随机化消除了样本选择偏误，那么就可以在限定条件下认为，未接受干预的那部分研究对象（Control Group）就是接受了干预的研究对象（Trentment Group）在未接受干预时（Untreated）的状态。但是，在政策研究中，通常无法人为对政策的执行施加影响，无法做到随机化，干预组和对照组的研究对象很难保证具有高度同质性，尽管如此，计量经济学家们探索出了利用干预组与对照组之间可能存在不随时间变化的系统性差异的特点，构造出了消除样本选择偏误而得到相对纯净政策效应的计量方法，其中就以DID为代表。
双重差分的基准模型如下：
\[ y_{it} = \beta_1 \cdot Treat*Time + \beta_2 Treat + \beta_3 Time + \beta_0+ \epsilon_{it} \]
各参数代表意义就不解释了，其采用ols估计就行了，该模型适用于前面所提到的第一种和第二种中的两期的数据类型。
但是，如果仔细看的话，该模型就是采用的双向固定效应模型，其中treat是个体效应，time是时间效应。因此，对于多期的面板数据，你又会看见下面的模型：
\[ y_{it} = \beta \cdot Treat*Time + \lambda t + \mu_i + \epsilon_{it} \]
其中没有了treat和time，取而代之是lambda和mu代表的个体效应和时间效应，因为在多期面板中，如果仍然用treat和time会过于粗糙，time只有0和1，无法识别出完整的时间效应，而treat同样存在这个问题。
当然模型还有其他变化，比如当treat为连续变量时，政策冲击不是在同一时间单一完成的，但是归根结底都是双向固定效应模型的变换形式。
为什么通过两次差分就能估计出净效应的原理不详述，这里只讨论一个问题：为什么RCT的研究在估计干预效应时只diff一次，而这里要diff两次？
这就涉及到了双重差分方法的一个最重要的前提假设：共同趋势（Common Trends）假设，即干预组和对照组在政策实施之前必须具有相同的发展趋势。
那么为什么要有这个假设呢，可以看下面的示意图：
左边的图就是最常见的双重差分原理的示意图，共同趋势假设的目的就是为了保证基线期（before）的diff，如果没有treatment，干预组和对照组之间的差异仍然是这个diff，所以在第二次差分时扣掉这个diff就能得到净干预效应。
所以这个共同趋势假设是为了保证干预和对照组在干预前后基础的变化都是一致的，它包含两层意思：一是必须趋势不变，也就是都是上升或下降或不变；二是必须平行。因此这个共同趋势假设只是一个基础，是在这个假设成立的基础上，才能推断在干预后，干预组和对照组也会（是也会，不是一定，所以这就是双重差分并不完美的地方）按照这个趋势变化，才能利用第二次差分。而有的公众号中写的是只要共同趋势满足，就可以大胆使用双重差分，其实并不是。
上面只说了共同趋势假设是为了保证反事实情况下的变化一致，但是并没有说出有这个假设的原因。应该也有人会问，为什么所有的双重差分的示意图在干预前的两条线都是平行的，但是有间隙。原因就是前面所讲的，这种自然试验无法保证干预组和对照组是随机分配的，换句话说就是这两组人在没有施加干预时，就是两个人群，所以会有差异，也就是示意的间隙。
因为这个“间隙”（就是组间差异）比较大，所以就需要在干预后把它扣除掉，如果不扣除直接对干预后的两组进行差值比较，那么这个差值中时包含了这个“间隙”的，如果这个“间隙”比真实干预效应还大，那么得到的两组差值就失去意义了，所以才会有第二次diff（第一次diff两组分别是前后减，第二次diff是在第一次的结果上干预和对照减）。
那么，如果这个“间隙”足够小呢，小到随机误差那么大，是不是就可以不用减的。那这就是为什么RCT只diff一次的原因，因为在试验设计阶段，通过严格的随机化已经保证的干预组和对照组个体的一致性，也就是已经将这个“间隙”降到足够小了，就像右边的图示意的一样。
DID对数据有比较高的要求，通常有四种情况下可以考虑双重差分。
两个时间点的个体横断面数据（微观数据）：
两个横断面个体不是同一批人，也就是常说的两个混合截面，很多较为早期开展的survey数据是这种形式。只要能找到一个变量，能够在第一个横断面中识别出与第二个横断面中受到干预的个体相对应的个体即可，也就需要找到第一个横断面中哪些个体属于干预组。通常可以通过地区变量来分辨，对于试点政策多采取这种方法，因为追踪调查比较费时费力。缺点：无法避免两次横断面调查带来的抽样误差，因为两个横断面的个体并不是同一批人，虽然抽样误差会在组间diff的时候被减去，但这种情况需要抽样误差在对照和干预个体中是均匀分布的。
两期及更多期的追踪数据，本质就是面板数据：
和上面情况一样，都是survey数据，但是为追踪调查，也就是每个时间点调查的个体都是同一批人。这种数据最优，可以进行的差分方法变化也最多，但是较难获得，好在最近些年开展的大规模人群调查都在采取追踪调查的方式。
两期及更多期的机构或地区数据，就是常见的机构或地区面板数据（宏观数据）：
比第二种情况更容易获得数据，但是缺点在于很难收集到样本量足够大、分析单元足够小的数据，比如较为容易获得的省级面板数据实际缺点较大。
单一的横断面数据：严格来讲是无法进行双重差分设计的，但是也有采用横断面进行差分研究的实例，比如研究中国大饥荒的通过不同出生队列构建的差分模型。
正如前面所将，双重差分法并不要求干预组和对照组是完全一致的，可以存在一定的差异，但是要求这种差异不随着时间产生变化，这也是在共同趋势假设成立的基础上，能够进一步推断干预后这个差异仍然是基线的差异的基础，因为前面几年是平行变化的，并不代表在政策施加之后也是平行变化的。但是这个不随时间变化的要求是无法检验的，所以这就成了双重差分法的一个无法避免的weakness。
那么，有人会问如果共同趋势不满足怎么办，当然会有其他的方法：匹配+双重差分，最常见的是psm+did，当然也可以是cem+did
对于两期的混合截面或者面板数据而言，是无法进行共同趋势检验的，所以有的论文中在这种情况下就不提共同趋势的事，对于多期面板虽然可以检验共同趋势，但是也会存在这个假设不满足的情况。
所以，为了研究的严谨性和可继续下去，有的学者就提出提出对样本进行匹配，因为共同趋势不满足的根本原因就是干预组和对照组之间差异太大，所以就通过匹配的办法让两组个体更加可比一些，将“间隙”缩小一些。
那么有的人又会问，既然前面说的横断面数据可以用匹配的方法进行效应估计，那既然都匹配了，也就是说趋近随机分组了，那么为什么还要diff两次。解释就是想要真正做到随机分组，匹配的变量需要很多，匹配的效果与RCT还是有差别，两组间的“间隙”仍然存在，所以还是需要用双重差分估计更干净的效应。
这里需要注意的一点是，应用匹配+双重差分的方法，不同的数据情况采用不同的匹配策略。
断点回归（Regression Discontinuity）适用于以下情形：人群是否接受干预（Treatment）是依据某一数值变量（rating variable）是否高于或低于某一确定的阈值（threshold）或者分割点（cut-point），例如在研究是否上大学会影响收入时，数值变量（rating variable，也叫 assignment variable，score，running variable，forcing variable， or index）就是高考分数，阈值或者分割点就是本科录取分数线。以下内容主要参考：(Calonico et al. 2022; D. Cattaneo, Idroboy, and Titiunik 2017; Jacob, Zhu, and Marie-Andrée 2012; 谢谦, 薛仙玲, and 付明卫 2019)。
Regression Discontinuity最早由社会学家Thistlethwaite and Campbell 在1960年提出的，用于评估社会项目，但是他们的研究虽然引起一些影响但是没有得到广泛的注意，后来，在1972-2009年期间，被一系列经济学家(Goldberger，van der Klaauw，Imbens and Kalyanaraman等)在方法学方面进行了完善，最终在2008年达到顶峰，标志是在2008年Journal of Econometrics出了一期RD分析的Special Issue，光看期刊名字就知道是经济学顶刊了(Journal of Econometrics是公认的计量经济学顶尖期刊，是教育部认可的12本经济学国际顶级期刊之一 )。
RD有两个特点：
RD有两种类型：
精确断点(sharp design): 即所有个体(All subject)在明确的cut-point之后全部接受干预(treatment)
模糊断点(fuzzy design): 即在cut-point前后，存在 no-shows (treatment group members who do not receive the treatment)或者crossovers (control group members who do receive the treatment) ，换句话说就无法找到一个明确的cut-point完全区分干预和对照组。更严格的分法是，fuzzy也可以分成两类：type I是no-shows，只存在处理组有未接受处理的个体，type II是同时存在no-shows和crossovers。实际就是RCT中常说的沾染问题。
由于RD仍然属于非实验方法，尽管也被成为类实验(quasi-experimental)，但本质还是非实验方法(nonexperimental)， 所以它必须满足一系列前期条件，才能提供无偏估计和更可能的接近RCT的严格情形。
一、 Rating variable (中文的翻译很多，但比较通用的翻译为驱动变量或者配置变量)不能被干预(treatment)所影响，Rating variable 的值必须是在干预(treatment)之前就已经确定了，或者是不可再更改的变量。也就是只能是驱动变量的值决定是否接受干预，不能是干预决定驱动变量的值。比如：高考成绩(Rating variable)是在判断是否能够上大学(treatment)之前就已经确定了，并且分数不会再发生变化，如果存在严重的根据分数线修改高考成绩的情况，则高考分数就不能作为Rating variable。 检验方法是：画Rating variable的频率分布直方图或者密度图，再进行McCrary 检验
二、 断点(cut-point)是完全外生的，并且干预与否完全取决于驱动变量和断点。即断点不受其他因素影响而改变干预的判断，即把本不该接受干预的对象划入了干预，或者把该接受干预的对象划入了对照。 比如，如果高考录取线(cut-point)的确定是完全已经严格的录取指标划定的， 而不存在为了使得某部分人群能够获得上大学的机会，而修改录取线。
三、 在驱动变量的cut-point前后，除了干预与否(treatment)是不连续或者跳跃的(比如断点前为0，断点后为1)，其他任何因素都不能出现不连续或者跳跃(Nothing other than treatment status is discontinuous in the analysis interval)，因为要保证outcome的跳跃只能是干预的不连续或者跳跃导致的。比如，如果高考分数线除了决定是否能上大学外，还决定是否享有创业扶持机会，那么如果研究上大学对某种结局的影响，就无法分离出到底是大学教育导致，还是获得创业扶持机会导致的。
四、 驱动变量(rating variable)与结局变量(outcome)之间的函数关系应该是连续的，严格的说，应该是在进行RD分析的这段区间内(interval)是连续的。换句话说就是，如果不存在干预的情况，驱动变量与结局变量之间不会出现不连续或者跳跃，才能在结局变量出现不连续的情况下，反推出是干预引起的，因为只有干预是不连续的。注意：此条件只需要在选择参数估计方法是要满足(applies only to parametric estimators)
图形分析是进行RD的第一步，也是非常重要的组成部分
通常，进行RD至少需要画4种图：
示列软件: R version 4.0.5 (2021-03-31)
示列分析包：在R语言中一共有三个package可以进行RD分析：
rdrobust简介：目前该包的开发项目可在https://sites.google.com/site/rdpackages/home主页上查看，该包共有三个函数： + rdplot：用来画图 + rdrobust：用来进行局部非参数估计 + rdbwselect：用来进行bandwidth选择 另外还有一个协同的包rddensity，需要从CRAN上单独安装，用来进行对驱动变量是否被操控进行检验，因为rdrobust包未包含McCrary检验
示列数据集：采用rdrobust包中自带的rdrobust_RDsenate数据集：
根据研究设计，判断是否存在采用RD进行因果识别或者效应估计的可能，即是否可以找到合适的驱动变量和明确的断点来识别施加干预的状态。
通常情况下时间（具体到天）、年龄（如研究退休，代表性的文章: “退休影响会健康吗”）、地理距离（研究雾霾，代表性文章: “冬季供暖导致雾霾？来自华北 城市面板的证据”）比较容易作为驱动变量。
1) 条件一：通过画驱动变量的频率分布直方图或者密度分布图，以及McCrary检验
可以利用graphics包或者ggplot2包画histogram图，也可以利用rddenstiy包中的rddenstiy函数画密度图，也可利用rdd包进行McCrary检验。
rddenstiy函数的其他参数可以用默认值，输出结果中，最下面一行robust的p值大于0.05，就说明驱动变量未受到操纵。
library(rdrobust)
library(rddensity)
library(magrittr)
data(rdrobust_RDsenate)
hist(x = rdrobust_RDsenate$margin,
breaks = 30,
col = 'red',
xlab = 'margin',
ylab = 'Frequance',
main = 'Histogram of Rating Variable')
abline(v = 0)
library(ggplot2)
ggplot(rdrobust_RDsenate) +
geom_histogram(aes(x = margin,
fill = margin < 0),
color = 'black',
binwidth = 4, boundary = 4) +
geom_vline(xintercept = 0, size = 1) +
labs(title = 'Histogram of the Rating varibale',
y = 'Number of Observation',
x = 'Margin') +
scale_fill_manual(name = '',
values = c('red', 'blue'),
labels = c('Control', 'Treat')) +
theme_bw() +
theme(legend.position = c(.9, .9),
legend.background = element_blank())
rdplotdensity(rddensity(rdrobust_RDsenate$margin),
X = rdrobust_RDsenate$margin)$Estplot[1]$data
list()
attr(,"class")
[1] "waiver"
2) 条件二三四根据研究设计及背景资料进行判断。
画驱动变量与干预的散点图，判断是否为sharp或者fuzzy类型。可以看出本例中，干预在cut-point前后，有明确的界限，为sharp类型。
rdrobust_RDsenate$treatment <-
ifelse(rdrobust_RDsenate$margin < 0, 0, 1)
rdrobust_RDsenate$color <-
ifelse(rdrobust_RDsenate$margin < 0,
'blue',
'red')
plot(x = rdrobust_RDsenate$margin,
y = rdrobust_RDsenate$treatment,
col = rdrobust_RDsenate$color,
type = 'p',
pch = 16,
xlab = 'Margin',
ylab = 'Treatment',
main = 'Relationship between
rating variable and treatment')
abline(v = 0)
一、画散点图
横轴为驱动变量，纵轴为结果变量，画散点图，初步判断两者之间关系，但是散点图噪音(noise)太大，不利于发现跳跃点，需要对驱动变量分不同的段(intervals or bin)将使得关系取线更平顺。
# 第一步：散点图
plot(x = rdrobust_RDsenate$margin,
y = rdrobust_RDsenate$vote,
type = 'p',
col = rdrobust_RDsenate$color,
pch = 16,
cex = 0.8,
xlab = 'Margin',
ylab = 'Vote')
abline(v = 0)
二、将驱动变量分段的步骤
将驱动变量分段后再做关系图，大致有四步：
三、确定bins数量的方法
如何设定bin的数量或者binwidth的方法比较复杂，bin的数量不宜过多或过少，过少起不到将噪音的作用，过多会导致jump不明显，在Stata和R中有打包的好的方法去判断bin的数量，本例中采用rdrobust包中的rdplot函数进行判断。
rdplot函数的常用参数有：
判断方法常用一共4种：Evenly-spaced (es)，Quantile-spaced (qs)，以及结合mimicking-variance (MV)的esmv和qsmv，推荐使用qs方法。p和nbin如果不确定也可缺失，让binselect方法自行设定，从结果可以看出总的样本量，以及cut-point左右的样本量，以及设定的bin的数量。
bins <- rdplot(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 0, p = 4,
nbins = c(20, 20),
binselect = 'esmv',
kernel = 'uniform')
summary(bins)Call: rdplot
Number of Obs. 1297
Kernel Uniform
Number of Obs. 595 702
Eff. Number of Obs. 595 702
Order poly. fit (p) 4 4
BW poly. fit (h) 100.000 100.000
Number of bins scale 1 1
Bins Selected 20 20
Average Bin Length 5.000 5.000
Median Bin Length 5.000 5.000
IMSE-optimal bins 8 9
Mimicking Variance bins 15 35
Relative to IMSE-optimal:
Implied scale 2.500 2.222
WIMSE variance weight 0.060 0.084
WIMSE bias weight 0.940 0.916
1) 效应大小的估计方法一共有两种：
\[ \begin{aligned} Y_i = & \alpha + \beta_0 \cdot Treatment_i + \beta_1 \cdot (rating_i - cutpoint)^k + \\ & \beta_3 \cdot Treatment_i \cdot (rating_i - cutpoint)^k + \\ & \beta_4 \cdot Confouders_i + \epsilon_i \end{aligned} \]
关于上述两种方法的比较或者权衡，主要是关于精度和误差，全局参数估计因为利用了全部样本，通常对效应的估计精度更高，然而其也存在一个问题，就是驱动变量和结局变量的关系函数在越大的区间下越难以准确识别；而非参数估计可以通过局部拟合加权线性或者多项式回归模型而减少这种偏误。
以上方法主要针对sharp类型，fuzzy类型的估计方法见最后一节。
2) 全局参数估计
分为一次线性、一次线性加交互项、以及二次、三次线性及交互项，共6个模型，比较AIC或者回归残差，较小者模型较优。
rdrobust_RDsenate$margin_del <- rdrobust_RDsenate$margin - 0
# linear
fit_1 <- lm(vote ~ margin_del + treatment,
data = rdrobust_RDsenate)
# linear interaction
fit_2 <- lm(vote ~ margin_del * treatment,
data = rdrobust_RDsenate)
# quadratic
fit_3 <- lm(vote ~ margin_del +
I(margin_del ^ 2) + treatment,
data = rdrobust_RDsenate)
# quadratic interaction
fit_4 <- lm(vote ~ (margin_del + I(margin_del ^ 2)) *
treatment,
data = rdrobust_RDsenate)
# cubic
fit_5 <- lm(vote ~ margin_del + I(margin_del ^ 2) +
I(margin_del ^ 3) + treatment,
data = rdrobust_RDsenate)
# cubic interaction
fit_6 <- lm(vote ~ (margin_del + I(margin_del ^ 2) +
I(margin_del ^ 3)) * treatment,
data = rdrobust_RDsenate) library(modelsummary)
modelsummary(list("Linear" = fit_1,
"Linear interaction" = fit_2,
"Quadratic" = fit_3),
output = "kableExtra",
stars = TRUE)| Linear | Linear interaction | Quadratic | |
|---|---|---|---|
| (Intercept) | 47.331*** | 44.904*** | 45.486*** |
| (0.542) | (0.699) | (0.616) | |
| margin_del | 0.348*** | 0.216*** | 0.285*** |
| (0.013) | (0.028) | (0.017) | |
| treatment | 4.785*** | 6.044*** | 6.663*** |
| (0.923) | (0.942) | (0.963) | |
| margin_del × treatment | 0.170*** | ||
| (0.032) | |||
| I(margin_del^2) | 0.001*** | ||
| (0.0002) | |||
| Num.Obs. | 1297 | 1297 | 1297 |
| R2 | 0.578 | 0.587 | 0.590 |
| R2 Adj. | 0.577 | 0.586 | 0.589 |
| AIC | 10083.7 | 10056.7 | 10049.9 |
| BIC | 10104.4 | 10082.5 | 10075.7 |
| Log.Lik. | −5037.848 | −5023.354 | −5019.946 |
| F | 886.433 | 613.586 | 619.091 |
| RMSE | 11.77 | 11.64 | 11.61 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
3) 局部非参数估计
利用rdrobust包的rdrobust函数进行拟合。由于该包是由Sebastian Calonico, Matias D. Cattaneo and Rocío Titiunik开发，所以该包中所提供的带宽选择方法可认为就是CCT法，因为CCT本来就是Calonico, Cattaneo, and Titiunik三人姓氏的缩写。
有两种确定带宽的方法：
rdrobust包中一共提供了10种选择最优带宽的方法，如下:
其中，共包含MSE = Mean Square Error和CER = Coverage Error Rate两大类，MSE更适用于进行点估计的带宽选择， 而CER更适合区间估计的带宽选择。
以rd结尾是表示选择的带宽在cut-point左右相等，而以two结尾是表示选择的带宽在cut-point左右不相等。
h即为选择的左右最优带宽，而b给出的是用来进行敏感性分析时应该考虑的带宽。
rdbwselect(y = rdrobust_RDsenate$vote,
x = rdrobust_RDsenate$margin,
all = TRUE) %>% summary()Call: rdbwselect
Number of Obs. 1297
BW type All
Kernel Triangular
VCE method NN
Number of Obs. 595 702
Order est. (p) 1 1
Order bias (q) 2 2
Unique Obs. 595 665
=======================================================
BW est. (h) BW bias (b)
Left of c Right of c Left of c Right of c
=======================================================
mserd 17.754 17.754 28.028 28.028
msetwo 16.170 18.126 27.104 29.344
msesum 18.365 18.365 31.319 31.319
msecomb1 17.754 17.754 28.028 28.028
msecomb2 17.754 18.126 28.028 29.344
cerrd 12.407 12.407 28.028 28.028
certwo 11.299 12.667 27.104 29.344
cersum 12.834 12.834 31.319 31.319
cercomb1 12.407 12.407 28.028 28.028
cercomb2 12.407 12.667 28.028 29.344
=======================================================虽然在此文中rdrobust: An R Package for Robust Nonparametric Inference in Regression-Discontinuity Designs提到也可以通过bdselect = ’CV’或者bdselect = ’IK’来采用CV和IK法，但是在写这篇笔记时，已无法使用这两种方法。
输出的结果的最后一张表的第一行Conventional即为干预的处理效应。
# 如果有混杂因素需要控制,可以用covs = c('var1', 'var2')
loc_fit_1 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 0, p = 1,
kernel = 'triangular',
bwselect = 'msetwo')
# c用来指定cut-point，p用来指定局部加权回归的多项式幂次
loc_fit_2 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 0, p = 2,
kernel = 'triangular',
bwselect = 'msetwo')
loc_fit_3 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 0, p = 1,
kernel = 'triangular',
bwselect = 'cerrd')
loc_fit_4 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 0, p = 2,
kernel = 'triangular',
bwselect = 'certwo')modelsummary(list(loc_fit_1, loc_fit_2,
loc_fit_3, loc_fit_4),
output = "kableExtra",
stars = TRUE,
estimate = "std.error")| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| Conventional | 1.497*** | 1.852*** | 1.680*** | 2.209*** |
| (1.497) | (1.852) | (1.680) | (2.209) | |
| Bias-Corrected | 1.497*** | 1.852*** | 1.680*** | 2.209*** |
| (1.497) | (1.852) | (1.680) | (2.209) | |
| Robust | 1.759*** | 2.056*** | 1.841*** | 2.282*** |
| (1.759) | (2.056) | (1.841) | (2.282) | |
| Kernel | Triangular | Triangular | Triangular | Triangular |
| Bandwidth | msetwo | msetwo | cerrd | certwo |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
敏感性分析主要用来检验模型估计结果的稳健性，RD分析主要有四种敏感性分析方式;
对前面所述的Internal Validity条件三的检验，即除了干预(Treatment)变量外，其他的任何非结局变量(包括混杂因素在内)都不能在cut-point处出现不连续或称为跳跃，如果存在，则无法推断结局变量的跳跃是由干预(Treatment)引起的，检验方式就是将所有的混杂因素作为结局，利用rdrobust函数进行检验，输出结果的coef应该很小，且p值应该统计学不显著 (Continuity-Based Analysis for Covariates)。
对cut-point的敏感性分析，即是更换cut-point，检验是否左右还存在处理效应，如果更换断点后，仍然存在处理效应，则无法说明本研究的干预措施是有效的，因为在研究设定的断点处识别到的处理效应，有可能是由其他因素引起的。
对cut-point附近个体的敏感性分析，在Internal Validity条件一中提到，驱动变量不可被操控，但是这个无法直接检验，因此，换个思路，如果驱动变量被超控，那自然是在cut-point左右离得最近的值被操纵的可能性最大，所以如果将这部分个体剔除，若仍然能观测到处理效应，则说明这种效应是真实由干预所导致。这种方法被称为甜甜圈(donut hole)法，同样也适用于个体过多的堆积与cut-point附近的RD分析。
对带宽bandwidth的敏感性分析，即是cut-point不改变，而是更换选择的最优带宽，进行多次局部非参数检验，如果更换带宽之后，仍然能在断点处识别的处理效应，说明研究的干预措施是有效的，因为在研究设定的断点处识别到的处理效应不是由于某一特点的带宽下才观测到的，说明处理效应稳健。
1) 对cut-point的敏感性分析
选择一个虚拟的断点用来替换真正的断点，但是真实的干预情况不改变。关于虚拟断点如何选择以及选择多少个并没有明确的标准，通常在真实断点的附近，左右对称选取即可
# 除c以外的其他参数应该与前面分析保持一致
sen_cut_1 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = 1, p = 1,
kernel = 'triangular',
bwselect = 'msetwo')
# 除c以外的其他参数应该与前面分析保持一致
sen_cut_2 <- rdrobust(rdrobust_RDsenate$vote,
rdrobust_RDsenate$margin,
c = -1, p = 1,
kernel = 'triangular',
bwselect = 'msetwo') modelsummary(list(sen_cut_1, sen_cut_2),
output = "kableExtra",
stars = TRUE,
estimate = "std.error")| Model 1 | Model 2 | |
|---|---|---|
| Conventional | 1.406*** | 1.614 |
| (1.406) | (1.614) | |
| Bias-Corrected | 1.406*** | 1.614 |
| (1.406) | (1.614) | |
| Robust | 1.640** | 1.865 |
| (1.640) | (1.865) | |
| Kernel | Triangular | Triangular |
| Bandwidth | msetwo | msetwo |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
2) 对cut-point附近个体的敏感性分析
对称剔除真实断点左右一定较小范围内的个体。关于范围的大小同样没有明确的标准，通常选取多个范围比较即可。
rdrobust_RD_hole_1 <- subset(rdrobust_RDsenate,
abs(margin) > 0.3)
rdrobust_RD_hole_2 <- subset(rdrobust_RDsenate,
abs(margin) > 0.4)
rdrobust_RD_hole_3 <- subset(rdrobust_RDsenate,
abs(margin) > 0.5)
# 除c以外的其他参数应该与前面分析保持一致
sen_hole_1 <- rdrobust(rdrobust_RD_hole_1$vote,
rdrobust_RD_hole_1$margin,
c = 0, p = 1,
kernel = 'triangular',
bwselect = 'msetwo')
# 除c以外的其他参数应该与前面分析保持一致
sen_hole_2 <- rdrobust(rdrobust_RD_hole_2$vote,
rdrobust_RD_hole_2$margin,
c = 0, p = 1,
kernel = 'triangular',
bwselect = 'msetwo')
# 除c以外的其他参数应该与前面分析保持一致
sen_hole_3 <- rdrobust(rdrobust_RD_hole_3$vote,
rdrobust_RD_hole_3$margin,
c = 0, p = 1,
kernel = 'triangular',
bwselect = 'msetwo') modelsummary(list(sen_hole_1, sen_hole_2, sen_hole_3),
output = "kableExtra",
stars = TRUE,
estimate = "std.error")| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| Conventional | 1.558*** | 1.558*** | 1.591*** |
| (1.558) | (1.558) | (1.591) | |
| Bias-Corrected | 1.558*** | 1.558*** | 1.591*** |
| (1.558) | (1.558) | (1.591) | |
| Robust | 1.839*** | 1.860*** | 1.918*** |
| (1.839) | (1.860) | (1.918) | |
| Kernel | Triangular | Triangular | Triangular |
| Bandwidth | msetwo | msetwo | msetwo |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
3) 对带宽bandwidth的敏感性分析
即是更换不同的带宽即可，带宽的变化实际是cut-point区间范围尾部的个体在发生变化，可以通过选取不同的带宽选择方式，真实存在的处理效应不会随着带宽的改变而发生变化。
Fuzzy类型的RD可以根据驱动变量与干预的实际情况采取不同的处理方式:
如果只是no-shows的情况，可以结合研究目的，决定是否采取意向性分析法(intent-to-treat)
采用两阶段回归two-stage least squares (2SLS) ，两阶段模型设定如下:
第一阶段：
\[ Treatment_i = \alpha_1 + \gamma_0 \cdot treat-by-cutpoint_i + \epsilon_i \]
第二阶段：
\[ Y_i = \alpha_i + \beta_0 \cdot \dot Treatment + \mu_i \]
通过两阶段回归，第二阶段中的Standard errors是经过了调整的。在rdrobust包的rdrobust函数中，有一个fuzzy的逻辑参数，用来指定RD类型是否为fuzzy，可以很方便的进行fuzzy类型的模型估计。
由于真实世界的复杂性，对因果推断的分析方法的要求也越来越高，没有一种方法可以适用所有的情形，因此RD的方法也衍生发展出很多类型，包括多个断点(Multi-Cutoff RD Design)、多个驱动变量(Multi-Score RD Design)、离散驱动变量(Multi-Score RD Design)等等。