Added the possibility to simultaneously fit dose-response curves with different

units in drfit. Added cat output from drfit to keep track of which substance is currently
being worked on.



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

1 files changed, 19 insertions, 3 deletions

diff --git a/R/drfit.R b/R/drfit.R
index ff19cca..59deb6a 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -41,15 +41,16 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
         linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA,
         b0 = 2, f0 = 0)
 {
-    if(!is.null(data$ok)) data <- subset(data,ok!="no fit")
+    if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data where ok 
+                                                            # was set to "no fit"
     substances <- levels(data$substance)
-    unit <- levels(as.factor(data$unit))
 
     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
     rn <- vector()                        # number of dose-response curves 
+    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
@@ -63,6 +64,17 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
     for (i in substances) {
         tmp <- splitted[[i]]
         fit <- FALSE
+        if (length(tmp) != 0) {
+            unit <- levels(as.factor(as.vector(tmp$unit)))
+            cat("\n",i,": Fitting data...\n",sep="")
+        } else {
+            unit <- ""
+            cat("\n",i,": No data\n",sep="")
+        }
+        if (length(unit) > 1) {
+            cat("More than one unit for substance ",i,", halting\n\n",sep="")
+            break
+        }
         n <- round(length(tmp$response)/9)
         if (length(tmp$response) == 0) {
             nodata = TRUE
@@ -96,6 +108,7 @@ 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]] <- "linearlogis"
@@ -158,6 +171,7 @@ 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]] <- "lognorm"
@@ -187,9 +201,11 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
             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 {
@@ -204,7 +220,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
             f[[ri]] <- NA
         }
     }
-    results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, unit, sigma)
+    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