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

5 files changed, 149 insertions, 85 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 3eee40b..6b1b85b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,10 +1,10 @@
 Package: drfit
-Version: 0.03-28
+Version: 0.03-29
 Date: 2005-07-16
 Title: Dose-response data evaluation
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
-Depends: R (>=2.1.0),stats,RODBC
+Depends: R (>= 2.1.0),stats,RODBC
 Description: drfit provides basic functions for fitting dose-response curves to
 	dose-response data, calculating some (eco)toxicological parameters and
 	plotting the results. Functions that are fitted are the cumulative densitiy
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)
                         }
diff --git a/man/drfit.Rd b/man/drfit.Rd
index 2ec4a73..0d9ee6f 100644
--- a/man/drfit.Rd
+++ b/man/drfit.Rd
@@ -6,58 +6,84 @@
   biometric results for (eco)toxicity evaluation
 }
 \usage{
-  drfit(data, startlogEC50 = NA, chooseone = TRUE, lognorm = TRUE, logis = FALSE, 
-    linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA, b0 = 2, f0 = 0)
+  drfit(data, startlogEC50 = NA, chooseone = TRUE, probit = TRUE, logit = FALSE,
+    weibull = FALSE, linlogit = FALSE, linlogitWrong = NA, allWrong = NA, 
+    b0 = 2, f0 = 0)
 }
 \arguments{
   \item{data}{
-    A data frame as returned from \code{\link{drdata}}.  The data frame has to
+    A data frame containing dose-response data. The data frame has to
     contain at least a factor called "substance", a vector called "unit"
     containing the unit used for the dose, a column "response" with the
     response values of the test system normalized between 0 and 1 and a column
-    "dose" with the numeric dose values. For later use of the
-    \code{\link{drplot}} function, a factor called "dosefactor" also has to be
-    present, containing the dose as a factor. If there is a column called "ok"
-    and it is set to "no fit", then the corresponding data point will be excluded
-    from the fitting procedure.
+    "dose" with the numeric dose values. Such a data frame can be easily obtained
+    if a compliant RODBC data source is available for use in conjunction with 
+    the function \code{\link{drdata}}.  
+    
+    If there is a column called "ok" and it is set to "no fit" in a specific
+    line, then the corresponding data point will be excluded from the fitting
+    procedure, although it will be plotted.
     }
   \item{startlogEC50}{
-    Especially for the linearlogis model, a suitable log10 of the EC50 has to be given 
+    Especially for the linlogit model, a suitable log10 of the EC50 has to be given 
     by the user, since it is not correctly estimated for data showing hormesis with
     the default estimation method.}
-  \item{lognorm}{
+  \item{probit}{
     A boolean defining if cumulative density curves of normal distributions
-    are fitted to the data. Default ist TRUE.} 
-  \item{logis}{
-    A boolean defining if cumulative densitiy curves of logistic distributions
-    are fitted to the data. Default is FALSE.} 
-  \item{linearlogis}{
+    \code{\link{pnorm}} are fitted against the decadic logarithm of the dose.
+    Default ist TRUE.} 
+  \item{logit}{
+    A boolean defining if cumulative density curves of logistic distributions
+    \code{\link{plogis}} are fitted to the decadic logarithm of the dose.
+    Default is FALSE.} 
+  \item{weibull}{
+    A boolean defining if the cumulative density curves of weibull distributions
+    (\code{\link{pweibull}} with additionall location parameter and scale=1)
+    are fitted to the decadic logarithm of the dose. Default is FALSE.}
+  \item{linlogit}{
     A boolean defining if the linear-logistic function
-    \code{\link{linearlogisf}} as defined by van Ewijk and Hoekstra 1993 is
+    \code{\link{linlogitf}} as defined by van Ewijk and Hoekstra 1993 is
     fitted to the data. Default is FALSE.}
-  \item{linearlogisWrong}{
+  \item{linlogitWrong}{
     An optional vector containing the names of the substances for which the
-    linearlogis function produces a wrong fit.}
+    linlogit function produces a wrong fit.}
   \item{allWrong}{
     An optional vector containing the names of the substances for which all
-    functions produces a wrong fit.}
+    functions produce a wrong fit.}
   \item{chooseone}{
-    If TRUE (default), the models are tried in the order linearlogis, logis and lognorm,
-    and the first model that produces a valid fit is used. Usually this will be the one 
-    with the lowest residual standard deviation. If FALSE, all models that are set to TRUE
-    and that can be fitted will be reported.}
+    If TRUE (default), the models are tried in the order linlogit, weibull,
+    logit, probit, and the first model that produces a valid fit is used.
+    If FALSE, all models that are set to TRUE and that can be fitted will be
+    reported.}
   \item{b0,f0}{
     If the linearlogistic model is fitted, b0 and f0 give the possibility to
     adapt the starting values for the parameters b and f.}
 }
 \value{
   \item{results}{
-      A data frame containing at least one line for each substance. If the data did not
-      show a mean response < 0.5 at the highest dose level, the modeltype is set to "none".
-      Every successful fit is reported in one line. Parameters of the fitted curves are only
-      reported if the fitted EC50 is not higher than the highest dose.}
-    
-} 
+      A data frame containing at least one line for each substance. If the data
+      did not show a mean response < 0.5 at the highest dose level, the
+      modeltype is set to "none".
+      Every successful fit is reported in one line. Parameters of the fitted
+      curves are only reported if the fitted EC50 is not higher than the
+      highest dose. 
+      \code{n} is the number of dose-response curves in the raw data (repetitions
+      in each point), \code{lld} is the decadic logarithm of the lowest dose and
+      \code{lhd} is the decadic logarithm of the highest dose.
+      For the "linlogit", "logit" and "probit" models, the parameter
+      \code{a} that is reported coincides with the logEC50, i.e the logEC50 is 
+      one of the model parameters that is being fitted, and therefore
+      a standard deviation \code{std} is reported for the logEC50. In the 
+      case of the "weibull" model, \code{a} is a location parameter.
+      Parameter \code{b} in the case of the "linlogit" fit is the variable
+      b from the \code{\link{linlogitf}} function. In the case of "probit" fit
+      it is the standard deviation of the fitted normal distribution, in the case
+      of the "logit" fit it is the \code{scale} parameter in the \code{\link{plogis}}
+      function, and in the "weibull" fit it is the \code{shape} parameter of the
+      fitted \code{\link{pweibull}} function. Only the "linlogit" fit produces a 
+      third parameter \code{c} which is the variable f from the
+      \code{\link{linlogitf}} function.}
+  } 
 \examples{
 data(antifoul)
 r <- drfit(antifoul)
diff --git a/man/drplot.Rd b/man/drplot.Rd
index ecd6895..e5201a2 100644
--- a/man/drplot.Rd
+++ b/man/drplot.Rd
@@ -80,9 +80,8 @@
     
 } 
 \note{
-  Treatment of legends is not really well solved. Be sure all devices are
-  closed (e.g. by calling \code{dev.off()}) before calling \code{drplot}
-  again after a failure.
+  Be sure all devices are closed (e.g. by calling \code{dev.off()}) before
+  calling \code{drplot} again after a failure.
 }
 \examples{
 data(antifoul)
diff --git a/man/linearlogisf.Rd b/man/linlogitf.Rd
index 1ebe8be..2eca8c4 100644..100755
--- a/man/linearlogisf.Rd
+++ b/man/linlogitf.Rd
@@ -1,12 +1,12 @@
-\name{linearlogisf}
-\alias{linearlogisf}
+\name{linlogitf}
+\alias{linlogitf}
 \title{Linear-logistic function}
 \description{
   Helper function describing a special type of dose-response curves, showing a stimulus
   at subtoxic doses.
 }
 \usage{
-  linearlogisf(x,k,f,mu,b)
+  linlogitf(x,k,f,mu,b)
 }
 \arguments{
   \item{x}{