- Some updates to the packaging

- Add the possibility to calculate EDx values
- see ChangeLog for a full description



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

5 files changed, 62 insertions, 7 deletions

diff --git a/R/checkexperiment.R b/R/checkexperiment.R
index 58fecad..2bf0ce7 100644
--- a/R/checkexperiment.R
+++ b/R/checkexperiment.R
@@ -1,3 +1,4 @@
+if(getRversion() >= '2.15.1') utils::globalVariables(c("type", "conc"))
 checkexperiment <- function(id, db = "ecotox", endpoint = "%")
 {
     databases <- data.frame(
@@ -8,7 +9,6 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%")
 
     if (!(db %in% rownames(databases))) stop("Database is not supported")
 
-    library(RODBC) 
     channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
 
 
diff --git a/R/checksubstance.R b/R/checksubstance.R
index 96be999..b50c3da 100644
--- a/R/checksubstance.R
+++ b/R/checksubstance.R
@@ -11,7 +11,6 @@ checksubstance <- function(substance, db = "cytotox", experimentator = "%",
 
     if (!(db %in% rownames(databases))) stop("Database is not supported")
 
-    library(RODBC) 
     channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
 
     responsename = as.character(databases[db,1])
diff --git a/R/drdata.R b/R/drdata.R
index 9da7d96..15c61ac 100644
--- a/R/drdata.R
+++ b/R/drdata.R
@@ -3,7 +3,6 @@ drdata <- function(substances, experimentator = "%", db = "cytotox",
     organism="Vibrio fischeri",endpoint="Luminescence",whereClause="1",
     ok="'ok','no fit'")
 {
-    library(RODBC) 
     channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
     slist <- paste(substances,collapse="','")
     if (db == "cytotox") {
diff --git a/R/drfit.R b/R/drfit.R
index 462bc43..fc73632 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,11 +1,13 @@
+if(getRversion() >= '2.15.1') utils::globalVariables(c("ok", "dose"))
 drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
         probit = TRUE, logit = FALSE, weibull = FALSE,
         linlogit = FALSE, level = 0.95,
         linlogitWrong = NA, allWrong = NA,
         ps0 = 1, ls0 = 0.5, ws0 = 0.5,
-        b0 = 2, f0 = 0)
+        b0 = 2, f0 = 0,
+        showED50 = FALSE,
+        EDx = NULL, EDx.tolerance = 1e-4)
 {
-    require(MASS)
     if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data
                                                             # with ok set to
                                                             # "no fit"
@@ -282,14 +284,68 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     }
     results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, 
         logED50, logED50low, logED50high, runit, sigma, a, b)
+    lower_level_percent = paste(100 * (1 - level)/2, "%", sep = "")
+    upper_level_percent = paste(100 * (1 + level)/2, "%", sep = "")
     names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50",
-        paste(100*(1-level)/2,"%",sep=""),
-        paste(100*(1+level)/2,"%",sep=""),
+        lower_level_percent, upper_level_percent,
         "unit","sigma","a","b")
 
     if (linlogit) {
         results$c <- c
     }
+    if (showED50) {
+        results[c("ED50", paste("ED50", c(lower_level_percent, upper_level_percent)))] <-
+          10^results[7:9]
+    }
+    if (!is.null(EDx)) {
+        for (row.i in 1:ri) {
+            logED50 <- results[row.i, "logED50"]
+            mtype <- as.character(results[row.i, "mtype"])
+            if (mtype == "probit") {
+                scale <- results[row.i, "b"]
+                drfunction <- function(x) pnorm(-x, -logED50, scale)
+            }
+            if (mtype == "logit") {
+                scale <- results[row.i, "b"]
+                drfunction <- function(x) plogis(-x, -logED50, scale)
+            }
+            if (mtype == "weibull") {
+                location <- results[row.i, "a"]
+                shape <- results[row.i, "b"]
+                drfunction <- function(x) pweibull(-x + location, shape)
+            }
+            if (mtype == "linlogit") {
+                drfunction <- function(x) {
+                       r <- linlogitf(10^x, 1, 
+                                      results[row.i, "c"], 
+                                      results[row.i, "logED50"],
+                                      results[row.i, "b"])
+                       # Do not return response values above 1 to avoid trapping
+                       # the optimization algorithm in a local minimum later on
+                       if (r < 1) return(r)
+                       else return(2 - 0.001 * x) # Start with 2 and decrease slowly to 
+                                                  # guide to the interesting part of the curve
+                }
+            }
+            if (mtype %in% c("probit", "logit", "weibull", "linlogit")) {
+                for (ED in EDx) {
+                    of <- function(x) abs(drfunction(x) - (1 - (ED/100)))
+                    
+                    # Search over interval starting an order of magnitude below 
+                    # the lowest dose up to one order of magnitude above the
+                    # highest dose
+                    o = optimize(of, 
+                                 results[row.i, c("lld", "lhd")] + c(-1, 1))
+                    # Only keep results within the tolerance
+                    if ((o$objective) < EDx.tolerance) {
+                        logdose.ED = o$minimum
+                        results[row.i, paste0("EDx", ED)] <- 10^logdose.ED 
+                    }
+                }
+            }
+        }
+    }
+
     rownames(results) <- 1:ri
     return(results)
 }
diff --git a/R/drplot.R b/R/drplot.R
index 55b9da4..149337c 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -1,3 +1,4 @@
+if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance"))
 drplot <- function(drresults, data, 
         dtype = "std", alpha = 0.95, ctype = "none",
         path = "./", fileprefix = "drplot", overlay = FALSE,