diff options
author | Johannes Ranke <jranke@uni-bremen.de> | 2017-03-24 15:28:44 +0100 |
---|---|---|
committer | Johannes Ranke <jranke@uni-bremen.de> | 2017-03-24 15:28:44 +0100 |
commit | baf8e4341b58c09532b4ca37ebf80c9068934db4 (patch) | |
tree | a6abd8cbdfa849a672d2abcfde5af4bfc180bcf6 /R/drcfit.R | |
parent | 22c11d5dbc118881f9cfd7c9c2ac9711295a0f16 (diff) |
Import version 0.6.7 from CRAN
I had deleted the git commits between 0.6.4 and 0.6.7 while removing
drfit from github, while the latest changes were not pushed to jrwb.de
yet.
Diffstat (limited to 'R/drcfit.R')
-rw-r--r-- | R/drcfit.R | 68 |
1 files changed, 37 insertions, 31 deletions
@@ -14,7 +14,7 @@ drcfit <- function(data, chooseone=TRUE, # model result was appended rsubstance <- array() # the substance names in the results rndl <- vector() # number of dose levels - rn <- vector() # mean number of replicates + rn <- vector() # mean number of replicates # in each dose level runit <- vector() # vector of units for each result row rlhd <- rlld <- vector() # highest and lowest doses tested @@ -27,6 +27,16 @@ drcfit <- function(data, chooseone=TRUE, models <- list() # a list containing the dose-response models splitted <- split(data,data$substance) + + # The indexing of the results of the ED50 function changed with drc 3.0.1 + if (packageVersion("drc") > 3) { + ED50_row_index = "e:1:50" + EDx_row_index_prefix = "e:1" + } else { + ED50_row_index = "1:50" + EDx_row_index_prefix = "1" + } + for (i in substances) { tmp <- splitted[[i]] fit <- FALSE @@ -69,8 +79,7 @@ drcfit <- function(data, chooseone=TRUE, active <- FALSE if (linlogit) { - m <- try(drm(response ~ dose, data = tmp, fct = BC.4(fixed = c(NA, 1, NA, NA))), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = BC.4(fixed = c(NA, 1, NA, NA)))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -91,26 +100,25 @@ drcfit <- function(data, chooseone=TRUE, a[[ri]] <- coef(m)[[2]] b[[ri]] <- coef(m)[[1]] c[[ri]] <- coef(m)[[3]] - ED50 <- try(ED(m, 50, interval = "delta", + ED50 <- try(ED(m, 50, interval = "delta", lower = lowestdose / 10, upper = highestdose * 10, - display = FALSE), - silent = TRUE) + display = FALSE)) if (!inherits(ED50, "try-error")) { - logED50[[ri]] <- log10(ED50["1:50", "Estimate"]) - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50[[ri]] <- log10(ED50[ED50_row_index, "Estimate"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" } - } + } } } } if (probit) { - m <- try(drm(response ~ dose, data = tmp, - fct = LN.2()), silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, + fct = LN.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -136,16 +144,15 @@ drcfit <- function(data, chooseone=TRUE, } else { mtype[[ri]] <- "probit" ED50 <- ED(m, 50, interval = "delta", display = FALSE) - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) } } } } if (logit) { - m <- try(drm(response ~ dose, data = tmp, fct = LL.2()), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = LL.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -173,8 +180,8 @@ drcfit <- function(data, chooseone=TRUE, } else { mtype[[ri]] <- "logit" ED50 <- ED(m, 50, interval = "delta", display = FALSE) - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) } } } @@ -182,8 +189,7 @@ drcfit <- function(data, chooseone=TRUE, } if (weibull) { - m <- try(drm(response ~ dose, data = tmp, fct = W1.2()), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = W1.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -199,7 +205,7 @@ drcfit <- function(data, chooseone=TRUE, rlhd[[ri]] <- log10(highestdose) c[[ri]] <- NA ED50 <- ED(m, 50, interval = "delta", display = FALSE) - logED50[[ri]] <- log10(ED50["1:50", "Estimate"]) + logED50[[ri]] <- log10(ED50[ED50_row_index, "Estimate"]) if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA @@ -209,8 +215,8 @@ drcfit <- function(data, chooseone=TRUE, b[[ri]] <- NA } else { mtype[[ri]] <- "weibull" - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) a[[ri]] <- logED50[[ri]] b[[ri]] <- coef(m)[[1]] } @@ -256,7 +262,7 @@ drcfit <- function(data, chooseone=TRUE, } } - results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, + 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 = "") @@ -279,14 +285,14 @@ drcfit <- function(data, chooseone=TRUE, mtype <- as.character(results[row.i, "mtype"]) if (mtype %in% c("probit", "logit", "weibull", "linlogit")) { for (EDi in EDx) { - EDx.drc = try(ED(m, EDi, interval = "delta", display = FALSE, level = level), - silent = TRUE) + EDx.drc = try(ED(m, EDi, interval = "delta", display = FALSE, level = level)) if (!inherits(EDx.drc, "try-error")) { - results[row.i, paste0("EDx", EDi)] <- EDx.drc[paste0("1:", EDi), "Estimate"] - results[row.i, paste0("EDx", EDi, " ", lower_level_percent)] <- EDx.drc[paste0("1:", EDi), - "Lower"] - results[row.i, paste0("EDx", EDi, " ", upper_level_percent)] <- EDx.drc[paste0("1:", EDi), - "Upper"] + results[row.i, paste0("EDx", EDi)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Estimate"] + results[row.i, paste0("EDx", EDi, " ", lower_level_percent)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Lower"] + results[row.i, paste0("EDx", EDi, " ", upper_level_percent)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Upper"] } } } |