aboutsummaryrefslogtreecommitdiff
path: root/R/drcfit.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/drcfit.R')
-rw-r--r--R/drcfit.R68
1 files changed, 37 insertions, 31 deletions
diff --git a/R/drcfit.R b/R/drcfit.R
index 0d27c1c..008c53c 100644
--- a/R/drcfit.R
+++ b/R/drcfit.R
@@ -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"]
}
}
}

Contact - Imprint