- The drplot function now accepts custom xlim and ylim values

- Confidence intervals for the EC50 can now (again) be generated by 
  drfit by giving the argument conf=TRUE
- Update of the INDEX



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

6 files changed, 63 insertions, 23 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 4fd580b..fb4798a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: drfit
-Version: 0.04-60
-Date: 2006-04-01
+Version: 0.04-62
+Date: 2006-04-03
 Title: Dose-response data evaluation
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
diff --git a/INDEX b/INDEX
index c02adfe..12b1206 100644
--- a/INDEX
+++ b/INDEX
@@ -4,6 +4,9 @@ drdata                  Get dose-response data
 drfit                   Fit dose-response models
 drfit-package           Dose-response data evaluation
 drplot                  Plot dose-response models
+IM1xIPC81               Dose-Response data for
+                        1-methyl-3-alkylimidazolium tetrafluoroborates
+                        in IPC-81 cells
 IM1xVibrio              Dose-Response data for
                         1-methyl-3-alkylimidazolium tetrafluoroborates
                         in V. fischeri
diff --git a/R/drfit.R b/R/drfit.R
index ff03d90..1bbffa8 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,6 +1,7 @@
 drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
         probit = TRUE, logit = FALSE, weibull = FALSE,
-        linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
+        linlogit = FALSE, conf = FALSE,
+        linlogitWrong = NA, allWrong = NA,
         s0 = 0.5, b0 = 2, f0 = 0)
 {
     if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data with 
@@ -20,6 +21,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     sigma <- array()                      # the standard deviation of the residuals
     logED50 <- vector()
     stderrlogED50 <- vector()
+    conflogED50 <- vector()
     a <- b <- c <- vector()
 
     splitted <- split(data,data$substance)
@@ -47,7 +49,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
             n <- ndl <- 0
         } else {
             ndl <- length(levels(factor(tmp$dose)))
-            n <- round(length(tmp$response)/ndl)
+            n <- length(tmp$response)
             if (is.na(startlogED50[i])){
                 w <- 1/abs(tmp$response - 0.3)
                 startlogED50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
@@ -86,12 +88,14 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                 mtype[[ri]] <- "no fit"
                                 logED50[[ri]] <- NA
                                 stderrlogED50[[ri]] <- NA
+                                conflogED50[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
                                 c[[ri]] <- NA
                             } else {
                                 mtype[[ri]] <- "linlogit"
                                 stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 3)
                                 a[[ri]] <- coef(m)[["logED50"]]
                                 b[[ri]] <- coef(m)[["b"]]
                                 c[[ri]] <- coef(m)[["f"]]
@@ -122,11 +126,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
                                     stderrlogED50[[ri]] <- NA
+                                    conflogED50[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "probit"
                                     stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                     a[[ri]] <- coef(m)[["logED50"]]
                                     b[[ri]] <- coef(m)[["scale"]]
                                 }
@@ -158,11 +164,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
                                     stderrlogED50[[ri]] <- NA
+                                    conflogED50[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "logit"
                                     stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                 }
                             }
                         }
@@ -196,11 +204,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
                                     stderrlogED50[[ri]] <- NA
+                                    conflogED50[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "weibull"
                                     stderrlogED50[[ri]] <- NA
+                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                 }
                             }
                         }
@@ -238,13 +248,21 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
             sigma[[ri]] <- NA
             logED50[[ri]] <- NA
             stderrlogED50[[ri]] <- NA
+            conflogED50[[ri]] <- NA
             a[[ri]] <- NA
             b[[ri]] <- NA
             c[[ri]] <- NA
         }
     }
-    results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b)
-    names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b")
+    if (conf)
+    {
+        results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, conflogED50, runit, sigma, a, b)
+        names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","conf","unit","sigma","a","b")
+    } else {
+        results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b)
+        names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b")
+    } 
+
     if (linlogit) {
         results$c <- c
     }
diff --git a/R/drplot.R b/R/drplot.R
index baec69d..8faaa0f 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -1,6 +1,7 @@
 drplot <- function(drresults, data, 
         dtype = "std", alpha = 0.95, ctype = "none",
         path = "./", fileprefix = "drplot", overlay = FALSE,
+        xlim = c("auto","auto"), ylim = c("auto","auto"),
         postscript = FALSE, pdf = FALSE, png = FALSE, 
         bw = TRUE,
         pointsize = 12,
@@ -14,14 +15,15 @@ drplot <- function(drresults, data,
         unit <- "different units"
     }
 
-    # Get the plot limits on the x-axis (log of the dose) and y axis
+    # Determine the plot limits on the x-axis and y axis
     if(is.data.frame(data)) {
+        # Get rid of pseudo substance names of controls
         nonzerodata <- subset(data,dose!=0)
-        nonzerodata$substance <- factor(nonzerodata$substance) # Get rid of pseudo substance names of controls
+        nonzerodata$substance <- factor(nonzerodata$substance) 
         zerodata <- subset(data,dose==0)
         nc <- length(zerodata$dose)     # Number of control points
         if (nc > 0) {
-            sdc <- sd(zerodata$response)    # Standard deviation of control responses
+            sdc <- sd(zerodata$response)
             controlconf <- sdc * qt((1 + alpha)/2, nc - 1) / sqrt(nc)
             cat("There are ",nc,"data points with dose 0 (control values)\n")
             cat("with a standard deviation of",sdc,"\n")
@@ -52,6 +54,12 @@ drplot <- function(drresults, data,
             hr <- 1.0
         }
     }
+    if (xlim[1] == "auto") xlim[1] <- lld - 0.5
+    if (xlim[2] == "auto") xlim[2] <- lhd + 1
+    if (ylim[1] == "auto") ylim[1] <- -0.1
+    if (ylim[2] == "auto") ylim[2] <- hr + 0.2
+    xlim <- as.numeric(xlim)
+    ylim <- as.numeric(ylim)
 
     # Prepare overlay plot if requested
     if (overlay)
@@ -79,10 +87,10 @@ drplot <- function(drresults, data,
         }
             
         plot(0,type="n",
-            xlim=c(lld - 0.5, lhd + 1),
-            ylim= c(-0.1, hr + 0.2),
-            xlab=paste("Decadic Logarithm of the dose in ", unit),    
-            ylab="Normalized response")
+            xlim = xlim,
+            ylim = ylim,
+            xlab = paste("Decadic Logarithm of the dose in ", unit),    
+            ylab = "Normalized response")
     }
 
     # Plot the data either as raw data or as error bars
@@ -123,10 +131,10 @@ drplot <- function(drresults, data,
                     }
                         
                     plot(0,type="n",
-                        xlim=c(lld - 0.5, lhd + 2),
-                        ylim= c(-0.1, hr + 0.2),
-                        xlab=paste("Decadic Logarithm of the dose in ", unit),    
-                        ylab="Normalized response")
+                        xlim = xlim,
+                        ylim = ylim,
+                        xlab = paste("Decadic Logarithm of the dose in ", unit),    
+                        ylab = "Normalized response")
                 }
                 if (!overlay) legend(lpos, i,lty = 1, col = color, inset=0.05)
                 tmp$dosefactor <- factor(tmp$dose)  # necessary because the old
diff --git a/man/drfit.Rd b/man/drfit.Rd
index 6e9ef47..9d87b4c 100644
--- a/man/drfit.Rd
+++ b/man/drfit.Rd
@@ -7,8 +7,8 @@
 }
 \usage{
   drfit(data, startlogED50 = NA, chooseone = TRUE, probit = TRUE, logit = FALSE,
-    weibull = FALSE, linlogit = FALSE, linlogitWrong = NA, allWrong = NA, 
-    s0 = 0.5, b0 = 2, f0 = 0)
+	weibull = FALSE, linlogit = FALSE, conf = FALSE, linlogitWrong = NA,
+	allWrong = NA, s0 = 0.5, b0 = 2, f0 = 0)
 }
 \arguments{
   \item{data}{
@@ -44,6 +44,10 @@
     A boolean defining if the linear-logistic function
     \code{\link{linlogitf}} as defined by van Ewijk and Hoekstra 1993 is
     fitted to the data. Default is FALSE.}
+  \item{conf}{
+	A boolean defining if a confidence interval of the log ED50 is to be
+	listed for alpha = 5 \% (two-sided). If FALSE (default), the standard
+	deviation is listed in the output of drfit.}
   \item{linlogitWrong}{
     An optional vector containing the names of the substances for which the
     linlogit function produces a wrong fit.}
@@ -73,8 +77,8 @@
 	  if the fitted ED50 is not higher than the highest dose. 
 
 	  \code{ndl} is the number of dose levels in the raw data, \code{n} is the
-	  rounded mean of the number of replicates at each dose level in the raw
-	  data, \code{lld} is the decadic logarithm of the lowest dose and
+	  total number of data points in the raw data used for the fit
+	  \code{lld} is the decadic logarithm of the lowest dose and
 	  \code{lhd} is the decadic logarithm of the highest dose.  For the
 	  \dQuote{linlogit}, \dQuote{logit} and \dQuote{probit} models, the
 	  parameter \code{a} that is reported coincides with the logED50, i.e the
diff --git a/man/drplot.Rd b/man/drplot.Rd
index 4c9e298..0249476 100644
--- a/man/drplot.Rd
+++ b/man/drplot.Rd
@@ -6,8 +6,9 @@
   either combined or separately, for one or more substances.
 }
 \usage{
-  drplot(drresults, data, dtype, alpha, ctype, path, fileprefix, overlay,
-    postscript, pdf, png, bw, pointsize, colors, devoff, lpos)
+  drplot(drresults, data, dtype, alpha, ctype, path,
+     fileprefix, overlay, xlim, ylim, postscript, pdf, png, bw,
+     pointsize, colors, devoff, lpos)
 }
 \arguments{
   \item{drresults}{
@@ -51,6 +52,12 @@
     (postscript=FALSE) only works correctly for data plots. Dose-response models
     will all be put into the last graph in this case.
     }
+  \item{xlim}{
+    The plot limits (min,max) on the dose axis.
+  }
+  \item{ylim}{
+    The plot limits (min,max) on the response axis.
+  }
   \item{postscript}{
     If TRUE, (a) postscript graph(s) will be created. Otherwise, and if 
     the pdf and png arguments are also FALSE, graphics will be