Importing and setting up data

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

EngData - England only

EngScotContro data setup - Eng v Scotland as control

EngWalContro data setup - Eng vs Wales as control

Initial exploratory models of England-only data

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

EngMod <- lm(Value ~ Time + Cat1 + Trend1, data = EngData)
summary(EngMod)$coefficients
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 40.800000  2.4243796 16.829047 1.147315e-13
## Time         0.750000  0.5421078  1.383489 1.810468e-01
## Cat1         4.016667  2.4104749  1.666338 1.104935e-01
## Trend1      -2.295614  0.5575524 -4.117306 4.908456e-04

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

## [1] "Durbin Watson Test"
##  lag Autocorrelation D-W Statistic p-value
##    1      0.77051355     0.3725459   0.000
##    2      0.50992122     0.8417204   0.002
##    3      0.17667011     1.4588131   0.236
##    4     -0.02768092     1.8009124   0.970
##    5     -0.25166863     2.2352735   0.208
##    6     -0.42708533     2.5828226   0.018
##    7     -0.56415048     2.8136162   0.000
##    8     -0.56958102     2.6501264   0.000
##    9     -0.43563602     2.1664036   0.018
##   10     -0.25030550     1.5869374   0.256
##   11     -0.09536613     1.2216464   0.760
##   12      0.01388605     0.9759968   0.812
##   13      0.07924590     0.8377393   0.768
##   14      0.12616731     0.7167512   0.682
##   15      0.14731591     0.6190418   0.612
##  Alternative hypothesis: rho[lag] != 0

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:

## [1] "Durbin Watson Test"
##  lag Autocorrelation D-W Statistic p-value
##    1      0.18335719     1.4090132   0.004
##    2     -0.21037978     2.1848548   0.954
##    3     -0.40649345     2.4644492   0.244
##    4     -0.18427433     1.6585720   0.766
##    5     -0.18777350     1.5415165   0.874
##    6      0.05529295     1.0427883   0.250
##    7      0.13015466     0.8708738   0.162
##    8      0.07780085     0.9297698   0.312
##    9      0.07435073     0.9254829   0.326
##   10      0.05638279     0.9025666   0.380
##   11      0.14285767     0.7206965   0.196
##   12     -0.02674649     1.0574924   0.954
##   13     -0.11503662     1.2260738   0.396
##   14     -0.21481723     1.4232226   0.046
##   15     -0.14384408     1.2723560   0.086
##  Alternative hypothesis: rho[lag] != 0

EngMod2_q3 <- gls(
  Value ~ Time + Cat1 + Trend1 + Cat2 + Trend2,
  data = EngData,
  correlation = corARMA(q = 3, form =  ~ Time),
  method = "ML"
)
printCoefficients(EngMod2_q3)
## # A tibble: 6 x 4
##   Coefficient Value  Std.Error `p-value`
##   <chr>       <chr>  <chr>     <chr>    
## 1 (Intercept) 40.931 0.749     0.000    
## 2 Time        0.693  0.178     0.001    
## 3 Cat1        -1.143 0.944     0.241    
## 4 Trend1      -1.140 0.191     0.000    
## 5 Cat2        1.230  0.869     0.173    
## 6 Trend2      -2.246 0.136     0.000

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

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)
## # A tibble: 12 x 4
##    Coefficient Value  Std.Error `p-value`
##    <chr>       <chr>  <chr>     <chr>    
##  1 (Intercept) 38.330 1.786     0.000    
##  2 Time        0.950  0.344     0.009    
##  3 England     2.470  2.009     0.227    
##  4 Time_Eng    -0.200 0.400     0.621    
##  5 Cat1        -4.372 1.154     0.001    
##  6 Trend1      -0.885 0.371     0.023    
##  7 Cat1_Eng    2.799  1.583     0.085    
##  8 Trend1_Eng  -0.298 0.447     0.509    
##  9 Cat2        0.849  1.035     0.417    
## 10 Trend2      -2.805 0.198     0.000    
## 11 Cat2_Eng    0.503  1.463     0.733    
## 12 Trend2_Eng  0.544  0.281     0.061

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

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)
## # A tibble: 8 x 4
##   Coefficient Value   Std.Error `p-value`
##   <chr>       <chr>   <chr>     <chr>    
## 1 (Intercept) 38.330  1.912     0.000    
## 2 Time        0.950   0.368     0.017    
## 3 England     2.470   2.150     0.263    
## 4 Time_Eng    -0.200  0.429     0.645    
## 5 Cat2        -11.487 4.168     0.012    
## 6 Trend2      -3.690  0.398     0.000    
## 7 Cat2_Eng    0.618   5.025     0.903    
## 8 Trend2_Eng  0.245   0.479     0.613

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

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)
## # A tibble: 12 x 4
##    Coefficient Value  Std.Error `p-value`
##    <chr>       <chr>  <chr>     <chr>    
##  1 (Intercept) 38.330 1.623     0.000    
##  2 Time        0.950  0.312     0.005    
##  3 England     2.470  1.825     0.186    
##  4 Time_Eng    -0.200 0.364     0.586    
##  5 Cat1        -7.966 1.566     0.000    
##  6 Trend1      -0.504 0.364     0.176    
##  7 Cat1_Eng    3.387  2.040     0.107    
##  8 Trend1_Eng  -0.546 0.450     0.233    
##  9 Cat2        0.003  0.984     0.997    
## 10 Trend2      -3.187 0.226     0.000    
## 11 Cat2_Eng    1.056  1.392     0.454    
## 12 Trend2_Eng  0.792  0.320     0.019

modWal99_07 - Comparing with Wales for changes at 1999 and 2007

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)
## # A tibble: 12 x 4
##    Coefficient Value  Std.Error `p-value`
##    <chr>       <chr>  <chr>     <chr>    
##  1 (Intercept) 44.729 1.031     0.000    
##  2 Time        1.339  0.231     0.000    
##  3 England     -3.929 1.458     0.010    
##  4 Time_Eng    -0.589 0.326     0.079    
##  5 Cat1        -4.692 1.215     0.000    
##  6 Trend1      -1.973 0.279     0.000    
##  7 Cat1_Eng    3.120  1.719     0.077    
##  8 Trend1_Eng  0.789  0.395     0.053    
##  9 Cat2        1.553  1.161     0.189    
## 10 Trend2      -2.135 0.223     0.000    
## 11 Cat2_Eng    -0.200 1.642     0.904    
## 12 Trend2_Eng  -0.127 0.315     0.690

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

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)
## # A tibble: 8 x 4
##   Coefficient Value   Std.Error `p-value`
##   <chr>       <chr>   <chr>     <chr>    
## 1 (Intercept) 44.729  1.120     0.000    
## 2 Time        1.339   0.250     0.000    
## 3 England     -3.929  1.583     0.020    
## 4 Time_Eng    -0.589  0.354     0.109    
## 5 Cat2        -20.893 3.194     0.000    
## 6 Trend2      -4.108  0.303     0.000    
## 7 Cat2_Eng    10.024  4.517     0.036    
## 8 Trend2_Eng  0.663   0.429     0.135

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

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)
## # A tibble: 12 x 4
##    Coefficient Value   Std.Error `p-value`
##    <chr>       <chr>   <chr>     <chr>    
##  1 (Intercept) 44.729  0.972     0.000    
##  2 Time        1.339   0.217     0.000    
##  3 England     -3.929  1.375     0.007    
##  4 Time_Eng    -0.589  0.307     0.064    
##  5 Cat1        -10.525 1.522     0.000    
##  6 Trend1      -1.579  0.307     0.000    
##  7 Cat1_Eng    5.946   2.152     0.009    
##  8 Trend1_Eng  0.529   0.435     0.232    
##  9 Cat2        0.682   1.146     0.556    
## 10 Trend2      -2.529  0.263     0.000    
## 11 Cat2_Eng    0.377   1.620     0.817    
## 12 Trend2_Eng  0.134   0.372     0.721