从指定的 x 值显示 geom_smooth() 趋势线

Question

假设一个数据集包含多个时间段和多个组的计数数据，格式如下：

set.seed(123)
df <- data.frame(group = as.factor(rep(1:3, each = 50)),
                 week = rep(1:50, 3),
                 rate = c(round(700 - rnorm(50, 100, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0)))

    group week rate
1       1    1  604
2       1    2  598
3       1    3  578
4       1    4  591
5       1    5  589
6       1    6  571
7       1    7  581
8       1    8  597
9       1    9  589
10      1   10  584

我有兴趣为每组拟合基于模型的趋势线，但是，我希望仅从某个 x 值显示此趋势线。使用所有数据点可视化趋势线（需要ggplot2）：

df %>%
 ggplot(aes(x = week,
            y = rate,
            group = group,
            lty = group)) + 
 geom_line() +
 geom_point() +
 geom_smooth(method = "glm", 
             method.args = list(family = "quasipoisson"),
             se = FALSE) 

或者根据特定范围的值拟合模型（需要ggplot2 和dplyr）：

df %>%
 group_by(group) %>%
 mutate(rate2 = ifelse(week < 35, NA, rate)) %>%
 ggplot(aes(x = week,
            y = rate,
            group = group,
            lty = group)) + 
 geom_line() +
 geom_point() +
 geom_smooth(aes(y = rate2),
             method = "glm", 
             method.args = list(family = "quasipoisson"),
             se = FALSE)

但是，我找不到使用所有数据拟合模型的方法，而是仅显示特定 x 值（比如 35+）的趋势线。因此，我本质上想要为图一计算的趋势线，但根据第二个图显示它，使用ggplot2，理想情况下只有一条管道。

Answer 1

我去看了@tjebo提到的after_stat函数。看看以下内容是否适合您？

df %>%
  ggplot(aes(x = week,
             y = rate,
             lty = group)) + 
  geom_line() +
  geom_point() +
  geom_smooth(method = "glm", 
              aes(group = after_stat(interaction(group, x > 35)),
                  colour = after_scale(alpha(colour, as.numeric(x > 35)))),
              method.args = list(family = "quasipoisson"),
              se = F)

这通过将与每条线相关联的点分成两组来工作，即 x 35 区域中的点，因为线的颜色不应变化，并为每个新组。因此，只有 x > 35 区域中的线条可见。

使用时，代码会触发警告，指出 after_scale 修改未应用于图例。不过我认为这不是问题，因为无论如何我们都不需要它出现在图例中。

Answer 2

如果您可以容忍警告，则可以使用 stage() 与示例代码相差 1 行来解决此问题。

library(tidyverse)

set.seed(123)
df <- data.frame(group = as.factor(rep(1:3, each = 50)),
                 week = rep(1:50, 3),
                 rate = c(round(700 - rnorm(50, 100, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0)))

df %>%
  ggplot(aes(x = week,
             y = rate,
             group = group,
             lty = group)) + 
  geom_line() +
  geom_point() +
  geom_smooth(method = "glm", 
              method.args = list(family = "quasipoisson"),
              aes(x = stage(week, after_stat = ifelse(x > 35, x, NA))),
              se = FALSE) 
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 165 rows containing missing values (geom_smooth).

Answer 3

一种方法是在ggplot 之外构造拟合值，以便您可以控制它们：

df$fit <- glm(rate ~ week + group, data = df, family = "quasipoisson")$fitted.values

library(dplyr)
library(ggplot2)

ggplot(df, aes(x = week, group = group, lty = group)) + 
  geom_line(aes(y = rate)) +
  geom_point(aes(y = rate)) +
  geom_line(data = df %>% filter(week >= 35), aes(y = fit), color = "blue", size = 1.25)

Answer 4

1. 我不确定在时间序列中使用线性模型是否一般正确。时间序列的全部意义在于它们需要特定的统计数据，因为它们具有预期的自相关。您可能需要类似平均滚动模型的东西。

2. 我不确定您的可视化是否会让人感到困惑，更危险的是，是否具有误导性。

此外，还有一个有趣的问题。我认为新的after_stat 可能会有所帮助，但我无法让它发挥作用。

所以，这里有一个快速破解。更改 geom-s 的顺序并在其间绘制一个矩形。我正在厚颜无耻地使用不同的主题，但如果你真的想使用theme_grey()，你也可以伪造轴线。

library(tidyverse)
set.seed(123)
df <- data.frame(group = as.factor(rep(1:3, each = 50)),
                 week = rep(1:50, 3),
                 rate = c(round(700 - rnorm(50, 100, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0),
                          round(1000 - rnorm(50, 200, 10) - 1:50 * 2, 0)))
df %>%
  ggplot(aes(x = week, y = rate, group = group, lty = group)) + 
  stat_smooth(se = FALSE) +
  geom_rect(xmin = -Inf, xmax = 35, ymin = -Inf, ymax = Inf, 
            fill = "white") +
  geom_line() +
  geom_point() +
  theme_classic()
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

^{由reprex package (v1.0.0) 于 2021 年 2 月 9 日创建}

附：我已经删除了代码中一些不必要的部分来重现这一点，比如模型规格。

Answer 5

您可以使用ggplot_build 来获取情节的结构：

p <- ggplot(df, aes(x = week,
                        y = rate,
                        group = group,
                        lty = group)) + 
                  geom_line() +
                  geom_point() +
                  geom_smooth(method = "glm", 
                              method.args = list(family = "quasipoisson"),
                              se = FALSE) 
p_build <- ggplot_build(p)

然后您可以修改内部数据，这里是数据列表的第三个元素 (geom_smooth)：

p_build$data[[3]]$x <- sapply(p_build$data[[3]]$x,function(x) {ifelse(x<35,NA,x)})

并使用ggplot_gtable 重新生成绘图（lm 计算仍适用于整个数据集）：

plot(ggplot_gtable(p_build))