aboutsummaryrefslogtreecommitdiff
path: root/R/drfit.R
diff options
context:
space:
mode:
authorranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2014-02-03 22:04:37 +0000
committerranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2014-02-03 22:04:37 +0000
commit184aacf1ad5a28b2428633cd1966d6fb881eb3b0 (patch)
tree8317bf1364476fb1681c0afdd090fc2a00265c0f /R/drfit.R
parent88937b957b30b9e197858ebc3a3219079fd01f96 (diff)
- 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
Diffstat (limited to 'R/drfit.R')
-rw-r--r--R/drfit.R64
1 files changed, 60 insertions, 4 deletions
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)
}

Contact - Imprint