自然实验（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())

2023-03-03