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

Question

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



摘自:Lindstrand A、Bennet R、Galanis I 等。引入肺炎球菌结合疫苗后鼻窦炎和肺炎住院。儿科。 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 lineplot 中创建的线条)

有持续的长期趋势 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="")