r - 绘制相同数据的 GLM 和 LM

Question

我想绘制相同数据的线性模型 (LM) 和非线性 (GLM) 模型。

16% - 84% 之间的范围应该在 LM 和 GLM 之间，Citation: section 3.5

我包含了更完整的代码块，因为我不确定应该在哪一点尝试削减线性模型。或者在哪一点我搞砸了-我认为是线性模型。

下面的代码会产生下图:

我的目标(取自之前的引用链接)。

这是我的数据:

mydata3 <- structure(list(
               dose = c(0, 0, 0, 3, 3, 3, 7.5, 7.5, 7.5, 10,     10, 10, 25, 25, 25, 50, 50, 50), 
               total = c(25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L), 
               affected = c(1, 0, 1.2, 2.8, 4.8, 9, 2.8, 12.8, 8.6, 4.8, 4.4, 10.2, 6, 20, 14, 12.8, 23.4, 21.6), 
               probability = c(0.04, 0, 0.048, 0.112, 0.192, 0.36, 0.112, 0.512, 0.344, 0.192, 0.176, 0.408, 0.24, 0.8, 0.56, 0.512, 0.936, 0.864)), 
               .Names = c("dose", "total", "affected", "probability"), 
               row.names = c(NA, -18L), 
               class = "data.frame")

我的脚本:

#load libraries

library(ggplot2)
library(drc)  # glm model
library(plyr) # rename function
library(scales) #log plot scale

#Creating linear model
mod_linear <- lm(probability ~ (dose), weights = total, data = mydata3)

#Creating data.frame: note values 3 and 120 refer to 16% and 84% response in sigmoidal plot
line_df <-expand.grid(dose=exp(seq(log(3),log(120),length=200)))

#Extracting values from linear model
p_line_df <- as.data.frame(cbind(dose = line_df,
                               predict(mod_linear, newdata=data.frame(dose = line_df),
                                       interval="confidence",level=0.95)))

#Renaming linear df columns
p_line_df <-rename(p_line_df, c("fit"="probability"))
p_line_df <-rename(p_line_df, c("lwr"="Lower"))
p_line_df <-rename(p_line_df, c("upr"="Upper"))
p_line_df$model <-"Linear"


#Create sigmoidal dose-response curve using drc package
    mod3 <- drm(probability ~ (dose), weights = total, data = mydata3, type ="binomial", fct=LL.2(names=c("Slope:b","ED50:e")))

    #data frame for ggplot2 
    base_DF_3 <-expand.grid(dose=exp(seq(log(1.0000001),log(10000),length=200)))

    #extract data from model
    p_df3 <- as.data.frame(cbind(dose = base_DF_3,
                                 predict(mod3, newdata=data.frame(dose = base_DF_3),
                                         interval="confidence", level=.95)))

#renaming columns
p_df3 <-rename(p_df3, c("Prediction"="probability"))
p_df3$model <-"Sigmoidal" 

#combining Both DataFames
 p_df_all <- rbind(p_df3, p_line_df)

#plotting
ggplot(p_df_all, aes(x=dose,y=probability, group=model))+
    geom_line(aes(x=dose,y=probability,group=model,linetype=model),show.legend = TRUE)+
  scale_x_log10(breaks = c(0.000001, 10^(0:10)),labels = c(0, math_format()(0:10)))

Answer 1

查看您提供的引用资料，作者描述的是使用线性模型来近似(S 形)逻辑函数的中心部分。实现这一点的线性模型是一条直线，它穿过逻辑曲线的拐点，并且与该拐点处的逻辑函数具有相同的斜率。我们可以使用一些基本的代数和微积分来解决这个问题。

从?LL.2，我们看到drm拟合的logistic函数的形式是

f(x) = 1 / {1 + exp(b(log(x) - log(e)))}

我们可以得到这个方程中的系数值

b = mod3$coefficients[1]
e = mod3$coefficients[2]

现在，通过微分，逻辑函数的斜率由下式给出

dy/dx = -(b * exp((log(x)-log(e))*b)) / (1+exp((log(x)-log(e))*b))^2

在拐点处，剂量 (x) 等于系数 e，因此拐点处的斜率(大大)简化为

sl50 = -b/4

由于我们还知道拐点出现在 probability = 0.5 和 dose = e 的点，我们可以构造直线(在对数变换坐标中)像这样:

linear_probability = sl50 * (log(p_df3$dose) - log(e)) + 0.5

现在，一起绘制逻辑函数和线性函数:

p_df3_lin = p_df3
p_df3_lin$model = 'linear'
p_df3_lin$probability = linear_probability

p_df_all <- rbind(p_df3, p_df3_lin)

ggplot(p_df_all, aes(x=dose,y=probability, group=model))+
  geom_line(aes(x=dose,y=probability,group=model,linetype=model),show.legend = TRUE)+
  scale_x_log10(breaks = c(0.000001, 10^(0:10)),labels = c(0, math_format()(0:10))) +
  scale_y_continuous(limits = c(0,1))