Feature-complete version of the new drfit package that can make use

of the drc package in order to obtain confidence intervals for EDx values.

See ChangeLog for details.



git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@98 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc

3 files changed, 328 insertions, 65 deletions

diff --git a/R/drcfit.R b/R/drcfit.R
new file mode 100644
index 0000000..ff9e6bb
--- /dev/null
+++ b/R/drcfit.R
@@ -0,0 +1,292 @@
+drcfit <- function(data, chooseone=TRUE,
+        probit = TRUE, logit = FALSE, weibull = FALSE,
+        linlogit = FALSE, level = 0.95,
+        showED50 = FALSE,
+        EDx = NULL)
+{
+    if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data
+                                                            # with ok set to
+                                                            # "no fit"
+    substances <- levels(data$substance)
+
+    ri <- rix <- 0                  # ri is the index over the result rows
+                                    # rix is used later to check if any
+                                    # model result was appended
+    rsubstance <- array()           # the substance names in the results
+    rndl <- vector()                # number of dose levels
+    rn <- vector()                  # mean number of replicates 
+                                    # in each dose level
+    runit <- vector()               # vector of units for each result row
+    rlhd <- rlld <- vector()        # highest and lowest doses tested
+    mtype <- array()                # the modeltypes
+    sigma <- array()                # the standard deviation of the residuals
+    logED50 <- vector()
+    logED50low <- logED50high <- vector()
+    a <- b <- c <- vector()
+
+    models <- list()                  # a list containing the dose-response models
+
+    splitted <- split(data,data$substance)
+    for (i in substances) {
+        tmp <- splitted[[i]]
+        fit <- FALSE
+        if (length(tmp) != 0) {
+            unit <- levels(as.factor(as.vector(tmp$unit)))
+            message("\n",i,": Fitting data...\n")
+        } else {
+            unit <- ""
+            message("\n",i,": No data\n")
+        }
+        if (length(unit) > 1) {
+            message("More than one unit for substance ",i,", halting\n\n")
+            break
+        }
+        if (length(tmp$response) == 0) {
+            nodata = TRUE
+        } else {
+            nodata = FALSE
+        }
+        rix <- ri
+        if (nodata) {
+            n <- ndl <- 0
+        } else {
+            ndl <- length(levels(factor(tmp$dose)))
+            n <- length(tmp$response)
+            highestdose <- max(tmp$dose)
+            lowestdose <- min(tmp$dose)
+            lhd <- log10(highestdose)
+            lld <- log10(lowestdose)
+            responseathighestdose <- mean(subset(tmp,dose==highestdose)$response)
+            responseatlowestdose <- mean(subset(tmp,dose==lowestdose)$response)
+            if (responseathighestdose < 0.5) {
+                inactive <- FALSE
+                if (responseatlowestdose < 0.5) {
+                    active <- TRUE
+                } else {
+                    active <- FALSE
+                    if (linlogit)
+                    {
+                        m <- try(drm(response ~ dose, data = tmp, fct = BC.4(fixed = c(NA, 1, NA, NA))))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                mtype[[ri]] <- "linlogit"
+                                models[[ri]] <- m
+                                s <- summary(m)
+                                sigma[[ri]] <- s$rseMat[1, "rse"]
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                ED50 <- try(ED(m, 50, interval = "delta", display = FALSE))
+                                if (!inherits(ED50, "try-error")) {
+                                    logED50[[ri]] <- log10(ED50["1:50", "Estimate"])
+                                    logED50low[[ri]] <- log10(ED50["1:50", "Lower"])
+                                    logED50high[[ri]] <- log10(ED50["1:50", "Upper"])
+                                    if (logED50[[ri]] > rlhd[[ri]]) {
+                                        mtype[[ri]] <- "no fit"
+                                        logED50[[ri]] <- NA
+                                        logED50low[[ri]] <- NA
+                                        logED50high[[ri]] <- NA
+                                    }
+                                } else {
+                                    mtype[[ri]] <- "no ED50"
+                                }
+                                a[[ri]] <- coef(m)[[2]]
+                                b[[ri]] <- coef(m)[[1]]
+                                c[[ri]] <- coef(m)[[3]]
+                            }
+                        }
+                    }
+                    if (probit)
+                    {
+                        m <- try(drm(response ~ dose, data = tmp, 
+                                     fct = LN.2()))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                models[[ri]] <- m
+                                s <- summary(m)
+                                sigma[[ri]] <- s$rseMat[1, "rse"]
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- log10(coef(m)[[2]])
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "probit"
+                                    ED50 <- ED(m, 50, interval = "delta", display = FALSE)
+                                    logED50low[[ri]] <- log10(ED50["1:50", "Lower"])
+                                    logED50high[[ri]] <- log10(ED50["1:50", "Upper"])
+                                    a[[ri]] <- coef(m)[[2]]
+                                    b[[ri]] <- coef(m)[[1]]
+                                }
+                            }
+                        }
+                    }
+                    if (logit)
+                    {
+                        m <- try(drm(response ~ dose, data = tmp, fct = LL.2()))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                models[[ri]] <- m
+                                s <- summary(m)
+                                sigma[[ri]] <- s$rseMat[1, "rse"]
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- log10(coef(m)[[2]])
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "logit"
+                                    ED50 <- ED(m, 50, interval = "delta", display = FALSE)
+                                    logED50low[[ri]] <- log10(ED50["1:50", "Lower"])
+                                    logED50high[[ri]] <- log10(ED50["1:50", "Upper"])
+                                    a[[ri]] <- coef(m)[[2]]
+                                    b[[ri]] <- coef(m)[[1]]
+                                }
+                            }
+                        }
+
+                    }
+                    if (weibull)
+                    {
+                        m <- try(drm(response ~ dose, data = tmp, fct = W1.2()))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                models[[ri]] <- m
+                                s <- summary(m)
+                                sigma[[ri]] <- s$rseMat[1, "rse"]
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                c[[ri]] <- NA
+                                ED50 <- ED(m, 50, interval = "delta", display = FALSE)
+                                logED50[[ri]] <- log10(ED50["1:50", "Estimate"])
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "weibull"
+                                    logED50low[[ri]] <- log10(ED50["1:50", "Lower"])
+                                    logED50high[[ri]] <- log10(ED50["1:50", "Upper"])
+                                    a[[ri]] <- logED50[[ri]]
+                                    b[[ri]] <- coef(m)[[1]]
+                                }
+                            }
+                        }
+
+                    }
+                }
+            } else {
+                inactive <- TRUE
+            }
+        }
+        if (ri == rix) {          # if no entry was appended for this substance
+            ri <- ri + 1
+            rsubstance[[ri]] <- i
+            rndl[[ri]] <- ndl
+            rn[[ri]] <- n
+            if (nodata) {
+                rlld[[ri]] <- rlhd[[i]] <- NA
+                mtype[[ri]] <- "no data"
+                runit[[ri]] <- NA
+            } else {
+                rlld[[ri]] <- log10(lowestdose)
+                rlhd[[i]] <- log10(highestdose)
+                runit[[ri]] <- unit
+                if (inactive) {
+                    mtype[[ri]] <- "inactive"
+                } else {
+                    if (active) {
+                        mtype[[ri]] <- "active"
+                    } else {
+                        mtype[[ri]] <- "no fit"
+                    }
+                }
+            }
+            sigma[[ri]] <- NA
+            logED50[[ri]] <- NA
+            logED50low[[ri]] <- NA
+            logED50high[[ri]] <- NA
+            a[[ri]] <- NA
+            b[[ri]] <- NA
+            c[[ri]] <- NA
+        }
+    }
+
+    results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, 
+        logED50, logED50low, logED50high, runit, sigma, a, b)
+    lower_level_percent = paste(100 * (1 - level)/2, "%", sep = "")
+    upper_level_percent = paste(100 * (1 + level)/2, "%", sep = "")
+    names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50",
+        lower_level_percent, upper_level_percent,
+        "unit","sigma","a","b")
+
+    if (linlogit) {
+        results$c <- c
+    }
+
+    if (showED50) {
+        results[c("ED50", paste("ED50", c(lower_level_percent, upper_level_percent)))] <-
+          10^results[7:9]
+    }
+
+    if (!is.null(EDx)) {
+        for (row.i in 1:ri) {
+            m <- models[[row.i]]
+            mtype <- as.character(results[row.i, "mtype"])
+            if (mtype %in% c("probit", "logit", "weibull", "linlogit")) {
+                for (EDi in EDx) {
+                    EDx.drc = try(ED(m, EDi, interval = "delta", display = FALSE, level = level))
+                    if (!inherits(EDx.drc, "try-error")) {
+                        results[row.i, paste0("EDx", EDi)] <- EDx.drc[paste0("1:", EDi), "Estimate"]
+                        results[row.i, paste0("EDx", EDi, " ", lower_level_percent)] <- EDx.drc[paste0("1:", EDi), 
+                                                                                                "Lower"]
+                        results[row.i, paste0("EDx", EDi, " ", upper_level_percent)] <- EDx.drc[paste0("1:", EDi), 
+                                                                                                "Upper"]
+                    }
+                }
+            }
+        }
+    }
+
+    attr(results, "models") <- models
+    return(results)
+}
+# vim: set ts=4 sw=4 expandtab:
diff --git a/R/drfit.R b/R/drfit.R
index fc73632..ad7d16c 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -28,6 +28,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     logED50low <- logED50high <- vector()
     a <- b <- c <- vector()
 
+    models <- list()                  # a list containing the dose-response models
+
     splitted <- split(data,data$substance)
     for (i in substances) {
         tmp <- splitted[[i]]
@@ -79,6 +81,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
                             ri <- ri + 1
+                            models[[ri]] <- m
                             s <- summary(m)
                             sigma[[ri]] <- s$sigma
                             rsubstance[[ri]] <- i
@@ -122,6 +125,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                             if (!inherits(m, "try-error")) {
                                 fit <- TRUE
                                 ri <- ri + 1
+                                models[[ri]] <- m
                                 s <- summary(m)
                                 sigma[[ri]] <- s$sigma
                                 rsubstance[[ri]] <- i
@@ -165,6 +169,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                             if (!inherits(m, "try-error")) {
                                 fit <- TRUE
                                 ri <- ri + 1
+                                models[[ri]] <- m
                                 s <- summary(m)
                                 sigma[[ri]] <- s$sigma
                                 rsubstance[[ri]] <- i
@@ -206,6 +211,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                         if (chooseone==FALSE || fit==FALSE) {
                             if (!inherits(m, "try-error")) {
                                 ri <- ri + 1
+                                models[[ri]] <- m
                                 a[[ri]] <- coef(m)[["location"]]
                                 b[[ri]] <- coef(m)[["shape"]]
                                 sqrdev <- function(logdose) {
@@ -299,38 +305,12 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     }
     if (!is.null(EDx)) {
         for (row.i in 1:ri) {
-            logED50 <- results[row.i, "logED50"]
-            mtype <- as.character(results[row.i, "mtype"])
-            if (mtype == "probit") {
-                scale <- results[row.i, "b"]
-                drfunction <- function(x) pnorm(-x, -logED50, scale)
-            }
-            if (mtype == "logit") {
-                scale <- results[row.i, "b"]
-                drfunction <- function(x) plogis(-x, -logED50, scale)
-            }
-            if (mtype == "weibull") {
-                location <- results[row.i, "a"]
-                shape <- results[row.i, "b"]
-                drfunction <- function(x) pweibull(-x + location, shape)
-            }
-            if (mtype == "linlogit") {
-                drfunction <- function(x) {
-                       r <- linlogitf(10^x, 1, 
-                                      results[row.i, "c"], 
-                                      results[row.i, "logED50"],
-                                      results[row.i, "b"])
-                       # Do not return response values above 1 to avoid trapping
-                       # the optimization algorithm in a local minimum later on
-                       if (r < 1) return(r)
-                       else return(2 - 0.001 * x) # Start with 2 and decrease slowly to 
-                                                  # guide to the interesting part of the curve
-                }
-            }
-            if (mtype %in% c("probit", "logit", "weibull", "linlogit")) {
+            if (mtype[[row.i]] %in% c("probit", "logit", "weibull", "linlogit")) {
                 for (ED in EDx) {
-                    of <- function(x) abs(drfunction(x) - (1 - (ED/100)))
-                    
+                    of <- function(x) {
+                        abs(predict(models[[row.i]], data.frame(dose = 10^x)) - 
+                            (1 - (ED/100)))
+                    } 
                     # Search over interval starting an order of magnitude below 
                     # the lowest dose up to one order of magnitude above the
                     # highest dose
@@ -347,5 +327,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     }
 
     rownames(results) <- 1:ri
+    attr(results, "models") <- models
     return(results)
 }
+# vim: set ts=4 sw=4 expandtab:
diff --git a/R/drplot.R b/R/drplot.R
index 149337c..ded7aca 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -1,4 +1,4 @@
-if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance"))
+if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance", "mtype"))
 drplot <- function(drresults, data, 
         dtype = "std", alpha = 0.95, ctype = "none",
         path = "./", fileprefix = "drplot", overlay = FALSE,
@@ -9,7 +9,8 @@ drplot <- function(drresults, data,
         postscript = FALSE, pdf = FALSE, png = FALSE, 
         bw = TRUE,
         pointsize = 12,
-        colors = 1:8, ltys = 1:8, devoff=TRUE, lpos="topright")
+        colors = 1:8, ltys = 1:8, pchs = "auto",
+        devoff=TRUE, lpos="topright")
 {
     # Check if all data have the same unit
     unitlevels <- levels(as.factor(drresults$unit))
@@ -63,6 +64,7 @@ drplot <- function(drresults, data,
     # Prepare overlay plot if requested
     if (overlay)
     {
+        if (pchs == "auto") pchs = 1:length(dsubstances)
         if (postscript) {
             filename = paste(path,fileprefix,".eps",sep="")
             postscript(file=filename,
@@ -91,6 +93,7 @@ drplot <- function(drresults, data,
             frame.plot = frame.plot)
     } else {
         # If overlay plot is not requested, ask before showing multiple plots on the screen
+        if (pchs == "auto") pchs = rep(1, length(dsubstances))
         if (!postscript && !png && !pdf && length(dsubstances) > 1) {
             op <- par(ask=TRUE)
             on.exit(par(op))
@@ -111,6 +114,7 @@ drplot <- function(drresults, data,
             tmp <- splitted[[i]]
             if (length(tmp$response) != 0) {
                 color <- colors[[n]]
+                pch <- pchs[[n]]
                 # Prepare the single graphs if an overlay is not requested
                 if (!overlay)
                 {
@@ -141,7 +145,7 @@ drplot <- function(drresults, data,
                         axes = axes,
                         frame.plot = frame.plot)
                 }
-                if (!overlay) legend(lpos, i, lty = 1, col = color, inset=0.05)
+                if (!overlay) legend(lpos, i, lty = 1, col = color, pch = pch, inset=0.05)
                 tmp$dosefactor <- factor(tmp$dose)  # necessary because the old
                                                     # factor has all levels, not 
                                                     # only the ones tested with
@@ -160,7 +164,7 @@ drplot <- function(drresults, data,
                 # Plot the data, if requested
                 if (dtype != "none") {
                     if (dtype == "raw") {
-                        points(log10(tmp$dose),tmp$response,col=color)
+                        points(log10(tmp$dose), tmp$response, col = color, pch = pch)
                     } else {
                         splitresponses <- split(tmp$response,tmp$dosefactor)
                         means <- sapply(splitresponses,mean)
@@ -183,7 +187,7 @@ drplot <- function(drresults, data,
                     {
                         x <- log10(as.numeric(levels(tmp$dosefactor)))
                         segments(x,bottoms,x,tops,col=color)
-                        points(x,means,col=color)
+                        points(x, means, col = color, pch = pch)
                         smidge <- 0.05
                         segments(x - smidge,bottoms,x + smidge,bottoms,col=color)
                         segments(x - smidge,tops,x + smidge,tops,col=color)
@@ -192,37 +196,22 @@ drplot <- function(drresults, data,
 
                 # Plot the fits for this substance, if there are any
                 fits <- subset(drresults,Substance == i)
-                nf <-  length(fits$Substance)  # number of fits to plot for this substance
+                fit.rows <- rownames(fits)
+                nf <-  length(fit.rows)  # number of fits to plot for this substance
                 if (nf > 0) {
                     for (j in 1:nf)
                     {
-                        logED50 <- fits[j,"logED50"]
-                        mtype <- as.character(fits[j, "mtype"])
-                        if (mtype == "probit") {
-                            if (overlay) nl <- nl + 1 else nl = j
-                            lty <- ltys[nl]
-                            scale <- fits[j,"b"]
-                            plot(function(x) pnorm(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty)
-                        }
-                        if (mtype == "logit") {
-                            if (overlay) nl <- nl + 1 else nl = j
-                            lty <- ltys[nl]
-                            scale <- fits[j,"b"]
-                            plot(function(x) plogis(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty)
-                        }
-                        if (mtype == "weibull") {
-                            if (overlay) nl <- nl + 1 else nl = j
-                            lty <- ltys[nl]
-                            location <- fits[j,"a"]
-                            shape <- fits[j,"b"]
-                            plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty)
-                        }
-                        if (mtype == "linlogit") {
-                            if (overlay) nl <- nl + 1 else nl = j
-                            lty <- ltys[nl]
-                            plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logED50"],fits[j,"b"]),
-                                lld - 0.5, lhd + 2,
-                                add=TRUE, col=color, lty=lty)
+                        fit.row = fit.rows[j]
+                        if (overlay) nl <- nl + 1 else nl = j
+                        lty <- ltys[nl]
+                        if (drresults[[fit.row, "mtype"]] %in% c("probit", "logit", "weibull", "linlogit")) 
+                        {
+                            m <- attr(drresults, "models")[[as.numeric(fit.row)]]
+                            of <- function(x) {
+                                predict(m, data.frame(dose = 10^x))
+                            } 
+                            plot(of, lld - 0.5, lhd + 2, 
+                                 add = TRUE, col = color, lty = lty)
                         }
                     }
                 }
@@ -232,7 +221,7 @@ drplot <- function(drresults, data,
             }
         }
     }
-    if (overlay) legend(lpos, dsubstances, col = colors, lty = ltys, inset=0.05)
+    if (overlay) legend(lpos, dsubstances, col = colors, pch = pchs, lty = ltys, inset=0.05)
     if (overlay && (postscript || png || pdf)) {
         if (devoff) {
             dev.off()