自然實驗（natural experiment）
隨著相關(guān)研究理論的不斷深入，其分析方法也在不斷拓展，目前常用的有雙重差分模型（difference-in-differences，DID）、工具變量法（instrumental variables，IV）、間斷時間序列模型（interrupted time series，ITS）、斷點回歸設(shè)計（regression discontinuity design，RDD）及其衍生的
方法等）
間斷時間序列分析（ITS）
ITS應(yīng)用
- 評估政策
- 時間范圍
- 因變量
- 模型
- 結(jié)果
R語言復(fù)現(xiàn)
1程序包加載
library(readxl)##讀入excel數(shù)據(jù)
library(car)##DWtest
library(orcutt)##Cochrane-Orcutt估計
library(prais)## prais_winsten估計
library(sandwich)##NeweyWest
library(lmtest)#線性回歸測試
library(ggplot2)##繪圖
1 數(shù)據(jù)準備
- ITS數(shù)據(jù)準備
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 自相關(guān)性檢驗
durbinWatsonTest(mod1) #durbinWatsontest
4 重新擬合模型
- 使用cochrane.orcutt解決自相關(guān)性
#cochrane.orcutt
mod2= cochrane.orcutt(mod1)
summary(mod2)
5 圖形可視化
- 繪圖
## 設(shè)置因子變量
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())
