Several changes in this revision:

- Addition of the weibull fit
- Renaming of lognorm to probit, logis to logit, and linearlogis to linlogit
- Fix of a bug that prevented the logis fit to work in recent version
- Different way of listing parameters in the output.
- Documentation of changes



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

1 files changed, 87 insertions, 48 deletions

diff --git a/R/drfit.R b/R/drfit.R
index 38f8e75..603cf94 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -31,14 +31,14 @@ drdata <- function(substances, experimentator = "%", db = "cytotox",
     return(data)
 }
     
-linearlogisf <- function(x,k,f,mu,b)
+linlogitf <- function(x,k,f,mu,b)
 {
     k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))
 }
 
 drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
-        lognorm = TRUE, logis = FALSE,
-        linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA,
+        probit = TRUE, logit = FALSE, weibull = FALSE,
+        linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
         b0 = 2, f0 = 0)
 {
     if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data where ok 
@@ -56,9 +56,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
     sigma <- array()                      # the standard deviation of the residuals
     logEC50 <- vector()
     stderrlogEC50 <- vector()
-    slope <- vector()
-    b <- vector()
-    f <- vector()
+    a <- b <- c <- vector()
 
     splitted <- split(data,data$substance)
     for (i in substances) {
@@ -95,10 +93,10 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
             if (responseathighestdose < 0.5) {
                 inactive <- FALSE
 
-                if (linearlogis && 
-                    length(subset(linearlogisWrong,linearlogisWrong == i))==0 &&
+                if (linlogit && 
+                    length(subset(linlogitWrong,linlogitWrong == i))==0 &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b),
+                    m <- try(nls(response ~ linlogitf(dose,1,f,logEC50,b),
                             data=tmp,
                             start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))
                     if (!inherits(m, "try-error")) {
@@ -111,27 +109,28 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                         runit[[ri]] <- unit
                         rlld[[ri]] <- log10(lowestdose)
                         rlhd[[ri]] <- log10(highestdose)
-                        mtype[[ri]] <- "linearlogis"
+                        mtype[[ri]] <- "linlogit"
                         logEC50[[ri]] <- coef(m)[["logEC50"]]
-                        slope[[ri]] <- NA
                         if (logEC50[[ri]] > rlhd[[ri]]) {
                             logEC50[[ri]] <- NA
                             stderrlogEC50[[ri]] <- NA
+                            a[[ri]] <- NA
                             b[[ri]] <- NA
-                            f[[ri]] <- NA
+                            c[[ri]] <- NA
                         } else {
                             stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+                            a[[ri]] <- coef(m)[["logEC50"]]
                             b[[ri]] <- coef(m)[["b"]]
-                            f[[ri]] <- coef(m)[["f"]]
+                            c[[ri]] <- coef(m)[["f"]]
                         }
                     }
                 }
 
-                if (logis &&
+                if (weibull &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope),
+                    m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
                             data=tmp,
-                            start=list(logEC50=startlogEC50[[i]],slope=1)))
+                            start=list(location=-0.16,shape=0.4659)))
                     if (chooseone==FALSE || fit==FALSE) {
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
@@ -140,29 +139,67 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                             sigma[[ri]] <- s$sigma
                             rsubstance[[ri]] <- i
                             rn[[ri]] <- n
+                            runit[[ri]] <- unit
                             rlld[[ri]] <- log10(lowestdose)
                             rlhd[[ri]] <- log10(highestdose)
-                            mtype[[ri]] <- "logis"
+                            mtype[[ri]] <- "weibull"
+                            a[[ri]] <- coef(m)[["location"]]
+                            b[[ri]] <- coef(m)[["shape"]]
+                            sqrdev <- function(logdose) {
+                                (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 
+                            }
+                            logEC50[[ri]] <- nlm(sqrdev,startlogEC50[[i]])$estimate
+                            c[[ri]] <- NA
+                            if (logEC50[[ri]] > rlhd[[ri]]) {
+                                logEC50[[ri]] <- NA
+                                stderrlogEC50[[ri]] <- NA
+                                a[[ri]] <- NA
+                                b[[ri]] <- NA
+                            } else {
+                                stderrlogEC50[[ri]] <- NA
+                            }
+                        }
+                    }
+                }
+
+                if (logit &&
+                    length(subset(allWrong,allWrong == i))==0) {
+                    m <- try(nls(response ~ plogis(-log10(dose),-logEC50,scale),
+                            data=tmp,
+                            start=list(logEC50=startlogEC50[[i]],scale=1)))
+                    if (chooseone==FALSE || fit==FALSE) {
+                        if (!inherits(m, "try-error")) {
+                            fit <- TRUE
+                            ri <- ri + 1
+                            s <- summary(m)
+                            sigma[[ri]] <- s$sigma
+                            rsubstance[[ri]] <- i
+                            rn[[ri]] <- n
+                            runit[[ri]] <- unit
+                            rlld[[ri]] <- log10(lowestdose)
+                            rlhd[[ri]] <- log10(highestdose)
+                            mtype[[ri]] <- "logit"
                             logEC50[[ri]] <- coef(m)[["logEC50"]]
-                            b[[ri]] <- NA
-                            f[[ri]] <- NA
+                            c[[ri]] <- NA
                             if (logEC50[[ri]] > rlhd[[ri]]) {
                                 logEC50[[ri]] <- NA
-                                slope[[ri]] <- NA
                                 stderrlogEC50[[ri]] <- NA
+                                a[[ri]] <- NA
+                                b[[ri]] <- NA
                             } else {
-                                slope[[ri]] <- coef(m)[["slope"]]
                                 stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+                                a[[ri]] <- coef(m)[["logEC50"]]
+                                b[[ri]] <- coef(m)[["scale"]]
                             }
                         }
                     }
                 }
 
-                if (lognorm &&
+                if (probit &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope),
+                    m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,scale),
                                 data=tmp,
-                                start=list(logEC50=startlogEC50[[i]],slope=1)))
+                                start=list(logEC50=startlogEC50[[i]],scale=1)))
                     if (chooseone==FALSE || fit==FALSE) {
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
@@ -174,17 +211,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                             runit[[ri]] <- unit
                             rlld[[ri]] <- log10(lowestdose)
                             rlhd[[ri]] <- log10(highestdose)
-                            mtype[[ri]] <- "lognorm"
+                            mtype[[ri]] <- "probit"
                             logEC50[[ri]] <- coef(m)[["logEC50"]]
-                            b[[ri]] <- NA
-                            f[[ri]] <- NA
+                            c[[ri]] <- NA
                             if (logEC50[[ri]] > rlhd[[ri]]) {
                                 logEC50[[ri]] <- NA
-                                slope[[ri]] <- NA
                                 stderrlogEC50[[ri]] <- NA
+                                a[[ri]] <- NA
+                                b[[ri]] <- NA
                             } else {
-                                slope[[ri]] <- coef(m)[["slope"]]
                                 stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+                                a[[ri]] <- coef(m)[["logEC50"]]
+                                b[[ri]] <- coef(m)[["scale"]]
                             }
                         }
                     }
@@ -215,19 +253,15 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
             sigma[[ri]] <- NA
             logEC50[[ri]] <- NA
             stderrlogEC50[[ri]] <- NA
-            slope[[ri]] <- NA
+            a[[ri]] <- NA
             b[[ri]] <- NA
-            f[[ri]] <- NA
+            c[[ri]] <- NA
         }
     }
-    results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma)
-    names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma")
-    if (lognorm || logis) {
-        results$slope <- slope
-    }
-    if (linearlogis) {
-        results$b <- b
-        results$f <- f
+    results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma, a, b)
+    names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma","a","b")
+    if (linlogit) {
+        results$c <- c
     }
     return(results)
 }
@@ -257,7 +291,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
     } else {
         lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2
         lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2
-        if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) {
+        if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) {
             hr <- 1.8 
         } else {
             hr <- 1.0
@@ -372,16 +406,21 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
                     {
                         logEC50 <- fits[j,"logEC50"]
                         mtype <- as.character(fits[j, "mtype"])
-                        if (mtype == "lognorm") {
-                            slope <- fits[j,"slope"]
-                            plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+                        if (mtype == "probit") {
+                            scale <- fits[j,"b"]
+                            plot(function(x) pnorm(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
+                        }
+                        if (mtype == "logit") {
+                            scale <- fits[j,"b"]
+                            plot(function(x) plogis(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
                         }
-                        if (mtype == "logis") {
-                            slope <- fits[j,"slope"]
-                            plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+                        if (mtype == "weibull") {
+                            location <- fits[j,"a"]
+                            shape <- fits[j,"b"]
+                            plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color)
                         }
-                        if (mtype == "linearlogis") {
-                            plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]),
+                        if (mtype == "linlogit") {
+                            plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logEC50"],fits[j,"b"]),
                                 lld - 0.5, lhd + 2,
                                 add=TRUE,col=color)
                         }