自然实验（natural experiment）
随着相关研究理论的不断深入，其分析方法也在不断拓展，目前常用的有双重差分模型（difference-in-differences，DID）、工具变量法（instrumental variables，IV）、间断时间序列模型（interrupted time series，ITS）、断点回归设计（regression discontinuity design，RDD）及其衍生的
方法等）
间断时间序列分析（ITS）
ITS应用
- 评估政策
- 时间范围
- 因变量
- 模型
- 结果
R语言复现
1程序包加载
library(readxl)##读入excel数据
library(car)##DWtest
library(orcutt)##Cochrane-Orcutt估计
library(prais)## prais_winsten估计
library(sandwich)##NeweyWest
library(lmtest)#线性回归测试
library(ggplot2)##绘图
1 数据准备
- ITS数据准备
COPD <- read_excel("COPD.xlsx")
head(COPD)
# A tibble: 6 x 5
month DDDs time level trend
<chr> <dbl> <dbl> <dbl> <dbl>
1 19-01 0.0929 1 0 0
2 19-02 0.0889 2 0 0
3 19-03 0.0943 3 0 0
4 19-04 0.0826 4 0 0
5 19-05 0.0888 5 0 0
6 19-06 0.0905 6 0 0
2 模型初拟合
- 分段线性拟合
mod1 <- lm(ACNY ~ time+level+trend, data = COPD)
summary(mod1)
3 自相关性检验
durbinWatsonTest(mod1) #durbinWatsontest
4 重新拟合模型
- 使用cochrane.orcutt解决自相关性
#cochrane.orcutt
mod2= cochrane.orcutt(mod1)
summary(mod2)
5 图形可视化
- 绘图
## 设置因子变量
COPD$predACNY<-fitted(mod2)
COPD$level=factor(COPD$level,labels = c( "政策前","政策后"))
## ggplot2
ggplot(data=COPD)+
geom_point(aes(x=month,y=ACNY,color=level),size=3,shape=19)+
scale_color_manual(values = c("#2A557A","#7FB7BB"))+
geom_vline(xintercept =25, linetype=4)+
geom_line(data = subset(COPD, time<= 24),aes(time,predACNY),size=1,color="#2A557A") +
geom_line(data = subset(COPD, time > 24),aes(time,predACNY),size=1,color="#7FB7BB") +
scale_x_discrete(breaks = c("15-01","15-06","15-12","16-06","17-01",
"17-06","17-12","18-06"))+
theme_classic()+
labs(x = "Time(months)",y = "ANCY")+
theme(legend.position = c(0.13,0.22),
legend.title=element_blank())
换个主题
ggplot(data=COPD)+
geom_point(aes(x=month,y=ACNY,color=level),size=3,shape=19)+
scale_color_manual(values = c("#2A557A","#7FB7BB"))+
geom_vline(xintercept =25, linetype=4)+
geom_line(data = subset(COPD, time<= 24),aes(time,predACNY),size=1,color="#2A557A") +
geom_line(data = subset(COPD, time > 24),aes(time,predACNY),size=1,color="#7FB7BB") +
scale_x_discrete(breaks = c("15-01","15-06","15-12","16-06","17-01",
"17-06","17-12","18-06"))+
theme_stata()+
labs(x = "Time(months)",y = "ANCY")+
theme(legend.position = c(0.13,0.22),
legend.title=element_blank())
