diff options
Diffstat (limited to 'R/drfit.R')
| -rw-r--r-- | R/drfit.R | 44 |
1 files changed, 13 insertions, 31 deletions
@@ -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: |