# Suitability of given AR1 model (using gls() in R) for nested time-series data?

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


Thank you

• How many times was the response or outcome variable Y measured pre- and post-intervention? – user158565 Apr 29 '17 at 6:48
• About 24 outcome observations/values are available per group for the pre- and post-intervention period (monthly observations made over ~2 years for each group-period). Hence probably no point looking beyond an AR1 model. – JupiterM104 Apr 29 '17 at 8:37
• Sorry I did not state the question clearly. I want to know is # of measurements on each subject (patient) pre- and post-intervention period. – user158565 Apr 29 '17 at 17:08
• Sorry. Each subject is technically a health facility, but its the same principle obviously: there were ~24 repeated observations made on each of the three health facilities in both the pre- and post-intervention periods. Edited question to use specific terminology, may help a little... – JupiterM104 Apr 30 '17 at 8:32

correlation = corAR1(form = ~ x | grp/period),

It seems grp/period separates the data into 6 groups, and assumes the independent between groups. Especially it means the independent between the observations at pre- and post-intervention from the SAME facility. Maybe this assumption is not true. So I prefer

correlation = corAR1(form = ~ x | grp),

I think others are correct.

• Thank you for your input. This certainly makes sense, and I agree it's more sensible to assume shared error structure within groups across time periods. However, are you able to provide any published exampled/textbook reference supporting this (at least in the general sense), as that was a condition of the bounty? – JupiterM104 May 4 '17 at 7:59
• If you are hungry, you need to eat something; if you are thirsty, you need to drink something. In your case, it is obvious that measurements from the same subject pre- and post-intervention are correlated (especially for control), so we need to incorporate this correlation into the model. So that is my suggestion to use grp only. For reference, I think it is based on common sense of the statistics and you do not need reference, If you insist on reference, any good book talking about mixed model can be added as reference. – user158565 May 4 '17 at 19:36
• Fair enough. Accepted but I have been totally distracted by work the last few days and the bounty has ended, and after a bit of reading I see that's it's lost and I can't assign it to your answer (or even get it back), so I'm sorry for that, but thank you for the help. – JupiterM104 May 8 '17 at 7:29

Can I ask why you are using the analysis strategy rather than another method, for example, why not use a spline model (or piece-wise, or split, or whatever your field calls a continuous time series with more than one slope)? Coding a continuous time variable along with your fixed levels (period one, period 2) seems like it would provide a little more clarity.

You have time period 1 (before treatment) and time period 2 (after treatment), you want to compare the slope and rate of change between the two so why not do the following model-

    model = y ~ X + T1 + T2 + I(T1)^2 + I(T2)^2 +
(X : T1) + (X : T2) + ( X:(I(T1)^2)) + (X:(I(T2)^2) +
(N [whatever random slope]|group), data=data)


The T1 and T2 terms will give you your slope, and the I(T1)^2 term will give you the rate of change.

Also, for reference, this coding using the lme4 package.

• Not sure I fully follow your explanation here. To me my model is a continuous time series with more than one slope, i.e. a slope for the pre-intervention period and the post-intervention period for each of the groups. Also I don't see how your suggested model accounts for serial autocorrelation in the errors, and it looks like it's treating the groups as random effects, when I want to treat them as fixed effects, because they are treatment groups, not simply repeated observation groups... – JupiterM104 May 4 '17 at 8:04