Plotting dr4pl dose response curves, and how to integrate them with ggplot2?

Question

I am trying to set up a high-throughput way of plotting out dose response curves from a large screening experiment. Prism obviously has the easiest way of doing dose-response curves well, but I can't copy and paste this much data.

Since CRAN removed drc, the package dr4pl seems the way to go, but there is very little instruction available yet.

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

The above will generate the type of curve I need using dr4pl, and get me the IC50. but how would I plot several datasets/curves on one plot?

Ideally I'd rather plot the data using ggplot2: plot+geom_point() and add in the dose response line by using the dr4pl summary as a +stat_smooth() model, if that makes sense? But I don't know how to do this.

Any help would be appreciated

Answer 1

I can get most of the way but not all the way. The main step is to write a predict() method for dr4pl objects:

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

I included a slightly hacky way to compute the confidence intervals via the delta method - this might not be too reliable (bootstrapping would be better ...)

It works OK (sort of) for your data (changed the name to dd because it's sometimes dicey to name the data 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")

• the confidence intervals are terrible, turn them off with se=FALSE
• dr4pl puts the x-axis on a log10-scale by default, but the standard scale_x_log10() screws this up because it is applied before the fitting and prediction, so I use coord_trans(x="log10") instead.
• However, coord_trans() doesn't play so nicely if the axes are on a very broad logarithmic scale - I tried the example above with the sample_data_1 data from the package and it didn't work.

But I'm afraid I've spent enough time on this for now.

It would more robust to use the predict method above to generate the values you want, over the range you want, separately, and then use geom_line() + geom_ribbon() to add the information to the plot ....

If you're willing to fit the model first (outside geom_smooth) you can do this (this is using sample_data_1 from dr4pl package - it's from the first example in ?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()

which is less sensitive to the order of scaling/unscaling the x axis.

---

Improved but slow bootstrap 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")
)