diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/drcfit.R | 292 | ||||
| -rw-r--r-- | R/drfit.R | 44 | ||||
| -rw-r--r-- | R/drplot.R | 57 |
3 files changed, 328 insertions, 65 deletions
diff --git a/R/drcfit.R b/R/drcfit.R new file mode 100644 index 0000000..ff9e6bb --- /dev/null +++ b/R/drcfit.R @@ -0,0 +1,292 @@ +drcfit <- function(data, chooseone=TRUE, + probit = TRUE, logit = FALSE, weibull = FALSE, + linlogit = FALSE, level = 0.95, + showED50 = FALSE, + EDx = NULL) +{ + if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data + # with ok set to + # "no fit" + substances <- levels(data$substance) + + ri <- rix <- 0 # ri is the index over the result rows + # rix is used later to check if any + # model result was appended + rsubstance <- array() # the substance names in the results + rndl <- vector() # number of dose levels + 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 + mtype <- array() # the modeltypes + sigma <- array() # the standard deviation of the residuals + logED50 <- vector() + 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]] + fit <- FALSE + if (length(tmp) != 0) { + unit <- levels(as.factor(as.vector(tmp$unit))) + message("\n",i,": Fitting data...\n") + } else { + unit <- "" + message("\n",i,": No data\n") + } + if (length(unit) > 1) { + message("More than one unit for substance ",i,", halting\n\n") + break + } + if (length(tmp$response) == 0) { + nodata = TRUE + } else { + nodata = FALSE + } + rix <- ri + if (nodata) { + n <- ndl <- 0 + } else { + ndl <- length(levels(factor(tmp$dose))) + n <- length(tmp$response) + highestdose <- max(tmp$dose) + lowestdose <- min(tmp$dose) + lhd <- log10(highestdose) + lld <- log10(lowestdose) + responseathighestdose <- mean(subset(tmp,dose==highestdose)$response) + responseatlowestdose <- mean(subset(tmp,dose==lowestdose)$response) + if (responseathighestdose < 0.5) { + inactive <- FALSE + if (responseatlowestdose < 0.5) { + active <- TRUE + } else { + active <- FALSE + if (linlogit) + { + 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 + ri <- ri + 1 + mtype[[ri]] <- "linlogit" + models[[ri]] <- m + s <- summary(m) + sigma[[ri]] <- s$rseMat[1, "rse"] + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + ED50 <- try(ED(m, 50, interval = "delta", 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"]) + if (logED50[[ri]] > rlhd[[ri]]) { + mtype[[ri]] <- "no fit" + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + } + } else { + mtype[[ri]] <- "no ED50" + } + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] + c[[ri]] <- coef(m)[[3]] + } + } + } + if (probit) + { + m <- try(drm(response ~ dose, data = tmp, + fct = LN.2())) + if (chooseone==FALSE || fit==FALSE) { + if (!inherits(m, "try-error")) { + fit <- TRUE + ri <- ri + 1 + models[[ri]] <- m + s <- summary(m) + sigma[[ri]] <- s$rseMat[1, "rse"] + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + logED50[[ri]] <- log10(coef(m)[[2]]) + c[[ri]] <- NA + if (logED50[[ri]] > rlhd[[ri]]) { + mtype[[ri]] <- "no fit" + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + } 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"]) + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] + } + } + } + } + if (logit) + { + m <- try(drm(response ~ dose, data = tmp, fct = LL.2())) + if (chooseone==FALSE || fit==FALSE) { + if (!inherits(m, "try-error")) { + fit <- TRUE + ri <- ri + 1 + models[[ri]] <- m + s <- summary(m) + sigma[[ri]] <- s$rseMat[1, "rse"] + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + logED50[[ri]] <- log10(coef(m)[[2]]) + c[[ri]] <- NA + if (logED50[[ri]] > rlhd[[ri]]) { + mtype[[ri]] <- "no fit" + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + } 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"]) + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] + } + } + } + + } + if (weibull) + { + m <- try(drm(response ~ dose, data = tmp, fct = W1.2())) + if (chooseone==FALSE || fit==FALSE) { + if (!inherits(m, "try-error")) { + fit <- TRUE + ri <- ri + 1 + models[[ri]] <- m + s <- summary(m) + sigma[[ri]] <- s$rseMat[1, "rse"] + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + c[[ri]] <- NA + ED50 <- ED(m, 50, interval = "delta", display = FALSE) + logED50[[ri]] <- log10(ED50["1:50", "Estimate"]) + if (logED50[[ri]] > rlhd[[ri]]) { + mtype[[ri]] <- "no fit" + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + } else { + mtype[[ri]] <- "weibull" + logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) + logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + a[[ri]] <- logED50[[ri]] + b[[ri]] <- coef(m)[[1]] + } + } + } + + } + } + } else { + inactive <- TRUE + } + } + if (ri == rix) { # if no entry was appended for this substance + ri <- ri + 1 + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + if (nodata) { + rlld[[ri]] <- rlhd[[i]] <- NA + mtype[[ri]] <- "no data" + runit[[ri]] <- NA + } else { + rlld[[ri]] <- log10(lowestdose) + rlhd[[i]] <- log10(highestdose) + runit[[ri]] <- unit + if (inactive) { + mtype[[ri]] <- "inactive" + } else { + if (active) { + mtype[[ri]] <- "active" + } else { + mtype[[ri]] <- "no fit" + } + } + } + sigma[[ri]] <- NA + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + c[[ri]] <- NA + } + } + + 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", + 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) { + m <- models[[row.i]] + 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)) + 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"] + } + } + } + } + } + + attr(results, "models") <- models + return(results) +} +# vim: set ts=4 sw=4 expandtab: @@ -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: @@ -1,4 +1,4 @@ -if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance")) +if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance", "mtype")) drplot <- function(drresults, data, dtype = "std", alpha = 0.95, ctype = "none", path = "./", fileprefix = "drplot", overlay = FALSE, @@ -9,7 +9,8 @@ drplot <- function(drresults, data, postscript = FALSE, pdf = FALSE, png = FALSE, bw = TRUE, pointsize = 12, - colors = 1:8, ltys = 1:8, devoff=TRUE, lpos="topright") + colors = 1:8, ltys = 1:8, pchs = "auto", + devoff=TRUE, lpos="topright") { # Check if all data have the same unit unitlevels <- levels(as.factor(drresults$unit)) @@ -63,6 +64,7 @@ drplot <- function(drresults, data, # Prepare overlay plot if requested if (overlay) { + if (pchs == "auto") pchs = 1:length(dsubstances) if (postscript) { filename = paste(path,fileprefix,".eps",sep="") postscript(file=filename, @@ -91,6 +93,7 @@ drplot <- function(drresults, data, frame.plot = frame.plot) } else { # If overlay plot is not requested, ask before showing multiple plots on the screen + if (pchs == "auto") pchs = rep(1, length(dsubstances)) if (!postscript && !png && !pdf && length(dsubstances) > 1) { op <- par(ask=TRUE) on.exit(par(op)) @@ -111,6 +114,7 @@ drplot <- function(drresults, data, tmp <- splitted[[i]] if (length(tmp$response) != 0) { color <- colors[[n]] + pch <- pchs[[n]] # Prepare the single graphs if an overlay is not requested if (!overlay) { @@ -141,7 +145,7 @@ drplot <- function(drresults, data, axes = axes, frame.plot = frame.plot) } - if (!overlay) legend(lpos, i, lty = 1, col = color, inset=0.05) + if (!overlay) legend(lpos, i, lty = 1, col = color, pch = pch, inset=0.05) tmp$dosefactor <- factor(tmp$dose) # necessary because the old # factor has all levels, not # only the ones tested with @@ -160,7 +164,7 @@ drplot <- function(drresults, data, # Plot the data, if requested if (dtype != "none") { if (dtype == "raw") { - points(log10(tmp$dose),tmp$response,col=color) + points(log10(tmp$dose), tmp$response, col = color, pch = pch) } else { splitresponses <- split(tmp$response,tmp$dosefactor) means <- sapply(splitresponses,mean) @@ -183,7 +187,7 @@ drplot <- function(drresults, data, { x <- log10(as.numeric(levels(tmp$dosefactor))) segments(x,bottoms,x,tops,col=color) - points(x,means,col=color) + points(x, means, col = color, pch = pch) smidge <- 0.05 segments(x - smidge,bottoms,x + smidge,bottoms,col=color) segments(x - smidge,tops,x + smidge,tops,col=color) @@ -192,37 +196,22 @@ drplot <- function(drresults, data, # Plot the fits for this substance, if there are any fits <- subset(drresults,Substance == i) - nf <- length(fits$Substance) # number of fits to plot for this substance + fit.rows <- rownames(fits) + nf <- length(fit.rows) # number of fits to plot for this substance if (nf > 0) { for (j in 1:nf) { - logED50 <- fits[j,"logED50"] - mtype <- as.character(fits[j, "mtype"]) - if (mtype == "probit") { - if (overlay) nl <- nl + 1 else nl = j - lty <- ltys[nl] - scale <- fits[j,"b"] - plot(function(x) pnorm(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty) - } - if (mtype == "logit") { - if (overlay) nl <- nl + 1 else nl = j - lty <- ltys[nl] - scale <- fits[j,"b"] - plot(function(x) plogis(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty) - } - if (mtype == "weibull") { - if (overlay) nl <- nl + 1 else nl = j - lty <- ltys[nl] - location <- fits[j,"a"] - shape <- fits[j,"b"] - plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE, col=color, lty=lty) - } - if (mtype == "linlogit") { - if (overlay) nl <- nl + 1 else nl = j - lty <- ltys[nl] - plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logED50"],fits[j,"b"]), - lld - 0.5, lhd + 2, - add=TRUE, col=color, lty=lty) + fit.row = fit.rows[j] + if (overlay) nl <- nl + 1 else nl = j + lty <- ltys[nl] + if (drresults[[fit.row, "mtype"]] %in% c("probit", "logit", "weibull", "linlogit")) + { + m <- attr(drresults, "models")[[as.numeric(fit.row)]] + of <- function(x) { + predict(m, data.frame(dose = 10^x)) + } + plot(of, lld - 0.5, lhd + 2, + add = TRUE, col = color, lty = lty) } } } @@ -232,7 +221,7 @@ drplot <- function(drresults, data, } } } - if (overlay) legend(lpos, dsubstances, col = colors, lty = ltys, inset=0.05) + if (overlay) legend(lpos, dsubstances, col = colors, pch = pchs, lty = ltys, inset=0.05) if (overlay && (postscript || png || pdf)) { if (devoff) { dev.off() |