Interrupted time series: multiple interventions with a control group

I’ve seen examples of Interrupted Time Series (ITS) analysis with multiple interventions. And I’ve seen examples with control groups. But I haven’t found an example with both. In particular, I’ve really enjoyed this EdX course.

I’m currently working on a project that needs ITS with both multiple interventions and a control group. So, I’m making my own example. For my data, I have 2 interventions. Below is an example of how I’m setting up the data for regression in R. Hopefully I’ll have a full example including regression and analysis soon!

###########################################
# 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")