绘制dr4pl剂量响应曲线，以及如何与ggplot2集成？

Question

我试图建立一种高通量的方法，从一个大型筛选实验中绘制出剂量反应曲线。棱镜显然有最简单的方法做剂量-反应曲线很好，但我不能复制和粘贴这么多的数据。

由于CRAN删除了drc，所以包dr4pl似乎是可行的，但目前可用的指令非常少。

#make data frame 
dose <- c("0.078125","0.156250","0.312500","0.625000","1.250000","2.500000","5.000000","10.000000","20.000000")
POC<-c("1.05637425", "0.87380081", "0.79171200", "0.83166848", "0.77361290", "0.35199288", "0.19404609",  "0.09079221", "0.09850658")
data<-data.frame(dose, POC)


#use the dr4pl pakcage to calculate curve and IC50 etc
model<-dr4pl(POC~dose, data)

summary.model <- summary(model)

summary.model$coefficients

#plot this
plot(dr4pl(POC~dose, data=data))

​

​

以上将生成我需要使用dr4pl的曲线类型，并给我IC50。但是，我如何在一幅图上绘制多个数据集/曲线？

理想情况下，我更愿意用ggplot2：plot+geom_point()绘制数据，并将dr4pl摘要作为+stat_smooth()模型添加到剂量响应行中，如果这有意义的话？但我不知道该怎么做。

如能提供任何帮助，将不胜感激。

Answer 1

我可以得到大部分的方式，但不是所有的方式。主要步骤是为predict()对象编写一个dr4pl方法：

predict.dr4pl <- function (object, newdata=NULL, se.fit=FALSE, level, interval) {
    xseq <- if (is.null(newdata)) object$data$Dose else newdata$x
    pred <- MeanResponse(xseq, object$parameters)
    if (!se.fit) {
        return(pred)
    }
    qq <- qnorm((1+level)/2)
    se <- sapply(xseq,
                 function(x) car::deltaMethod(object, 
    "UpperLimit + (LowerLimit - UpperLimit)/(1 + (x/IC50)^Slope)")[["Estimate"]])
    return(list(fit=data.frame(fit=pred,lwr=pred-qq*se,
                                upr=pred+qq*se), se.fit=se))
}

我包括了一种通过delta方法计算置信区间的有点麻烦的方法--这可能不太可靠(引导可能更好.)

它对您的数据工作正常(某种程度上)(将名称更改为dd，因为有时命名数据data (fortunes::fortune("dog"))很危险)。

dd <- data.frame(dose = c(0.078125,0.156250,0.312500,0.625000,1.25,
                          2.50,5.0,10.0,20.0),
                 POC = c(1.05637425, 0.87380081, 0.79171200,
                         0.83166848, 0.77361290, 0.35199288,
                         0.19404609,  0.09079221, 0.09850658))

library(dr4pl)
ggplot(dd, aes(dose,POC)) + geom_point() +
     geom_smooth(method="dr4pl",se=TRUE) + coord_trans(x="log10")

​

​

se=FALSE

• dr4pl将x轴默认设置为10比例尺，但是标准的scale_x_log10()错误了，因为它是在拟合和预测之前应用的，所以如果轴是在非常宽的对数尺度上，那么使用coord_trans(x="log10")

• coord_trans()就不会发挥那么好的作用--我尝试了上面的例子和软件包中的sample_data_1数据，但是没有起作用。

但恐怕我已经花了足够的时间在这上面了。

如果使用上面的predict方法分别在所需范围内生成所需的值，然后使用geom_line() + geom_ribbon()将信息添加到绘图中，则会更加健壮.

如果您愿意首先符合模型(geom_smooth之外)，您可以这样做(这是使用来自dr4pl包的sample_data_1 -它来自于?dr4pl中的第一个示例)

model2 <- dr4pl(dose = sample_data_1$Dose,
           response = sample_data_1$Response)

ggplot(sample_data_1, aes(Dose,Response)) + geom_point() +
    stat_function(fun=function(x) predict(model2,newdata=data.frame(x=x))) +
    scale_x_log10()

它对x轴的缩放/不缩放顺序不那么敏感。

改进的但缓慢的引导CIs：

predictdf.dr4pl <- function (model, xseq, se, level, nboot=200) {
    pred <- MeanResponse(xseq, model$parameters)
    if (!se) {
        return(base::data.frame(x=xseq, y=pred))
    }
    ## bootstrap residuals
    pred0 <- MeanResponse(model$data$Dose, model$parameters)
    res <- pred0-model$data$Response
    bootres <- matrix(nrow=length(xseq), ncol=nboot)
    pb <- txtProgressBar(max=nboot,style=3)
    for (i in seq(nboot)) {
        setTxtProgressBar(pb,i)
        mboot <- dr4pl(model$data$Dose,
                       pred0 + sample(res, size=length(pred0),
                                      replace=TRUE))
        bootres[,i] <- MeanResponse(xseq, mboot$parameters)
    }
    fit <- data.frame(x = xseq,
                      y=pred,
                      ymin=apply(bootres,1,quantile,(1-level)/2),
                      ymax=apply(bootres,1,quantile,(1+level)/2))
    return(fit)
}

print(ggplot(dd, aes(dose,POC))
    + geom_point()
    + geom_smooth(method="dr4pl",se=TRUE) + coord_trans(x="log10")
)

​

​