I have nested time series data where the outcome (number of visits per month) is available for 24 months (repeated measures) across two time periods (pre- and post-intervention), for three health facilities (a control, intervention 1 and intervention 2, where the interventions differ). I wish to see if once the interventions were applied whether the number of visits changed in the two intervention health facilities vs the control, both in terms of their level and trend.
I am happy setting up my fixed-effects regression model (three way interaction between 1) time (month), 2) pre- and post-intervention dummy indicator and 3) health facility indicator (1-3)), but I am not sure I am correctly dealing with the serial autocorrelation in the outcome, which appears to follow an AR1 process within each time period for each health facility.
I believe I should be able to deal with this using a generalized least squares model with a suitable AR1 error structure for each group applied within each time period, which can be done in R using the nlme package and gls function. My attempted code is below. Feedback/confirmation of suitability very welcome.
# Create lag-1 correlated outcome data and covariates n <- 100 # Number of observations within each group and time-period gn <- 3 # Number of groups tn <- 2 # Number of time periods grp <- rep(c("c", "i1", "i2"), each = n*tn) # Treatment group dummy covariate period <- rep(rep(c("pre", "pst"), each = n), gn) # Pre-post intervention period dummy covariate xt <- rep(1:n, times = gn*tn) # Time covariate # Create a serially correlated error term for each time period within each group, # and use as the outcome as no need to add any any fixed effects for this example y <- c(filter(rnorm(n), filter=rep(1,2), circular=TRUE), filter(rnorm(n), filter=rep(1,2), circular=TRUE), filter(rnorm(n), filter=rep(1,2), circular=TRUE), filter(rnorm(n), filter=rep(1,2), circular=TRUE), filter(rnorm(n), filter=rep(1,2), circular=TRUE), filter(rnorm(n), filter=rep(1,2), circular=TRUE)) d <- data.frame(y=y, grp=grp, period=period, xt=xt) # Create AR1 model with fixed effects of group-intervention period interaction, # and AR1 correlated errors for each time period within each group m1 <- gls(y ~ x*grp*period, correlation = corAR1(form = ~ x | grp/period), data = d) acf(resid(m1, type = "n")) # Check lagged autocorrelation