Initial ITS models.Rmd

---
title: "Initial ITS models - draft"
output:
  # html_document:
  #   df_print: paged
  #   theme: readable
  word_document:
    reference_docx: report_template.docx
    fig_width: 8
    fig_height: 6
  # github_document: default
    
---

```{r setup, include=FALSE}
library(tidyverse)
library(readxl)
library(broom)
library(nlme)
library(car)

# options(scipen=10, digits=3)

testAutocorr <- function(model, data=NULL, max.lag = 10, time.points = 25) {
  data <- eval(model$call$data)  # Only works if 'lm' call has dataframe named in bracket
  print("Durbin Watson Test")
  print(dwt(model, max.lag = max.lag, alternative = "two.sided"))
  par(cex = 0.7, mai = c(0.1, 0.1, 0.2, 0.1))
  par(fig = c(0.03, 1, 0.8, 1))
  plot(
    data$Time[1:time.points],
    residuals(model)[1:time.points],
    type = 'o',
    pch = 16,
    col = "red",
    main = "Plotting residuals by time, with Autocorrelation and Partial-autocorrelation plots"
  )
  
  par(fig = c(0.03, 0.5, 0.05, 0.75), new = TRUE)
  acf(residuals(model))
  
  par(fig = c(0.55, 1, 0.05, 0.75), new = TRUE)
  acf(residuals(model), type = 'partial')
}

constructCIRibbon <- function(newdata, model) {
  newdata <- newdata %>%
    mutate(Predict = predict(model, newdata = newdata))
  mm <- model.matrix(as.formula(paste0("~ ", model$call$model[3])),
                     data = newdata)
  vars <- mm %*% vcov(model) %*% t(mm)
  sds <- sqrt(diag(vars))
  newdata <- newdata %>% mutate(lowCI = Predict - 1.96 * sds,
                                HiCI = Predict + 1.96 * sds)
}

printCoefficients <- function(model){
as_tibble(trimws(format(round(summary(model)$tTable, 3), nsmall=3))) %>%
  mutate(Coefficient = rownames(summary(model)$tTable)) %>% 
  select(Coefficient, Value, Std.Error, 'p-value') %>% 
  print()
}
```
# Importing and setting up data

Conceptions to under-18s were extracted from the dataset and coded into three datasets:

```{r include=FALSE}
all.UK.rates <-
  read_xlsx("Conception rates by age and country.xlsx", sheet = "Under 18")
```
### EngData - England only
```{r echo=FALSE}

EngData <- all.UK.rates %>% filter(Country == "England") %>%
  select(Country, '1992':'2016') %>%
  gather("Year", "Value",-1) %>%
  mutate(
    Cat1 = ifelse(Year < 1999, 0, 1),
    Cat2 = ifelse(Year <= 2007, 0, 1),
    Year = as.numeric(Year),
    Time = 1:nrow(.),
    Trend1 = ifelse(Cat1 == 0, 0, Year - 1998),
    Trend2 = ifelse(Cat2 == 0, 0, Year - 2007)
  )

head(EngData)
```

### EngScotContro data setup - Eng v Scotland as control
```{r echo=FALSE}
EngScotContro <-
  all.UK.rates %>% filter(Country == "Scotland" |
                            Country == "England") %>%
  gather("Year", "Value",-1) %>%
  filter(Year > 1991,!is.na(Value)) %>%
  arrange(Country) %>%
  mutate(
    England = ifelse(Country == "England", 1, 0),
    Year = as.numeric(Year),
    Time = Year - min(Year) + 1,
    Cat1 = ifelse(Year < 1999, 0, 1),
    Cat2 = ifelse(Year <= 2007, 0, 1),
    Trend1 = ifelse(Cat1 == 0, 0, Year - 1998),
    Trend2 = ifelse(Cat2 == 0, 0, Year - 2007)
  ) %>%
  mutate_at(., colnames(.)[5:9], list(Eng = ~ .*England))

head(arrange(EngScotContro, Year), 10)
```

### EngWalContro data setup - Eng vs Wales as control
```{r echo=FALSE}
EngWalContro <-
  all.UK.rates %>% filter(Country == "Wales" |
                            Country == "England") %>%
  gather("Year", "Value",-1) %>%
  filter(Year > 1991,!is.na(Value)) %>%
  arrange(Country) %>%
  mutate(
    England = ifelse(Country == "England", 1, 0),
    Year = as.numeric(Year),
    Time = Year - min(Year) + 1,
    Cat1 = ifelse(Year < 1999, 0, 1),
    Cat2 = ifelse(Year <= 2007, 0, 1),
    Trend1 = ifelse(Cat1 == 0, 0, Year - 1998),
    Trend2 = ifelse(Cat2 == 0, 0, Year - 2007)
  ) %>%
  mutate_at(., colnames(.)[5:9], list(Eng = ~ .*England))

head(arrange(EngWalContro, Year), 10)
```

# Initial exploratory models of England-only data

I first constructed a model with a single intervention at 1999 to test England data:

```{r echo=FALSE}
EngData %>% ggplot(aes(x = Year, y = Value, group = Cat1)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```
```{r}
EngMod <- lm(Value ~ Time + Cat1 + Trend1, data = EngData)
summary(EngMod)$coefficients
```

A visual inspection of the model and a test for autocorrelation raise concerns about this model:

```{r echo=FALSE}
testAutocorr(EngMod, max.lag = 15)
```

# Adding second intervention to England model

A second intervention at 2007 makes a better-fit model. This corresponds with a reevaluation and adaptation of the model after 2005 or potentially a phase-in period:

``` {r echo=FALSE}
EngMod2 <-
  lm(Value ~ Time + Cat1 + Trend1 + Cat2 + Trend2, data = EngData)

## Test for autocorrelation

testAutocorr(EngMod2, max.lag = 15)  # Potential MA3?
```
```{r}
EngMod2_q3 <- gls(
  Value ~ Time + Cat1 + Trend1 + Cat2 + Trend2,
  data = EngData,
  correlation = corARMA(q = 3, form =  ~ Time),
  method = "ML"
)
printCoefficients(EngMod2_q3)
```
```{r echo=FALSE, message=FALSE}
EM2q3_co <- EngMod2_q3$coefficients

EngMod2_q3_cfac <- tibble(
  Time = 8:25,
  Cat = c(rep.int(0, 9), rep(1, 9)),
  Trend1 = c(1:18),
  Trend2 = c(rep(0, 9), c(1:9)),
  Value = (
    EM2q3_co[1] +
      EM2q3_co[2] * Time +
      EM2q3_co[3] * Cat +
      EM2q3_co[4] * Trend1 * Cat
  )
)

EngData %>%
  mutate(predict = predict(EngMod2_q3)) %>%
  ggplot(aes(Time, Value, group = interaction(Cat1, Cat2))) +
  geom_point() +
  geom_line(aes(y = predict), col = "red") +
  geom_line(
    data = EngMod2_q3_cfac,
    aes(x = Time, y = Value, group = Cat),
    linetype = "dashed",
    col = "#FC8D62",
    size = 1,
    inherit.aes = FALSE
  ) +
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  geom_vline(xintercept = 16.5,
             linetype = "dotted",
             col = "#000000CC")
```

# modScot99_07 - Comparing with Scotland for three stages, split at 1999 and 2007

```{r message=FALSE}
modScot99_07_null <- gls(
  Value ~ Time +
    England +
    Time_Eng +
    Cat1 +
    Trend1 +
    Cat1_Eng +
    Trend1_Eng +
    Cat2 +
    Trend2 +
    Cat2_Eng +
    Trend2_Eng,
  data = EngScotContro,
  correlation = NULL,
  method = "ML"
)

printCoefficients(modScot99_07_null)
```

```{r echo=FALSE, message=FALSE}
modScot99_07_null_cfac <- EngScotContro[8:25, ] %>%
  select(Time, England, Cat1, Cat2, Trend1, Trend2, Time_Eng) %>%
  mutate(
    Cat1_Eng = c(rep(0, 9), rep(1, 9)),
    Trend1_Eng = c(rep(0, 9), 10:18),
    Cat2_Eng = 0,
    Trend2_Eng = 0
  )

modScot99_07_null_cfac <-
  modScot99_07_null_cfac %>%
  left_join(constructCIRibbon(modScot99_07_null_cfac, modScot99_07_null))

#** Graphing final model -----

EngScotContro %>%
  left_join(constructCIRibbon((filter(., England==1, Year>1998)), modScot99_07_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modScot99_07_null)) %>%  # Add Predicts for non-England
  ggplot(aes(
    Time,
    Value,
    col = Country,
    fill = Country,
    group = interaction(Country, Cat1, Cat2)
  )) +
  # Show all data points
  geom_point(data=EngScotContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modScot99_07_null_cfac,
    aes(
      x = Time,
      y = Predict,
      group = Cat2,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modScot99_07_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      group = Cat2,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  geom_vline(xintercept = 16.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```

# modScotPI - Comparing Scotland pre-1999 and post-2007, proposing phase-in

```{r message=FALSE}
modScotPI_null <-
  gls(
    Value ~ Time + England + Time_Eng + Cat2 + Trend2 + Cat2_Eng + Trend2_Eng,
    data = EngScotContro %>% filter(Year < 1999 | Year > 2007),
    correlation = NULL,
    method = "ML"
  )

printCoefficients(modScotPI_null)
```
```{r echo=FALSE, message=FALSE}
modScotPI_null_cfac <- tibble(
  Time = c(17:25),
  England = 1,
  Cat2 = 1,
  Trend2 = c(1:9),
  Time_Eng = c(17:25),
  Cat2_Eng = 0,
  Trend2_Eng = 0
)  # Emlulating England with no added effect

modScotPI_null_cfac <-
  modScotPI_null_cfac %>%
  left_join(constructCIRibbon(modScotPI_null_cfac, modScotPI_null))

#** Graphing final model -----

EngScotContro %>%
  filter(Year < 1999 | Year > 2007) %>%
  left_join(constructCIRibbon((filter(., England==1, Year>2007)), modScotPI_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modScotPI_null)) %>%
  ggplot(aes(Time, Value, col = Country, group = interaction(Country, Cat2))) +
  # Show all data points
  geom_point(data=EngScotContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modScotPI_null_cfac,
    aes(
      x = Time,
      y = Predict,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modScotPI_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Phase-in period greyed out
  geom_rect(
    xmin = 7.5,
    xmax = 16.5,
    ymin = 0,
    ymax = 60,
    fill = "grey",
    alpha = 0.01,
    inherit.aes = FALSE
  ) +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```

# modScotPI_99_07 - Comparing Eng vs Scot with phase-in and two interventions

```{r message=FALSE}
modScotPI_99_07_null <- gls(
  Value ~ Time +
    England +
    Time_Eng +
    Cat1 +
    Trend1 +
    Cat1_Eng +
    Trend1_Eng +
    Cat2 +
    Trend2 +
    Cat2_Eng +
    Trend2_Eng,
  data = {EngScotContro %>% 
      filter(Year < 1999 | Year > 2000) %>% 
      mutate(Trend1=ifelse(Cat1==0,0,Trend1-2),
             Trend1_Eng=Trend1*England)},
  correlation = NULL,
  method = "ML"
)

printCoefficients(modScotPI_99_07_null)
```
```{r echo=FALSE, message=FALSE}
modScotPI_99_07_null_cfac <- EngScotContro[10:25, ] %>%
  select(Time, England, Cat1, Cat2, Trend2, Time_Eng) %>%
  mutate(
    Trend1 = c(1:16),
    Cat1_Eng = c(rep(0, 7), rep(1, 9)),
    Trend1_Eng = c(rep(0, 7), 8:16),
    Cat2_Eng = 0,
    Trend2_Eng = 0
  )

modScotPI_99_07_null_cfac <-
  modScotPI_99_07_null_cfac %>%
  left_join(constructCIRibbon(modScotPI_99_07_null_cfac, modScotPI_99_07_null))

#** Graphing final model -----

EngScotContro %>%
  filter(Year < 1999 | Year > 2000) %>%
  mutate(Trend1=ifelse(Cat1==0,0,Trend1-2),
         Trend1_Eng=Trend1*England) %>% 
  left_join(constructCIRibbon((filter(., England==1, Year>2000)), modScotPI_99_07_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modScotPI_99_07_null)) %>%  # Add Predicts for non-England
  ggplot(aes(
    Time,
    Value,
    col = Country,
    fill=Country,
    group = interaction(Country, Cat1, Cat2)
  )) +
  # Show all data points
  geom_point(data=EngScotContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modScotPI_99_07_null_cfac,
    aes(
      x = Time,
      y = Predict,
      group = Cat2,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modScotPI_99_07_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      group = Cat2,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  geom_vline(xintercept = 16.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Phase-in period greyed out
  geom_rect(
    xmin = 7.5,
    xmax = 9.5,
    ymin = 0,
    ymax = 60,
    fill = "grey",
    alpha = 0.01,
    inherit.aes = FALSE
  ) +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```

# modWal99_07 - Comparing with Wales for changes at 1999 and 2007 

```{r message=FALSE}
modWal99_07_null <- gls(
  Value ~ Time +
    England +
    Time_Eng +
    Cat1 +
    Trend1 +
    Cat1_Eng +
    Trend1_Eng +
    Cat2 +
    Trend2 +
    Cat2_Eng +
    Trend2_Eng,
  data = EngWalContro,
  correlation = NULL,
  method = "ML"
)

printCoefficients(modWal99_07_null)
```
```{r echo=FALSE, message=FALSE}
modWal99_07_null_cfac <- tibble(
  Time = c(8:25),
  England = 1,
  Time_Eng = c(8:25),
  Cat1 = 1,
  Trend1 = c(1:18),
  Cat2 = c(rep(0,9), rep(1,9)),
  Trend2 = c(rep(0,9), 1:9),
  Cat1_Eng = c(rep(0,9), rep(1,9)),
  Trend1_Eng = c(rep(0,9), 10:18),
  Cat2_Eng = 0,
  Trend2_Eng = 0
  )  # Remove _Eng interactions (retaining 1st intervention interactions for 2nd int)

modWal99_07_null_cfac <-
  modWal99_07_null_cfac %>%
  left_join(constructCIRibbon(modWal99_07_null_cfac, modWal99_07_null))

EngWalContro %>%
  left_join(constructCIRibbon((filter(., England==1, Year>1998)), modWal99_07_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modWal99_07_null)) %>%  # Add Predicts for non-England
  ggplot(aes(
    Time,
    Value,
    col = Country,
    fill = Country,
    group = interaction(Country, Cat1, Cat2)
  )) +
  # Show all data points
  geom_point(data=EngWalContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modWal99_07_null_cfac,
    aes(
      x = Time,
      y = Predict,
      group = Cat2,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modWal99_07_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      group = Cat2,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  geom_vline(xintercept = 16.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```

# modWalPI - Comparing pre-1999 and post 2007 (assuming phase-in)

```{r message=FALSE}
modWalPI_null <-
  gls(
    Value ~ Time + England + Time_Eng + Cat2 + Trend2 + Cat2_Eng + Trend2_Eng,
    data = EngWalContro %>% filter(Year < 1999 | Year > 2007),
    correlation = NULL,
    method = "ML"
  )

printCoefficients(modWalPI_null)
```
```{r echo=FALSE, message=FALSE}
modWalPI_null_cfac <- tibble(
  Time = c(17:25),
  England = 1,
  Cat2 = 1,
  Trend2 = c(1:9),
  Time_Eng = c(17:25),
  Cat2_Eng = 0,
  Trend2_Eng = 0
)  # Emlulating England with no added effect

modWalPI_null_cfac <-
  modWalPI_null_cfac %>%
  left_join(constructCIRibbon(modWalPI_null_cfac, modWalPI_null))

EngWalContro %>%
  filter(Year < 1999 | Year > 2007) %>%
  left_join(constructCIRibbon((filter(., England==1, Year>2007)), modWalPI_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modWalPI_null)) %>%
  ggplot(aes(Time, Value, col = Country, group = interaction(Country, Cat2))) +
  # Show all data points
  geom_point(data=EngWalContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modWalPI_null_cfac,
    aes(
      x = Time,
      y = Predict,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modWalPI_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Phase-in period greyed out
  geom_rect(
    xmin = 7.5,
    xmax = 16.5,
    ymin = 0,
    ymax = 60,
    fill = "grey",
    alpha = 0.01,
    inherit.aes = FALSE
  ) +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```

# modWalPI_99_07 - comparing Eng vs Wal with phase-in period and two interventions

```{r message=FALSE}
modWalPI_99_07_null <- gls(
  Value ~ Time +
    England +
    Time_Eng +
    Cat1 +
    Trend1 +
    Cat1_Eng +
    Trend1_Eng +
    Cat2 +
    Trend2 +
    Cat2_Eng +
    Trend2_Eng,
  data = {EngWalContro %>% 
      filter(Year < 1999 | Year > 2000) %>% 
      mutate(Trend1=ifelse(Cat1==0,0,Trend1-2),
             Trend1_Eng=Trend1*England)},
  correlation = NULL,
  method = "ML"
)

printCoefficients(modWalPI_99_07_null)
```
```{r echo=FALSE, message=FALSE}
modWalPI_99_07_null_cfac <- EngWalContro[10:25, ] %>%
  select(Time, England, Cat1, Cat2, Trend2, Time_Eng) %>%
  mutate(
    Trend1 = c(1:16),
    Cat1_Eng = c(rep(0, 7), rep(1, 9)),
    Trend1_Eng = c(rep(0, 7), 8:16),
    Cat2_Eng = 0,
    Trend2_Eng = 0
  )

modWalPI_99_07_null_cfac <-
  modWalPI_99_07_null_cfac %>%
  left_join(constructCIRibbon(modWalPI_99_07_null_cfac, modWalPI_99_07_null))

#** Graphing final model -----

EngWalContro %>%
  filter(Year < 1999 | Year > 2000) %>%
  mutate(Trend1=ifelse(Cat1==0,0,Trend1-2),
         Trend1_Eng=Trend1*England) %>% 
  left_join(constructCIRibbon((filter(., England==1, Year>2000)), modWalPI_99_07_null)) %>%  # England CI ribbon
  mutate(Predict = predict(modWalPI_99_07_null)) %>%  # Add Predicts for non-England
  ggplot(aes(
    Time,
    Value,
    col = Country,
    fill=Country,
    group = interaction(Country, Cat1, Cat2)
  )) +
  # Show all data points
  geom_point(data=EngWalContro, show.legend = FALSE) +
  # Counterfactual trend lines
  geom_line(
    data = modWalPI_99_07_null_cfac,
    aes(
      x = Time,
      y = Predict,
      group = Cat2,
      col = "Control",
      fill = NULL
    ),
    linetype = "longdash",
    size = 1,
    inherit.aes = FALSE
  ) +
  # Counterfactual confidence intervals (not shown in legend)
  geom_ribbon(
    data=modWalPI_99_07_null_cfac,
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      group = Cat2,
      col=NULL,
      fill="Control"
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE,
    inherit.aes = FALSE) +
  # Model trend lines
  geom_line(aes(y=Predict), size = 1) +
  # Confidence intervals (not shown in legend)
  geom_ribbon(
    aes(
      x=Time,
      ymin = lowCI,
      ymax=HiCI,
      col=NULL,
      fill=Country
    ),
    alpha=0.5,
    size = 1,
    show.legend = FALSE) +
  # Intervention time points
  geom_vline(xintercept = 7.5,
             linetype = "dotted",
             col = "#000000CC") +
  geom_vline(xintercept = 16.5,
             linetype = "dotted",
             col = "#000000CC") +
  # Phase-in period greyed out
  geom_rect(
    xmin = 7.5,
    xmax = 9.5,
    ymin = 0,
    ymax = 60,
    fill = "grey",
    alpha = 0.01,
    inherit.aes = FALSE
  ) +
  # Display parameters
  scale_x_continuous(breaks = c(4, 9, 14, 19, 24),
                     labels = seq(1995, 2015, by = 5)) +
  theme(panel.background = element_blank(),
        legend.key  = element_blank()) +
  ylab("Rate of pregnancies to under-18s, per 1,000") +
  xlab("Year") +
  coord_cartesian(ylim = c(0, 60)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_colour_manual(
    breaks = c("England", "Wales", "Scotland", "Control"),
    values = c("Wales" = "#00AB39",
               "Scotland" = "#0072C6",
               "England" = "#CF142B",
               "Control" = "#F7D917"),
    aesthetics = c("colour", "fill"))
```