使用 ggplot 和 glm 绘制剂量反应曲线的问题

Question

我目前正在尝试使用 glm 创建剂量反应曲线。我可以使用 glm 中的 bionominal family 和 probit 函数创建曲线，但想使用 ggplot 而不是 R 的基本绘图函数来绘制曲线。将基本图与 ggplot 进行比较时，ggplot 生成的曲线不正确，我不确定如何使其与基本图相同。另外，使用 ggplot 绘制曲线时，置信区间不正确。谢谢你的帮助。

library(ggplot2)
library(Hmisc)
library(plyr)
library(MASS)



create dataframe: 
#1) column is exposure concentration
#2) column is the total number of organism died over 12 h of exposure to the 
    corresponding concentration 
#3) column is the total number that survived over 12 h to the corresponding 
    concentration
#4) column is the total number of organism exposed to the corresponding 
    concentration
#5) fifth is the percentage of organism that survived exposure at the 
    corresponding concentration 

 conc <- c(0.02, 0.45, 0.46, 0.50, 0.78, 0.80, 0.80, 0.92, 0.93, 1.00, 1.16, 
   1.17, 1.17, 1.48,1.51, 1.55, 1.88, 1.90, 2.02)

 dead <- c(0, 0,  0,  0,  0,  0,  0,  0,  0,  1,  7, 11, 4, 14, 14, 12, 12, 18, 17)

 survive <- c(15, 16, 15, 15, 15, 16, 14, 14, 10, 15, 12,  5, 12,  0,  1,  3,  0,  0,  0)

 total <- c(15, 16, 15, 15, 15, 16, 14, 14, 10, 16, 19, 16, 16, 14, 15, 15, 12, 18, 17)

 perc <- c(1.00, 1.00, 1.00, 1.00, 1.00,1.00, 1.00, 1.00, 1.00, 0.94,0.63, 
      0.31,0.75,0.00, 0.07, 0.20, 0.00, 0.00,0.00)

 data<-data.frame(conc,dead,survive,total,perc)
 head(data)
 attach(data)
 #create matrix of dead and survival
 y = cbind(dead,survive)

 #create binomial glm (probit model)
 model.results = glm(data = data, y ~ conc,binomial(link="probit"))
 summary(model.results)



 #use function from MASS to calculate LC
 dose.p(model.results,p=0.5)
 dose.p(model.results,p=c(0.1,0.25,0.5,0.99))

 #plot curve 
 plot(conc,(survive/(survive+dead)), ylab = "Percent Survival", 
 xlab="Concentration ")


 #To make function use the Estimate parameters from the binomial glm 
  used above
  logisticline <- function(z) {eta = -6.7421 + 5.4468 * z;1 / (1 + 
  exp(eta))}
  x <- seq(0,200.02,0.01)
  lines(x,logisticline(x),new = TRUE)




  #plot using ggplot

  ggplot(data, aes(x = conc, y = perc)) +
  geom_point() +
  geom_smooth(method="glm",method.args = list(family = "binomial"))

Answer 1

您可以使用 ggplot2 绘制拟合线，方法是根据模型进行预测或直接使用 geom_smooth 拟合模型。要执行后者，您需要以死亡比例作为响应变量来拟合模型，以 total 作为权重，而不是使用成功和失败矩阵作为响应变量。

使用glm，用比例加权重拟合模型如下：

# Calculate proportion
data$prop = with(data, dead/total)

# create binomial glm (probit model)
model.results2 = glm(data = data, prop ~ conc, 
                    family = binomial(link="probit"), weights = total)

您可以使用您拥有的数据集进行预测，或者为了制作更平滑的线条，您可以创建一个新数据集进行预测，该数据集具有更多 conc 值，就像您所做的那样。

preddat = data.frame(conc = seq(0, 2.02, .01) )

现在您可以通过predict 从模型中进行预测，使用这个data.frame 作为newdata。如果您使用type = "response"，您将通过反向链接获得数据规模的预测。因为你拟合了一个概率模型，所以这将使用逆概率。在您的示例中，您使用逆 logit 进行预测。

# Predictions with inverse probit
preddat$pred = predict(model.results2, newdata = preddat, type = "response")
# Predictions with inverse logit (?)
preddat$pred2 = plogis( predict(model.results2, newdata = preddat) )

要在ggplot 中拟合概率模型，您需要将比例用作y 变量和weight = total。在这里，我添加了来自模型预测的线条，因此您可以看到适合 ggplot 的概率模型给出与拟合概率模型相同的估计线。使用逆对数会给你一些不同的东西，这并不奇怪。

ggplot(data, aes(conc, prop) ) +
     geom_smooth(method = "glm", method.args = list(family = binomial(link = "probit") ), 
                 aes(weight = total, color = "geom_smooth line"), se = FALSE) +
     geom_line(data = preddat, aes(y = pred, color = "Inverse probit") ) +
     geom_line(data = preddat, aes(y = pred2, color = "Inverse logit" ) )