Continued: ITS: multiple interventions and control

Continuing my example of an interrupted time series with multiple interventions and a control. Here is the completed example in R — including examining the model residuals for evidence of autoregressive (AR) and moving average (MA) processes.

###########################################
# Modifid by Ryan L. Melvin From...
# ITSx: Week 4 Multiple Interventions
# Michael Law (michael.law@ubc.ca)
# October 2015
###########################################

###################################
# Load the necessary libraries
###################################

library(readxl)


########################
# Read in the dataset
########################

# Data set needs to have a column with time labels (for me month),
# value of interest (for me case delay minutes),
# and a colum indicating control or treatment group. (for me, 1 for treatment, 0 for control)
data <- read_excel("Data.xlsx")

# Setup variables
# I have 116 months of data
data$time <- c(1:116,1:116)
# the first intervention happens after 61 months
# I have to set up the level and trend twice -- once for control and once for treatment

# the first intervention for my data set is the announcement of a policy
data$announce_level <- c(rep(0,61),rep(1,55),rep(0,61),rep(1,55))
data$announce_trend <- c(rep(0,61),1:55,rep(0,61),1:55 )

# the second intervention is the implementation of the policy
# This happens after the 72nd month in the data
data$implemented_level <- c(rep(0,72),rep(1,44),rep(0,72),rep(1,44))
data$implemented_trend <- c(rep(0,72),1:44,rep(0,72),1:44 )

data$treatmenttime <- data$treatment * data$time

data$treatment_announce_level <- data$treatment * data$announce_level
data$treatment_announce_trend<- data$treatment * data$announce_trend

data$treatment_implemented_level <- data$treatment * data$implemented_level
data$treatment_implemented_trend <- data$treatment * data$implemented_trend 

########################
# Initial Plot
########################

# Plot the time series for the treatment delays
plot(data$time[1:116],data$AvgDelay[1:116],
     ylab="Average Case Delay (minutes)",
     ylim=c(0,2000),
     xlab="Month",
     type="l",
     col="red",
     xaxt="n")

# Add in control group non-treatment delays
points(data$time[117:232],data$AvgDelay[117:232],
       type='l',
       col="blue")

# Add x-axis year labels
axis(1, at=1:116, labels=data$Month[1:116])

# Add in the points for the figure
points(data$time[1:116],data$AvgDelay[1:116],
       col="red",
       pch=20)

points(data$time[117:232],data$AvgDelay[117:232],
       col="blue",
       pch=20)

# Label the first  change
abline(v=61.5,lty=2)

# Label the second  change
abline(v=72.5,lty=2)

# Add in a legend
legend('topright', legend=c("treatment","No treatment"),
       col=c("red","blue"),pch=20)

title(main="Average Case Delays")


# Initial OLS Model
# A preliminary OLS regression
model_ols <- lm(AvgDelay ~ time + treatment + treatmenttime
                + announce_level + announce_trend + treatment_announce_level + treatment_announce_trend
                + implemented_level + implemented_trend + treatment_implemented_level + treatment_implemented_trend
                , data=data)
summary(model_ols)

# check autocorrelation

# Durbin-watson test, 12 time periods
dwt(model_ols,max.lag=12,alternative="two.sided")

# Graph the residuals from the OLS regression to check for serially correlated errors
plot(residuals(model_ols),
     type='o',
     pch=16,
     xlab='Time',
     ylab='OLS Residuals',
     col="red")
abline(h=0,lty=2)

# Plot ACF and PACF
# Set plotting to two records on one page
par(mfrow=c(1,2))

# Produce plots
acf(residuals(model_ols))
acf(residuals(model_ols),type='partial')

# Fit the GLS regression model using the AR and MA components determined
model_p4 <- gls(AvgDelay ~ time + CIED + ciedtime
                + announce_level + announce_trend + cied_announce_level + cied_announce_trend
                + implemented_level + implemented_trend + cied_implemented_level + cied_implemented_trend
                , data=data
                ,correlation=corARMA(p=4
                                     #,q=1
                                     ,form=~time|treatment)
                ,method="ML")
summary(model_p4)


########################
# Plot results
#########################

# Produce the plot, first plotting the raw data points
plot(data$time,data$strips_pt,
     ylab="Test Strips per 1,000 people",
     ylim=c(0,300),
     xlab="Month",
     pch=20,
     col="pink",
     xaxt="n")

# Add x axis with dates
axis(1, at=1:72, labels=data$yearmonth)

# Add line indicating the policy changes
abline(v=30.5,lty=2)
abline(v=56.5,lty=2)

# Plot the first line segment
lines(data$time[1:30], fitted(model_p4)[1:30], col="red",lwd=2)

# Plot the second line segment
lines(data$time[31:56], fitted(model_p4)[31:56], col="red",lwd=2)

# Plot the third line segment
lines(data$time[57:72], fitted(model_p4)[57:72], col="red",lwd=2)

# And the first counterfactual
segments(31, model_p4$coef[1]+model_p4$coef[2]*31,
         56, model_p4$coef[1]+model_p4$coef[2]*56,
         lty=2, lwd=2, col='red')

# And the second counterfactual
segments(57, model_p4$coef[1] + model_p4$coef[2]*57 +
                 model_p4$coef[3] + model_p4$coef[4]*27,
         72, model_p4$coef[1] + model_p4$coef[2]*72 +
                 model_p4$coef[3] + model_p4$coef[4]*42,
         lty=2, lwd=2, col='red')

# END