diff options
| -rw-r--r-- | ChangeLog | 24 | ||||
| -rw-r--r-- | DESCRIPTION | 9 | ||||
| -rw-r--r-- | NAMESPACE | 2 | ||||
| -rw-r--r-- | R/drcfit.R | 292 | ||||
| -rw-r--r-- | R/drfit.R | 44 | ||||
| -rw-r--r-- | R/drplot.R | 57 | ||||
| -rw-r--r-- | man/IM1xIPC81.Rd | 7 | ||||
| -rw-r--r-- | man/IM1xVibrio.Rd | 8 | ||||
| -rw-r--r-- | man/antifoul.Rd | 13 | ||||
| -rw-r--r-- | man/drcfit.Rd | 110 | ||||
| -rw-r--r-- | man/drfit.Rd | 8 | ||||
| -rw-r--r-- | man/drplot.Rd | 7 |
12 files changed, 508 insertions, 73 deletions
@@ -1,3 +1,27 @@ +2014-02-25 Johannes Ranke + + * DESCRIPTION, NAMESPACE: New verson 0.6.1, with the drc package as + dependency. Import needed functionality to avoid masking standard functions. + + * R/drcfit.R: Add function drcfit() that imitates the functionality of + drfit(), but internally uses the methods supplied by the drc package + + * R/{drfit,drcfit}.R: Return the list of fitted models in an attribute of + the resulting dataframe + + * R/drfit.R: Use the fitted nls() models directly to calculate EDx values + + * R/drplot.R: Use different point characters (pch) for the data, and add the + possibility to specify them using the argument pchs. Also make drplot() work + for the results from the drcfit function. + + * man/drcfit.Rd: Document the above + + * tests/drcfit.R: Add some tests for this new functionality + + * man/{antifoul,IM1xIPC81,IM1xVibrio}.Rd: Add more example code including + the use of drcfit(). + 2014-02-03 Johannes Ranke * DESCRIPTION: New version 0.5-97 (leading zeros in version numbers are diff --git a/DESCRIPTION b/DESCRIPTION index 7dac5df..f706688 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,10 @@ Package: drfit -Version: 0.5-97 -Date: 2014-02-03 +Version: 0.6.1 +Date: 2014-02-25 Title: Dose-response data evaluation Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre"), email = "jranke@uni-bremen.de")) -Depends: R (>= 2.1.0), MASS, RODBC +Imports: MASS, RODBC, drc Description: drfit provides basic and easy-to-use functions for fitting dose-response curves to dose-response data, calculating some (eco)toxicological parameters and plotting the results. Functions that are @@ -14,7 +14,8 @@ Description: drfit provides basic and easy-to-use functions for fitting derived from the latter, which is used to describe data showing stimulation at low doses (hormesis). In addition, functions checking, plotting and retrieving dose-response data retrieved from a database accessed via RODBC - are included. + are included. As an alternative to the original fitting methods, the + algorithms from the drc package can be used. Encoding: latin1 License: GPL (>= 2) LazyLoad: yes @@ -10,3 +10,5 @@ import( MASS, RODBC ) + +importFrom("drc", drm, ED, LN.2, LL.2, BC.4, W1.2) 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() diff --git a/man/IM1xIPC81.Rd b/man/IM1xIPC81.Rd index 4f5a1dd..dda0c3b 100644 --- a/man/IM1xIPC81.Rd +++ b/man/IM1xIPC81.Rd @@ -21,7 +21,12 @@ the tested organism (name of the cell line). } \examples{ - \dontrun{demo(IM1xIPC81)} + rIM1xIPC81 <- drfit(IM1xIPC81, linlogit = TRUE, showED50 = TRUE, EDx = 10) + + rIM1xIPC81.drc <- drcfit(IM1xIPC81, linlogit = TRUE, showED50 = TRUE, EDx = 10) + + print(rIM1xIPC81,digits=4) + print(rIM1xIPC81.drc,digits=4) } \source{ Ranke J, Mölter K, Stock F, Bottin-Weber U, Poczobutt J, diff --git a/man/IM1xVibrio.Rd b/man/IM1xVibrio.Rd index 5315432..cb0bd48 100644 --- a/man/IM1xVibrio.Rd +++ b/man/IM1xVibrio.Rd @@ -19,7 +19,13 @@ regarded valid (\code{ok}). } \examples{ - \dontrun{demo(IM1xVibrio)} + rIM1xVibrio <- drfit(IM1xVibrio, logit = TRUE, chooseone = FALSE, + showED50 = TRUE, EDx = c(10, 20)) + print(rIM1xVibrio, digits = 4) + + rIM1xVibrio.drc <- drcfit(IM1xVibrio, logit = TRUE, chooseone = FALSE, + showED50 = TRUE, EDx = c(10, 20)) + print(rIM1xVibrio.drc, digits = 4) } \source{ Ranke J, Mölter K, Stock F, Bottin-Weber U, Poczobutt J, Hoffmann J, diff --git a/man/antifoul.Rd b/man/antifoul.Rd index 2aa532e..21ae885 100644 --- a/man/antifoul.Rd +++ b/man/antifoul.Rd @@ -15,8 +15,19 @@ database are also present. } \examples{ - \dontrun{demo(antifoul)} +rantifoul.ED50 <- drfit(antifoul, + linlogit = TRUE, logit = TRUE, weibull = TRUE, + chooseone = FALSE, + showED50 = TRUE, EDx = c(10)) +print(rantifoul.ED50, digits = 5) + +rantifoul.drc <- drcfit(antifoul, + linlogit = TRUE, logit = TRUE, weibull = TRUE, + chooseone = FALSE, + showED50 = TRUE, EDx = c(10)) +print(rantifoul.drc, digits = 5) } + \source{ \url{http://www.uft.uni-bremen.de/chemie} } diff --git a/man/drcfit.Rd b/man/drcfit.Rd new file mode 100644 index 0000000..2da0418 --- /dev/null +++ b/man/drcfit.Rd @@ -0,0 +1,110 @@ +\name{drcfit} +\alias{drcfit} +\title{Fit dose-response models using the drc package} +\description{ + Fit dose-response relationships to dose-response data and calculate + biometric results for (eco)toxicity evaluation using the drc package +} +\usage{ + drcfit(data, chooseone = TRUE, probit = TRUE, logit = FALSE, + weibull = FALSE, linlogit = FALSE, level = 0.95, + showED50 = FALSE, EDx = NULL) +} +\arguments{ + \item{data}{ + A data frame containing dose-response data. The data frame has to contain + at least a factor called \dQuote{substance}, a numeric vector \dQuote{dose} + with the dose values, a vector called \dQuote{unit} containing the unit + used for the dose and a numeric vector \dQuote{response} with the response + values of the test system normalized between 0 and 1. Such a data frame can + be easily obtained if a compliant RODBC data source is available for use in + conjunction with the function \code{\link{drdata}}. + + If there is a column called \dQuote{ok} and it is set to \dQuote{no fit} in + a specific line, then the corresponding data point will be excluded from + the fitting procedure, although it will be plotted.} + \item{probit}{ + A boolean defining if cumulative density curves of normal distributions + \code{\link{pnorm}} are fitted against the decadic logarithm of the dose. + Default ist TRUE. + Note that the parameter definitions used in the model are different to the + ones used in \code{\link{drfit}}. Parameter e from \code{\link{LN.2}} is listed + as a here, and parameter b from LN.2 is listed as b.} + \item{logit}{ + A boolean defining if cumulative density curves of logistic distributions + \code{\link{plogis}} are fitted to the decadic logarithm of the dose. + Default is FALSE. + Again the parameter definitions used in the model are different to the + ones used in \code{\link{drfit}}. Parameter e from \code{\link{LL.2}} is listed + as a here, and parameter b from LL.2 is listed as b.} + \item{weibull}{ + A boolean defining if Weibull dose-response models + (\code{\link{W1.2}} are fitted to the untransformed dose. Default is FALSE. + Note that the results differ from the ones obtained with + \code{\link{drfit}}, due to a different model specification.} + \item{linlogit}{ + A boolean defining if the linear-logistic function + \code{\link{linlogitf}} as defined by van Ewijk and Hoekstra 1993 is + fitted to the data. Default is FALSE. Obtaining the ED50 (and EDx values + in general) uses \code{\link{ED}} internally and does not always give a + result. + } + \item{level}{ + The level for the confidence interval listed for the log ED50.} + \item{chooseone}{ + If TRUE (default), the models are tried in the order linlogit, probit, + logit, weibull, and the first model that produces a valid fit is used. + If FALSE, all models that are set to TRUE and that can be fitted will be + reported.} + \item{EDx}{ + A vector of inhibition values x in percent for which the corresponding doses + EDx should be reported. + } + \item{showED50}{ + If set to TRUE, the ED50 and its confidence interval on the original dose + scale (not log scale) is included in the output. + } +} +\value{ + \item{results}{ + A data frame containing at least one line for each substance. If the data + did not show a mean response < 0.5 at the highest dose level, the + modeltype is set to \dQuote{inactive}. If the mean response at the lowest + dose is smaller than 0.5, the modeltype is set to \dQuote{active}. In + both cases, no fitting procedure is carried out. Every successful fit is + reported in one line. Parameters of the fitted curves are only reported + if the fitted ED50 is not higher than the highest dose. + + \code{ndl} is the number of dose levels in the raw data, \code{n} is the + total number of data points in the raw data used for the fit + \code{lld} is the decadic logarithm of the lowest dose and + \code{lhd} is the decadic logarithm of the highest dose. + + If the parameter \code{showED50} was set to TRUE, the ED50 values and their + confidence intervals are also included on the original dose scale. + + If one or more response leves were specified in the argument \code{EDx}, + the corresponding dose levels are given, together with their confidence + intervals. + } +} +\examples{ +data(antifoul) +r <- drcfit(antifoul, showED50 = TRUE, EDx = c(5, 10, 20)) +format(r, digits = 2) +} +\note{There is a demo for each dataset that can be accessed by + \code{demo(dataset)}} +\seealso{ + Further examples are given in help pages to the datasets + \code{\link{antifoul}}, \code{\link{IM1xIPC81}} and + \code{\link{IM1xVibrio}}. +} +\author{ + Johannes Ranke + \email{jranke@uni-bremen.de} + \url{http://www.uft.uni-bremen.de/chemie/ranke} +} +\keyword{models} +\keyword{regression} +\keyword{nonlinear} diff --git a/man/drfit.Rd b/man/drfit.Rd index 7237404..fe84a91 100644 --- a/man/drfit.Rd +++ b/man/drfit.Rd @@ -114,6 +114,9 @@ If the parameter \code{showED50} was set to TRUE, the ED50 values and their confidence intervals are also included on the original dose scale. + + If one or more response leves were specified in the argument \code{EDx}, + the corresponding dose levels are given in addition. } } \examples{ @@ -123,6 +126,11 @@ format(r, digits = 2) } \note{There is a demo for each dataset that can be accessed by \code{demo(dataset)}} +\seealso{ + Further examples are given in help pages to the datasets + \code{\link{antifoul}}, \code{\link{IM1xIPC81}} and + \code{\link{IM1xVibrio}}. +} \author{ Johannes Ranke \email{jranke@uni-bremen.de} diff --git a/man/drplot.Rd b/man/drplot.Rd index b2c2745..e9e4dd1 100644 --- a/man/drplot.Rd +++ b/man/drplot.Rd @@ -8,7 +8,7 @@ \usage{ drplot(drresults, data, dtype, alpha, ctype, path, fileprefix, overlay, xlim, ylim, xlab, ylab, axes, frame.plot, postscript, - pdf, png, bw, pointsize, colors, ltys, devoff, lpos) + pdf, png, bw, pointsize, colors, ltys, pchs, devoff, lpos) } \arguments{ \item{drresults}{ @@ -99,6 +99,11 @@ \item{ltys}{ This is a vector of line types for the dose-response models, defaulting to 1:8. } + \item{pchs}{ + This is a vector of character types for the data. The default uses built-in + characters 1:n with n being the number of substances for which data are plotted + for overlays, or always 1 when no overlay plot is generated. + } \item{lpos}{ An optional argument defaulting to "topright" specifying the position of the legend by being passed to the legend function. See the help for the |