如何在中断的分段时间序列回归中向ggplot添加线性段

Question

我已经安装了一个中断的时间序列回归来计算数据,并希望显示与此类似的结果

取自:Lindstrand A,Bennet R,Galanis I,et al.引入肺炎球菌结合疫苗后鼻窦炎和肺炎住院治疗.儿科.2014; 134(6):e1528-36.DOI:10.1542/peds.2013-4177.

具体来说,我正在尝试(和失败)再现的是分别添加品红色和青色趋势线.我一直试图在ggplot中这样做.问题是我的模型是合适的,glm(family = poisson)因此系数不是原始尺度.更复杂的是,我提供了风险人口作为偏移,glm(count ~ ., offset(log(at_risk)), family = poisson, data = df)但是想要(count / at_risk)*1000在Y轴上显示数据.

set.seed(42)
int = 85
df <- data.frame(
    count = as.integer(rpois(132, 9) + rnorm(132, 1, 1)),
    time = 1:132,
    at_risk = rep(
        c(4305, 4251, 4478, 4535, 4758, 4843, 4893, 4673, 4522, 4454, 4351),
        each  = 12
    )
)
df$month <- factor(month.name, levels = month.name)
df$intv <- ifelse(df$time >= int, 1, 0)
df$intv_trend <- c(rep(0, (int - 1)),
                   1:(length(unique(df$time)) - (int - 1)))
df <-
    df %>%
    mutate(lag_count = dplyr::lag(count))

fit <- glm(
    count ~ month + time + intv + intv_trend +
        log(lag_count) + offset(log(at_risk)),
    family = "poisson",
    data = df
)
df$fit <- exp(c(NA, predict(fit)))


ggplot(df, aes(x = time, y = (fit / at_risk) * 1000)) +
    geom_line()

(我已经绘制了我希望能够创建到生成的ggplot行图中的行)

time伪方程给出了一个连续的长期趋势,count ~ intercept + B1 * time我想截断它,使其大致停止time = 72.这类似于上图中的洋红色线.干预intv发生在time = 85这导致在水平的变化intv斜率与变化intv_trend.intv效果线相对于时间的伪代码count ~ intercept + intv + B1 * time + B2* intv_trend类似于上面的青色线.

我试过geom_abline()不同版本的exp(coef(fit)[1] ...等等,但是我无法在剧情中看到这条线.

有什么想法吗？

Answer 1

正如我在评论中所说,如果你有办法识别变化点,你可以添加一个名为,比如group标记预测线的第一部分Control和第二部分Intervention(或者你喜欢的任何标签)的列.然后在你的情节中使用group作为颜色美学来获得两条不同的线条.在下面的代码中,我手动添加了分组变量.要获得有关数据规模的预测,请添加type="response"到predict.

首先,设置数据:

library(ggplot2)
library(dplyr)

int = 85
set.seed(42)
df <- data.frame(
  count = as.integer(rpois(132, 9) + rnorm(132, 1, 1)),
  time = 1:132,  
  at_risk = rep(
    c(4305, 4251, 4478, 4535, 4758, 4843, 4893, 4673, 4522, 4454, 4351),
    each  = 12
  )
)

df$month <- factor(month.name, levels = month.name)
df$intv <- ifelse(df$time >= int, 1, 0)
df$intv_trend <- c(rep(0, (int - 1)),
                   1:(length(unique(df$time)) - (int - 1)))
df <- df %>%
  mutate(lag_count = dplyr::lag(count))

创建模型并获得预测:

fit <- glm(
  count ~ month + time + intv + intv_trend +
    log(lag_count) + offset(log(at_risk)),
  family = "poisson",
  data = df
)

df$fit <- exp(c(NA, predict(fit)))

# Get predictions on the same scale as the data
df$fit2 = c(NA, predict(fit, type="response"))

# Add a grouping variable manually
df$group = rep(c("Control","Intervention"), c(72, 132 - 72))

情节:

ggplot(df, aes(x = time, y = fit2)) +
  geom_line() +
  geom_smooth(method="lm", se=FALSE, aes(colour=group)) +
  theme_bw() +
  labs(colour="")