2
$\begingroup$

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:

enter image description here

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

$\endgroup$
3
  • 2
    $\begingroup$ I don't think they would. You are asking how to model a particular relationship. The fact that you wish to implement it in R is kind of irrelevant (most users over at Cross Validated know how to specify a model in R). I know quite a bit about R, and know that the model glm(outcome ~ time*iv*location, data = d, family = poisson) will give you different breaks and slopes for each location, but I don't know if that's an appropriate model. There are a few folks here who could answer, but more over at CV (plus the ones here who could answer also hang out at CV). $\endgroup$ Feb 19 at 19:24
  • 1
    $\begingroup$ Incidentally, it's a decent question, I'm just trying to get you the best quality answers for it. $\endgroup$ Feb 19 at 19:30
  • 2
    $\begingroup$ I have to agree with @AllanCameron here. The question is clearly about what the calculation should be - not about how to write code that implements that calculation. $\endgroup$ Feb 20 at 5:10

2 Answers 2

5
$\begingroup$

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:

enter image description here

$\endgroup$
0
$\begingroup$

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)

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.