绘制相同数据的GLM和LM

Question

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

LM和GLM之间的范围应在16%-84%之间,引用:第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是

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))