Why single confidence interval when plotting drm model with multiple curves?

918 Views Asked by At

I'm doing nonlinear regression in R using drm function from drc package. The drm function takes a curveid argument, which makes drm to fit multiple curves and store the results in single model object. Next, I plot the curves and confidence intervals. I have a problem, however, when the model object contains multiple curves, the curves appear fine, but confidence interval is often plotted only for the first curve (though by the alpha level of it, it seems as if it is done several times).

Starting from the examples provided by drm, I found I get the desired behaviour, all confidence intervals appear, when I change the name of a variable fed into curveid argument—odd enough, only "CURVE" works, nothing else, even "curve" does not work (reproducible example given below). This made me think that maybe there is a bug in plot.drc (S3 method for class 'drc') so that "CURVE" is hardcoded in there. I eyeballed at the source code, but can't tell if it is true since I'm a beginner R user with little programming experience.

Most likely there is no bug and I am just missing something important.

library(drc)

# create some data
df <- data.frame(
  x=rep(c(0.003, 0.01, 0.03, 1, 3, 10, 100),2),
  y=c(3,3,3,1.5,-2,-3,-3.2, 3.5,3.5,3.2,1,-2.5,-2.8,-2.8),
  CURVE=rep(1:2, each=7)
)

# working as it should
mod1 <- drm(y~x, curveid=CURVE, data=df)
plot(mod1, type="confidence", main="working as it should")

# not working
names(df)[3] <- "curve"
mod2 <- drm(y~x, curveid=curve, data=df)
plot(mod2, type="confidence", main="not working")

# working again
names(df)[3] <- "CURVE"
mod3 <- drm(y~x, curveid=CURVE, data=df)
plot(mod3, type="confidence", main="working again")

Resulting three plots as an image

1

There are 1 best solutions below

1
On BEST ANSWER

Vallo,

I can confirm the behavior you describe. Only the column name 'CURVE' works. All other names I tired do not work.

I think i have tracked this behavior to the behavior (bug?) of the ciFct subfunction in the plot.drc function.

   ciFct <- function(level, ...)
    {
        newdata <- data.frame(DOSE=dosePts, CURVE=rep(level, length(dosePts)))
        predictMat <- predict(object, 
                              newdata=newdata,
                              interval = "confidence",
                              level=confidence.level)

        x <- c(dosePts, rev(dosePts))
        y <- c(predictMat[,"Lower"], rev(predictMat[,"Upper"]))
        polygon(x,y, border=NA, ...)
    }

This is the function that calls predict() and generates the confidence bands. Note that it uses the hard-coded value 'CURVE' to build newdata on line 3. This only works properly if you used 'CURVE' in your original dataframe -- Otherwise it fails to work properly.

To fix this, you can of course always choose the column name 'CURVE'. Otherwise you will need to replace the plot.drc function with a patched version of the code above. The code can be fixed by changing the ciFit function to the following (note, only one line has been added):

  ciFct <- function(level, ...)
  {
    newdata <- data.frame(DOSE=dosePts, CURVE=rep(level, length(dosePts)))
    names(newdata)[2] <- object$curveVarNam # this is a new line to fix the bug
    predictMat <- predict(object, 
                          newdata=newdata,
                          interval = "confidence",
                          level=confidence.level)

    x <- c(dosePts, rev(dosePts))
    y <- c(predictMat[,"Lower"], rev(predictMat[,"Upper"]))
    polygon(x,y, border=NA, ...)
  }