使用 lme4 分析时间序列数据时，如何同时包含阶跃变化和斜率变化？

Question

我有一些时间序列数据，包括四个不同的位置。在某个时间点有干预（每个地点不同）。

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))

我正在尝试拟合一个模型，该模型允许在每个位置进行不同的步长变化和坡度变化。目前我有一个允许改变步长的模型：

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)

特别是，您可以看到位置 2 的斜坡建模不佳：

您将如何修改

glmer

公式以使其包含斜率变化？

Answer 1

您的模型当前有一个固定的

time

参数，该参数对于所有观察结果都是相同的。该估计值非常轻微为负，因此您的所有预测都呈轻微下降趋势。请注意，干预时的跳跃主要是通过

lines

调用连接的这些点的产物；它显示了不同

iv

值的观测值之间的距离以及它们的随机截距之间的差异，但不是来自模型本身的“斜率”。

为了允许这些斜率不同，您必须通过包含额外的模型参数来放宽所有观测值都遵循相同

time

效应的假设。只有您可以判断什么是最有意义的，但一个明显的方法是添加与

iv

的交互，以便干预前后的斜率可以不同：

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

问题仍然是，所有

location

在干预前仍然具有相同的斜率，其中位置4具有最多的数据，但向上的斜率也没有那么陡。您可以为每个

time

添加随机

location

斜率，尽管您可能会开始遇到拟合如此广泛的随机效应设计的麻烦 - 至少

lme4

无法获得

glmmTMB

可以收敛的模型做作业：

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

正如评论中所建议的，在完整的三向交互中添加

location

作为固定效果也允许所有观察集群之间的斜率不同，尽管这不再是混合模型：

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

无论哪种方式，挑战都在于模型中现在有相当多的参数，因此解释变得更加复杂。我还纯粹关注如何添加此类参数，实际上检查模型拟合情况是一个完全不同的主题。这些模型的图如下所示：