Interrupted Time Series (Uncontrolled) with counterfactual and confidence intervals

library(tidyverse)
library(nlme)

# create a dummy dataframe
# Here the interruption takes place at week 51
df<-tibble(
  Time = 1:100,
  Intervention = c(rep(0,50),rep(1,50)),
  Post.intervention.time = c(rep(0,50),1:50),
  quantity.x = c(sort(sample(200:300,size = 50,replace = T),decreasing = T)+sample(-20:20,50,replace = T),c(sort(sample(20:170,size = 50,replace = T),decreasing = T)+sample(-40:40,50,replace = T)))
)


# Run a gls model of the interrupted time series.
model.a = model.a = gls(quantity.x ~ Time + Intervention + Post.intervention.time, data = df,correlation= corARMA(p=1, q=1, form = ~ Time),method="ML")

# Show a summary of the model
model.a.results<-summary(model.a)

# Make a new variable with predictions of the model
df<-df %>% mutate(
  model.a.predictions = predictSE.gls (model.a, df, se.fit=T)$fit,
  model.a.se = predictSE.gls (model.a, df, se.fit=T)$se
  )

# create counterfactual model
# here, the counterfactual is an extrapolation of the model up to week 50, created by filtering df and then dropping variables from the model
# if they are invariant before the intervention (i.e. Intervention == 0 and Post.intervention.time == 0 prior to week 51)
# So the original model a is
# quantity.x ~ Time + Intervention + Post.intervention.time
# but the counterfactual is
# quantity.x ~ Time

df2<-filter(df,Time<51)
model.b = gls(quantity.x ~ Time, data = df2, correlation= corARMA(p=1, q=1, form = ~ Time),method="ML")

# Make a new variable with predictions of the counterfactual model, providing the full 100 week data frame as 'newdata'
df<-df %>% mutate(
  model.b.predictions = predictSE.gls (model.b, df, se.fit=T)$fit,
  model.b.se = predictSE.gls (model.b, df, se.fit=T)$se
)

# plot the model with a counterfactual line

ggplot(df,aes(Time,quantity.x))+
  geom_ribbon(aes(ymin = model.b.predictions - (1.96*model.b.se), ymax = model.b.predictions + (1.96*model.b.se)), fill = "pink")+
  geom_line(aes(Time,model.b.predictions),color="red",lty=2)+
  geom_ribbon(aes(ymin = model.a.predictions - (1.96*model.a.se), ymax = model.a.predictions + (1.96*model.a.se)), fill = "lightgreen")+
  geom_line(aes(Time,model.a.predictions),color="black",lty=1)+
  geom_point(alpha=0.3)+
  ylim(0,300)