2

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

16% - 84% 之间的范围应该在 LM 和 GLM 之间排列，引文：第 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)))
4

1 回答 1

6

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

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

在此处输入图像描述

于 2017-02-27T04:36:12.137 回答