diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/checkexperiment.R | 2 | ||||
| -rw-r--r-- | R/checksubstance.R | 1 | ||||
| -rw-r--r-- | R/drdata.R | 1 | ||||
| -rw-r--r-- | R/drfit.R | 64 | ||||
| -rw-r--r-- | R/drplot.R | 1 |
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]) @@ -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") { @@ -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) } @@ -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, |