用R语言实现“中断时间序列分析” --- R code for ITS

#### Rcode for ITS ####
####用R语言实现"中断时间序列分析"###2023.03.28###
####联系方式:shaofengsui@163.com  ###
####推荐一篇参考文献,请大家下载学习###
####Bernal, J.L., Cummins, S. and Gasparrini, A. (2020) Corrigendum to interrupted time series regression for the evaluation of public health interventions: a tutorial. Int J Epidemiol 49(4), 1414.
############################################################################

###R code 编码如下###

# Install packages required for the analysis (uncomment if needed)
##install.packages("lmtest") ; install.packages("Epi")
##install.packages("tsModel"); install.packages("vcd")
#### load the packages  安装包####
library(foreign) ; 
library(tsModel) ;
library(zoo) ;
library("lmtest") ; 
library("Epi") ;
library("splines") ; 
library("grid") ; 
library("vcd")
setwd("D:/R")  ####建立工作文件路径######
### read data from csv file  读取资料###
data <- read.csv("L2.csv")   ###模拟数据可发邮件向作者索取##
head(data)   ##显示数据表的前6行
View(data)   ##看一下数据表

# 变量的含义
# year   年
# sason  季节
# time = 时间的序号
# thms = 三卤甲烷浓度（µg/L）
# period = 干预措施，干预之前为 0 , 干预后为1###
###本研究的干预措施为藻类监测与应对方案的实施####

###描述性分析#######################################
# 检查资料的特点是第一步，线性模型是否合适,是否有季节性变动##观察干预前的趋势是否稳定
## 绘制散点图###
# 开始画图，不包括点和X轴###
plot(data$thms,type="n",ylim=c(0,35),xlab="Year", ylab="concentrations(µg/L)",
     bty="c",xaxt="n") 

##box type 绘图区域边框类型可选 "7", "c", "u", or "]" ##
###xaxt ,"x axis type". x轴标签类型,默认是xaxt="s"###
### 方形阴影，将干预后的时期涂成灰色###
rect(12,0,40,35,col=grey(0.9),border=F)
###等价于下面语句###
### rect(xleft = 9,ybottom = 0,xright = 24,ytop = 12,col=grey(0.9),border=T)
# 绘制干预前期观察到的率###
points(data$thms[data$period==0],cex=0.7)   ###cex 控制文字或点的大小###
##points(data$chcl3,cex=0.7) 

###specify the x-axis (i.e. time units) 设置x 轴的时间单位##
axis(1,at=0:9*4,labels=F)    ###9年，每年4个季度####
axis(1,at=0:8*4+2,tick=F,labels=2011:2019)   
###at在哪里添加刻度，在6个月的地方标注年号，可将2改为3试试效果###
# 加标题
title("L water plant 2011-2019")

# 对数据表中的变量作描述性统计分析###
summary(data)

#水源地污染物污染物治理措施实施前后浓度描述#
summary(data$thms)  ###总体###
summary(data$thms[data$period==0])  ###干预前###
summary(data$thms[data$period==1])  ###干预后###

#######################################
线性模型(lm)的回归分析
#######################################
###本研究的模型：Yt=β0 + β1*T + β2*Xt + β3 (T-T0 )⋅Xt ##
### Yt为thms, T为time, Xt 为干预阶段，T0=12 ###
###level模型，不考虑斜率改变##
model1 <- lm(thms ~  time + period, data=data) 

### level + slope 模型，考虑斜率的改变###
model2 <- lm(thms ~  time + period + I(time-12):period, data=data) 
###I(time-12):period  ###  I:交互作用的英文缩写，T0=12##

AIC(model1,model2) 
###选择最小的AIC(Akaike information criterion)值,来选择进入方程的变量#
###可利用AIC准则来判断模型拟合的优劣###

###以下工作基于model2,考虑斜率的改变###
model2 <- lm(thms ~ time + period + I(time-12):period, data=data) 
##回归方程常见的参数###
summary(model2)

###***统计结果的解释***###结果因模拟数据的改变，可能会有所不同###
###intercept:干预前的回归线在Y轴的截距，本例7.7736###
###time：干预前的回归线的斜率β1=1.5614，每过一个单位时间（季度）Y值的变化量##
###period:干预措施引起的污染物水平leve的变化β2=-6.1799 ####
###I(time-12)：period,交互作用项，干预引起的回归线斜率slope的改变量β3=干预后回归线的斜率-β1=-2.1872 ##
###干预后回归线的斜率 = β1 + β3 = 1.5614 +（-2.1872）= -0.6258###

summary(model2)$dispersion
confint(model2)   ###回归方程中各自变量回归系数的可信区间

# 用 0.1 time units 产生一个新的数据框，提高图片的精度##
datanew <- data.frame(thms=mean(data$thms),period=rep(c(0,1),c(120,240)),
                      time= 1:360/10,season=rep(1:40/10,9))  ###

# 根据模型生成绘图需要的预测值#
pred1 <- predict(model2,type="response",datanew) 

#绘图###
plot(data$thms,type="n",ylim=c(00,40),xlab="年份",ylab="THMs浓度(µg/L)",
     bty="l",xaxt="n")

rect(12,00,38,42,col=grey(0.9),border=F)
points(data$thms,cex=0.7)
axis(1,at=0:9*4,labels=F)
axis(1,at=0:8*4+2,tick=F,labels=2011:2019)
lines((1:360/10),pred1,col=4)   ###加回归线###col=4，蓝色线条##
###title("Influence of Algae control on THMs levels, 2011-2019")##

###回归线画出来了###