From 35362cfc7c02f12fd9689f68a772a804d895535a Mon Sep 17 00:00:00 2001 From: ranke Date: Tue, 25 Feb 2014 17:30:01 +0000 Subject: Feature-complete version of the new drfit package that can make use of the drc package in order to obtain confidence intervals for EDx values. See ChangeLog for details. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@98 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drfit.R | 44 +++++++++++++------------------------------- 1 file changed, 13 insertions(+), 31 deletions(-) (limited to 'R/drfit.R') diff --git a/R/drfit.R b/R/drfit.R index fc73632..ad7d16c 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -28,6 +28,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, logED50low <- logED50high <- vector() a <- b <- c <- vector() + models <- list() # a list containing the dose-response models + splitted <- split(data,data$substance) for (i in substances) { tmp <- splitted[[i]] @@ -79,6 +81,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 + models[[ri]] <- m s <- summary(m) sigma[[ri]] <- s$sigma rsubstance[[ri]] <- i @@ -122,6 +125,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 + models[[ri]] <- m s <- summary(m) sigma[[ri]] <- s$sigma rsubstance[[ri]] <- i @@ -165,6 +169,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 + models[[ri]] <- m s <- summary(m) sigma[[ri]] <- s$sigma rsubstance[[ri]] <- i @@ -206,6 +211,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { ri <- ri + 1 + models[[ri]] <- m a[[ri]] <- coef(m)[["location"]] b[[ri]] <- coef(m)[["shape"]] sqrdev <- function(logdose) { @@ -299,38 +305,12 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, } 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")) { + if (mtype[[row.i]] %in% c("probit", "logit", "weibull", "linlogit")) { for (ED in EDx) { - of <- function(x) abs(drfunction(x) - (1 - (ED/100))) - + of <- function(x) { + abs(predict(models[[row.i]], data.frame(dose = 10^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 @@ -347,5 +327,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, } rownames(results) <- 1:ri + attr(results, "models") <- models return(results) } +# vim: set ts=4 sw=4 expandtab: -- cgit v1.2.1