我已经安装了一个 GLM，但未能使用 ggplot 绘制模型

Question

我制作了一个 GLM 并尝试使用以下代码使用 ggplot 绘制模型。我想我需要添加类型参数，这样我的模型就不会只产生一条水平的扁平线。但是我正在努力使用类型参数 as i 运行将其放入此错误消息中；

 Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = 
 object$xlevels) : 
 object is not a matrix

这是我的数据集和用于获取此数据的代码（如果有人知道此问题的修复方法）

我的数据（前 10 行）；

 aids
   cases quarter  date
 1      2       1 83.00
 2      6       2 83.25
 3     10       3 83.50
 4      8       4 83.75
 5     12       1 84.00
 6      9       2 84.25
 7     28       3 84.50
 8     28       4 84.75
 9     36       1 85.00
 10    32       2 85.25

我的代码用于创建模型和绘图

 model3 = glm(cases ~ date,
          data = aids,
          family = poisson(link='log'))


 #plotting the model (“link”, “response”, “terms”)
 plot_predictions <- function(model, type = 'response') {

   #make predictions
    preds <- predict(model, df)

   #plot
   df %>% 
    ggplot(aes(date, cases)) +
    geom_point() +
    geom_line(aes(date, preds), col = 'red') +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'))
 }

 plot_predictions(model3, aids)

dput输出

dput(head(aids, 10))
structure(list(cases = c(2, 6, 10, 8, 12, 9, 28, 28, 36, 32), 
quarter = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
2L), .Label = c("1", "2", "3", "4"), class = "factor"), date = 
c(83, 
83.25, 83.5, 83.75, 84, 84.25, 84.5, 84.75, 85, 85.25)), 
row.names = c(NA, 
10L), class = "data.frame")

Answer 1

这是一个将预测放在正确范围内的方法：

library(tidyverse)
aids <- structure(list(cases = c(2, 6, 10, 8, 12, 9, 28, 28, 36, 32), 
               quarter = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
                                     2L), .Label = c("1", "2", "3", "4"), class = "factor"), date = 
                 c(83, 
                   83.25, 83.5, 83.75, 84, 84.25, 84.5, 84.75, 85, 85.25)), 
          row.names = c(NA, 
                        10L), class = "data.frame")
model3 = glm(cases ~ date,
             data = aids,
             family = poisson(link='log'))


plot_predictions <- function(model, df, type = 'response') {
  require(tidyverse)
  #make predictions
  preds <- predict(model, df, type= type)
  
  #plot
  df %>% 
    ggplot(aes(date, cases)) +
    geom_point() +
    geom_line(aes(date, preds), col = 'red') +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'))
}

plot_predictions(model3, aids)

^{由 reprex package (v2.0.1)}

于 2022-05-25 创建

以下是您如何使用置信区间进行计算，尽管这仅适用于双变量模型。我已经包含了一个使用 ggeffects 包中的 ggpredict() 的选项，它可以在这种情况和其他情况下使用。这里的技巧是你必须在 link 尺度上进行预测，根据预测的标准误差做出置信区间，然后通过逆 link.

library(tidyverse)
aids <- structure(list(cases = c(2, 6, 10, 8, 12, 9, 28, 28, 36, 32), 
                       quarter = structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 
                                             2L), .Label = c("1", "2", "3", "4"), class = "factor"), date = 
                         c(83, 
                           83.25, 83.5, 83.75, 84, 84.25, 84.5, 84.75, 85, 85.25)), 
                  row.names = c(NA, 
                                10L), class = "data.frame")
model3 = glm(cases ~ date,
             data = aids,
             family = poisson(link='log'))


plot_predictions <- function(model, df, conf=.95, type = 'response') {
  require(tidyverse)
  #make predictions
  preds <- predict(model, df, type= "link", se.fit=TRUE)
  preds <- as.data.frame(preds[1:2])
  preds$x <- df$date
  preds <- preds %>% 
    mutate(lwr = fit - pnorm(1-(1-conf/2))*se.fit, 
           upr = fit + pnorm(1-(1-conf/2))*se.fit, 
           across(c(fit, lwr, upr), ~family(model)$linkinv(.x)))
  
  
  
  #plot
  ggplot(data=df, aes(date, cases)) +
    geom_ribbon(data=preds, aes(x=x, y=fit, ymin=lwr, ymax=upr), alpha=.25, fill="red", col="transparent") + 
    geom_line(data=preds, aes(x, fit), col = 'red') +
    geom_point() +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'))
}

plot_predictions(model3, aids)

这是 ggpredict() 选项：

library(ggeffects)
g <- ggpredict(model3, terms="date [all]")
plot(g, rawdata=TRUE)

^{由 reprex package (v2.0.1)}

于 2022-05-25 创建

Answer 2

你在这里有点重新发明轮子。这可以在 ggplot 中使用 geom_smooth

自动完成

ggplot(aids, aes(date, cases)) +
    geom_point() +
    geom_smooth(col = 'red', se= FALSE, method = glm,
                method.args = list(family = poisson(link = 'log'))) +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'))

您甚至可以通过省略 se = FALSE:

来显示标准错误

ggplot(aids, aes(date, cases)) +
    geom_segment(aes(xend = date, yend = 0), color = "deepskyblue4") +
    geom_point(size = 3) +
    geom_smooth(col = 'red3', fill = "red3", method = glm, alpha = 0.1,
                method.args = list(family = poisson(link = 'log'))) +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme_minimal(base_size = 16) +
    theme(plot.title = element_text(hjust = 0.5, face = 'bold'))

Answer 3

您也可以使用 broom::augment 将模型的预测添加到数据中并绘制：

model3 %>%
    broom::augment() %>%
    ggplot(aes(date, cases)) +
    geom_point() +
    geom_line(aes(date, exp(.fitted)), col = 'red') +
    ggtitle("Model 2 - Poisson GLM predicting cases") +
    theme(plot.title = element_text(hjust = 0.5, size = 12, face = 'bold'))

我已经安装了一个 GLM，但未能使用 ggplot 绘制模型

I have fitted a GLM and failing to plot the model using ggplot

r

ggplot2

glm