aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2014-02-25 17:30:01 +0000
committerranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2014-02-25 17:30:01 +0000
commit35362cfc7c02f12fd9689f68a772a804d895535a (patch)
tree7171ac4ee23fa4e56135970936e488b4836e49a8
parent184aacf1ad5a28b2428633cd1966d6fb881eb3b0 (diff)
Feature-complete version of the new drfit package that can make use
of the drc package in order to obtain confidence intervals for EDx values. See ChangeLog for details. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@98 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc
-rw-r--r--ChangeLog24
-rw-r--r--DESCRIPTION9
-rw-r--r--NAMESPACE2
-rw-r--r--R/drcfit.R292
-rw-r--r--R/drfit.R44
-rw-r--r--R/drplot.R57
-rw-r--r--man/IM1xIPC81.Rd7
-rw-r--r--man/IM1xVibrio.Rd8
-rw-r--r--man/antifoul.Rd13
-rw-r--r--man/drcfit.Rd110
-rw-r--r--man/drfit.Rd8
-rw-r--r--man/drplot.Rd7
12 files changed, 508 insertions, 73 deletions
diff --git a/ChangeLog b/ChangeLog
index 9163fbb..79ee431 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -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
diff --git a/NAMESPACE b/NAMESPACE
index 0b39883..4802ab7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -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:
diff --git a/R/drfit.R b/R/drfit.R
index fc73632..ad7d16c 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -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:
diff --git a/R/drplot.R b/R/drplot.R
index 149337c..ded7aca 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -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

Contact - Imprint