Improvements on the EDx calculations in drcfit. However, in cases that there is no

hormesis in the data, drfit often fails to fit the linlogit model, and drcfit
often fails to find the EDx values.



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

5 files changed, 29 insertions, 20 deletions

diff --git a/ChangeLog b/ChangeLog
index 79ee431..35a8a35 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,4 @@
-2014-02-25  Johannes Ranke
+2014-03-18  Johannes Ranke
 
 	* DESCRIPTION, NAMESPACE: New verson 0.6.1, with the drc package as
 	dependency. Import needed functionality to avoid masking standard functions.
@@ -7,7 +7,8 @@
 	drfit(), but internally uses the methods supplied by the drc package
 
 	* R/{drfit,drcfit}.R: Return the list of fitted models in an attribute of
-	the resulting dataframe
+	the resulting dataframe. Deal with the case of only NA values in the unit 
+	column of the data.
 
 	* R/drfit.R: Use the fitted nls() models directly to calculate EDx values
 
diff --git a/DESCRIPTION b/DESCRIPTION
index f706688..4bf9df4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: drfit
 Version: 0.6.1
-Date: 2014-02-25
+Date: 2014-03-18
 Title: Dose-response data evaluation
 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre"), 
                     email = "jranke@uni-bremen.de"))
diff --git a/R/drcfit.R b/R/drcfit.R
index ff9e6bb..8c52f6d 100644
--- a/R/drcfit.R
+++ b/R/drcfit.R
@@ -37,6 +37,9 @@ drcfit <- function(data, chooseone=TRUE,
             unit <- ""
             message("\n",i,": No data\n")
         }
+        if (length(unit) == 0) {
+            unit <- NA
+        }
         if (length(unit) > 1) {
             message("More than one unit for substance ",i,", halting\n\n")
             break
@@ -81,23 +84,24 @@ drcfit <- function(data, chooseone=TRUE,
                                 runit[[ri]] <- unit
                                 rlld[[ri]] <- log10(lowestdose)
                                 rlhd[[ri]] <- log10(highestdose)
-                                ED50 <- try(ED(m, 50, interval = "delta", display = FALSE))
+                                logED50[[ri]] <- NA
+                                logED50low[[ri]] <- NA
+                                logED50high[[ri]] <- NA
+                                a[[ri]] <- coef(m)[[2]]
+                                b[[ri]] <- coef(m)[[1]]
+                                c[[ri]] <- coef(m)[[3]]
+                                ED50 <- try(ED(m, 50, interval = "delta", 
+                                               lower = lowestdose / 10,
+                                               upper = highestdose * 10,
+                                               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]]
+                                } 
                             }
                         }
                     }
@@ -119,21 +123,19 @@ drcfit <- function(data, chooseone=TRUE,
                                 rlld[[ri]] <- log10(lowestdose)
                                 rlhd[[ri]] <- log10(highestdose)
                                 logED50[[ri]] <- log10(coef(m)[[2]])
+                                a[[ri]] <- coef(m)[[2]]
+                                b[[ri]] <- coef(m)[[1]]
                                 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]]
                                 }
                             }
                         }
@@ -155,6 +157,8 @@ drcfit <- function(data, chooseone=TRUE,
                                 rlld[[ri]] <- log10(lowestdose)
                                 rlhd[[ri]] <- log10(highestdose)
                                 logED50[[ri]] <- log10(coef(m)[[2]])
+                                a[[ri]] <- coef(m)[[2]]
+                                b[[ri]] <- coef(m)[[1]]
                                 c[[ri]] <- NA
                                 if (logED50[[ri]] > rlhd[[ri]]) {
                                     mtype[[ri]] <- "no fit"
@@ -168,8 +172,6 @@ drcfit <- function(data, chooseone=TRUE,
                                     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]]
                                 }
                             }
                         }
diff --git a/R/drfit.R b/R/drfit.R
index ad7d16c..7fa53bf 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -41,6 +41,9 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
             unit <- ""
             message("\n",i,": No data\n")
         }
+        if (length(unit) == 0) {
+            unit <- NA
+        }
         if (length(unit) > 1) {
             message("More than one unit for substance ",i,", halting\n\n")
             break
diff --git a/R/drplot.R b/R/drplot.R
index ded7aca..11c4d42 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -204,7 +204,10 @@ drplot <- function(drresults, data,
                         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")) 
+                        if (drresults[[fit.row, "mtype"]] %in% c("probit",
+                                                                 "logit",
+                                                                 "weibull",
+                                                                 "linlogit")) 
                         {
                             m <- attr(drresults, "models")[[as.numeric(fit.row)]]
                             of <- function(x) {