使用偏移变量时如何绘制负二项式回归图（拟合值和 95% 置信带）？

Question

我正在尝试在 R 中绘制负二项式模型图，但在使用偏移变量时我什至无法提取拟合值的置信限。有没有什么方法可以在不尝试手动获取 95% CI 的情况下绘制绘图？

y <- c(602,38,616,256,21,723,245,176,89,1614,31,27,313,251,345)
x <- c(31,35,21,30,37,26,45,21,74,27,37,37,31,37,25)
offset_1 <-c(72,50,31,30,16,25,75,16,78,40,68,25,71,52,17) 
newdata <- data.frame(y,x,offset_1)
nb.model <- MASS::glm.nb(y ~ x + offset(log(offset_1)), data=newdata)
summary(nb.model)
predict(nb.model, type="response")/newdata$offset_1

Answer 1

以下是如何使用

emmeans

+

ggplot

来完成此操作（

ggeffects

/

sjPlot

可能会进一步自动化该过程，但我喜欢自己完成的控制）。

emmeans

中的偏移处理在here进行了讨论。

library(emmeans)
## generated predicted values / CIs
ee <- emmeans(nb.model, ~x,
       at = list(x = seq(20,75, length.out = 51)), 
       type = "response", offset = log(1))
library(ggplot2)
ggplot(as.data.frame(ee), aes(x, response)) +
    geom_line() +
    geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL), 
                colour = NA, alpha = 0.2) +
    ## add observed data
    geom_point(data = newdata, aes(y=y/offset_1))

Answer 2

我不知道有任何直接方法（除了 ggeffects 包）可以从具有偏移量的负二项式模型中绘制置信区间，尽管我们可以使用 broom 返回拟合的标准误差并使用以下方法计算置信区间正态近似（Wald 置信区间）。

library(broom)
library(dplyr)
library(ggplot2)
conf_level <- 0.95
cval <- qnorm(1 - ( (1 - conf_level) / 2))
augment(nb.model, se_fit = TRUE) %>%
  mutate(lcl = .fitted - cval * .se.fit,
         ucl = .fitted + cval * .se.fit) %>%
  mutate(across(c(.fitted, lcl, ucl), exp)) %>%
  mutate(across(c(y, .fitted, lcl, ucl), ~ .x / exp(offset(log(offset_1))))) %>%
  ggplot(aes(x = x, y = y)) + 
  geom_point() +
  geom_line(aes(y = lcl), lty = 2) + 
  geom_line(aes(y = ucl), lty = 2)

Answer 3

正如 Ben 建议的那样，这也是使用 ggeffects 包的“自动化过程”。该包警告用户由于涉及转换后的偏移量，预测可能不准确。

library(ggeffects)
newdata <- data.frame(
  y = c(602, 38, 616, 256, 21, 723, 245, 176, 89, 1614, 31, 27, 313, 251, 345),
  x = c(31, 35, 21, 30, 37, 26, 45, 21, 74, 27, 37, 37, 31, 37, 25),
  offset_1 = c(72, 50, 31, 30, 16, 25, 75, 16, 78, 40, 68, 25, 71, 52, 17)
)

m <- MASS::glm.nb(y ~ x + offset(log(offset_1)), data = newdata)
pr <- ggemmeans(m, "x") # use "x [20:75]" for smoother plots
#> Model uses a transformed offset term. Predictions may not be correct.
#>   Please apply transformation of offset term to the data before fitting
#>   the model and use `offset(offset_1)` in the model formula.
#>   You could also fix the offset term using the `condition` argument, e.g.
#>   `condition = c(offset_1 = 1)`.
plot(pr, add.data = TRUE)
#> Data points may overlap. Use the `jitter` argument to add some amount of
#>   random variation to the location of data points and avoid overplotting.

您可以使用

condition

参数将偏移项固定为特定值。但是，对于绘图，原始数据未进行调整，因此在使用当前 CRAN 版本的 ggeffects 包时，预测和原始数据无法正确匹配。但这在 ggeffects 的最新版本中已得到修复，您现在可以重现 Ben 的帖子中显示的相同结果：

pr <- ggemmeans(m, "x [20:75]", condition = c(offset_1 = 1))
plot(pr, add.data = TRUE)
#> Data points may overlap. Use the `jitter` argument to add some amount of
#>   random variation to the location of data points and avoid overplotting.

^{创建于 2023-08-30，使用 reprex v2.0.2}

要安装最新的软件包版本，请使用

ggeffects::install_latest()

，它会安装最新的 GitHub 版本（来自类似 CRAN 的服务器，因此您不需要从源安装软件包时所需的任何工具）。