When analysing time series data with lme4, how do you include both a step-change and a slope-change?

Question

I have some time series data including four different locations. There is an intervention at a certain point in time (different in each location).

d <- structure(list(location = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), 
    time = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
    13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
    25L, 26L, 27L, 28L, 29L, 30L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
    8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 
    20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 1L, 
    2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
    15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 
    27L, 28L, 29L, 30L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
    11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 
    23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L), iv = c(0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L), outcome = c(27L, 
    23L, 13L, 31L, 27L, 31L, 27L, 24L, 31L, 20L, 13L, 14L, 12L, 
    14L, 12L, 16L, 18L, 14L, 20L, 15L, 9L, 15L, 13L, 24L, 13L, 
    12L, 14L, 17L, 12L, 9L, 22L, 21L, 22L, 30L, 28L, 28L, 32L, 
    30L, 51L, 42L, 43L, 32L, 45L, 39L, 43L, 26L, 22L, 25L, 14L, 
    21L, 22L, 17L, 8L, 12L, 14L, 13L, 14L, 11L, 7L, 6L, 20L, 
    24L, 22L, 27L, 27L, 22L, 23L, 27L, 27L, 24L, 26L, 35L, 32L, 
    26L, 22L, 29L, 26L, 38L, 24L, 15L, 13L, 15L, 9L, 12L, 9L, 
    4L, 8L, 7L, 8L, 4L, 37L, 22L, 27L, 24L, 33L, 20L, 28L, 26L, 
    23L, 21L, 29L, 28L, 26L, 24L, 31L, 27L, 24L, 24L, 18L, 18L, 
    24L, 19L, 24L, 27L, 30L, 13L, 23L, 15L, 13L, 16L)), class = "data.frame", row.names = c(NA, 
-120L))

I am trying to fit a model that allows both a different step change and slope change in each location. Currently I have a model that allows the step change:

m <- glmer(outcome ~ time  + (1 + iv|location), family = 'poisson', data = d)
d$p <- predict(m, newdata = d, type = 'response')

par(mfrow = c(4, 1), mar = c(1, 0, 0, 0), oma = c(5, 5, 0, 0))
for(i in 1:4) {
  plot(1, type = 'n', xlim = c(0, 30), ylim = c(0, 50), axes = F, xlab = NA, ylab = NA)
  with(d[d$location == i,], {
    rect(0, 0, 31, 50)
    points(time, outcome)
    lines(time, p)
    if(i == 4) axis(1, pos = 0)
    axis(2, las = 2, pos = 0)
    text(28, 50 * 0.9, paste0('Location ', i), font = 2)
  })
}
mtext('Time period', side = 1, outer = T, line = 2.7, cex = 0.8)
mtext('Event count', side = 2, outer = T, line = 3, cex = 0.8)

In particular, you can see the slopes for location 2 are poorly modelled:

How would you modify the glmer formula so it includes a change-in-slope?

Answer 1

Your model currently has a single fixed time parameter, which is the same for all observations. The estimate is very slightly negative, so all your predictions trend very slightly downwards. Note that the jump at the intervention is mostly an artifact of these points being connected by the lines call; it shows the distance between the observations of different iv values and the difference between their random intercepts but isn't a 'slope' coming from the model itself.

To allow these slopes to differ you have to relax the assumption that all observations share the same time effect, by including additional model parameters. Only you can judge what makes most sense here, but an obvious one would be to add the interaction with iv so that the slopes pre- and post-intervention can differ:

m <- glmer(outcome ~ time*iv + (iv|location), family = poisson, 
           data = d)

The problem remains that all location still have the same slope pre-intervention, where Location 4 has the most data but also not as steep an upward slope. You could add a random time slope for each location, though you might start running into trouble fitting more extensive random effect designs -- at least lme4 wasn't able to get convergence for the model that glmmTMB could make work:

m <- glmmTMB::glmmTMB(outcome ~ iv*time + (iv+time|location), 
                      family = poisson, data=d)

As suggested in the comments adding location as a fixed effect in a full 3-way interaction also allows slopes to differ across all observation clusters, though this is no longer a mixed model:

m <- glm(outcome ~ time*iv*location, family = poisson, data=d)

The challenge either way becomes that there's quite a few more parameters in your model now, so interpretation becomes a bit more complicated. I've also focused purely on how you might go about adding such parameters, actually checking the model fits is a whole different topic. Plots for these models are shown below:

Answer 2

Having looked into this, I think a good way to approach it is to combine the intervention and location variables.

We want the effect of time to vary by by location and before/after the intervention. So we can first combine the intervention and location variables:

d$location_iv <- paste0(d$location, '-', d$iv)

Then specify a random intercept and a random slope for time across levels of location/intervention:

m <- glmer(outcome ~ (1 + time|location_iv), family = 'poisson', data = d)