diff options
| -rw-r--r-- | DESCRIPTION | 4 | ||||
| -rw-r--r-- | INDEX | 3 | ||||
| -rw-r--r-- | R/drfit.R | 26 | ||||
| -rw-r--r-- | R/drplot.R | 30 | ||||
| -rw-r--r-- | man/drfit.Rd | 12 | ||||
| -rw-r--r-- | man/drplot.Rd | 11 |
6 files changed, 63 insertions, 23 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 4fd580b..fb4798a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drfit -Version: 0.04-60 -Date: 2006-04-01 +Version: 0.04-62 +Date: 2006-04-03 Title: Dose-response data evaluation Author: Johannes Ranke <jranke@uni-bremen.de> Maintainer: Johannes Ranke <jranke@uni-bremen.de> @@ -4,6 +4,9 @@ drdata Get dose-response data drfit Fit dose-response models drfit-package Dose-response data evaluation drplot Plot dose-response models +IM1xIPC81 Dose-Response data for + 1-methyl-3-alkylimidazolium tetrafluoroborates + in IPC-81 cells IM1xVibrio Dose-Response data for 1-methyl-3-alkylimidazolium tetrafluoroborates in V. fischeri @@ -1,6 +1,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, probit = TRUE, logit = FALSE, weibull = FALSE, - linlogit = FALSE, linlogitWrong = NA, allWrong = NA, + linlogit = FALSE, conf = FALSE, + linlogitWrong = NA, allWrong = NA, s0 = 0.5, b0 = 2, f0 = 0) { if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data with @@ -20,6 +21,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, sigma <- array() # the standard deviation of the residuals logED50 <- vector() stderrlogED50 <- vector() + conflogED50 <- vector() a <- b <- c <- vector() splitted <- split(data,data$substance) @@ -47,7 +49,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, n <- ndl <- 0 } else { ndl <- length(levels(factor(tmp$dose))) - n <- round(length(tmp$response)/ndl) + n <- length(tmp$response) if (is.na(startlogED50[i])){ w <- 1/abs(tmp$response - 0.3) startlogED50[[i]] <- sum(w * log10(tmp$dose))/sum(w) @@ -86,12 +88,14 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, mtype[[ri]] <- "no fit" logED50[[ri]] <- NA stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } else { mtype[[ri]] <- "linlogit" stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] + conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 3) a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["b"]] c[[ri]] <- coef(m)[["f"]] @@ -122,11 +126,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, mtype[[ri]] <- "no fit" logED50[[ri]] <- NA stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "probit" stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] + conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2) a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["scale"]] } @@ -158,11 +164,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, mtype[[ri]] <- "no fit" logED50[[ri]] <- NA stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "logit" stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] + conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2) } } } @@ -196,11 +204,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, mtype[[ri]] <- "no fit" logED50[[ri]] <- NA stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "weibull" stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2) } } } @@ -238,13 +248,21 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, sigma[[ri]] <- NA logED50[[ri]] <- NA stderrlogED50[[ri]] <- NA + conflogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } } - results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b) - names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b") + if (conf) + { + results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, conflogED50, runit, sigma, a, b) + names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","conf","unit","sigma","a","b") + } else { + results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b) + names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b") + } + if (linlogit) { results$c <- c } @@ -1,6 +1,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, ctype = "none", path = "./", fileprefix = "drplot", overlay = FALSE, + xlim = c("auto","auto"), ylim = c("auto","auto"), postscript = FALSE, pdf = FALSE, png = FALSE, bw = TRUE, pointsize = 12, @@ -14,14 +15,15 @@ drplot <- function(drresults, data, unit <- "different units" } - # Get the plot limits on the x-axis (log of the dose) and y axis + # Determine the plot limits on the x-axis and y axis if(is.data.frame(data)) { + # Get rid of pseudo substance names of controls nonzerodata <- subset(data,dose!=0) - nonzerodata$substance <- factor(nonzerodata$substance) # Get rid of pseudo substance names of controls + nonzerodata$substance <- factor(nonzerodata$substance) zerodata <- subset(data,dose==0) nc <- length(zerodata$dose) # Number of control points if (nc > 0) { - sdc <- sd(zerodata$response) # Standard deviation of control responses + sdc <- sd(zerodata$response) controlconf <- sdc * qt((1 + alpha)/2, nc - 1) / sqrt(nc) cat("There are ",nc,"data points with dose 0 (control values)\n") cat("with a standard deviation of",sdc,"\n") @@ -52,6 +54,12 @@ drplot <- function(drresults, data, hr <- 1.0 } } + if (xlim[1] == "auto") xlim[1] <- lld - 0.5 + if (xlim[2] == "auto") xlim[2] <- lhd + 1 + if (ylim[1] == "auto") ylim[1] <- -0.1 + if (ylim[2] == "auto") ylim[2] <- hr + 0.2 + xlim <- as.numeric(xlim) + ylim <- as.numeric(ylim) # Prepare overlay plot if requested if (overlay) @@ -79,10 +87,10 @@ drplot <- function(drresults, data, } plot(0,type="n", - xlim=c(lld - 0.5, lhd + 1), - ylim= c(-0.1, hr + 0.2), - xlab=paste("Decadic Logarithm of the dose in ", unit), - ylab="Normalized response") + xlim = xlim, + ylim = ylim, + xlab = paste("Decadic Logarithm of the dose in ", unit), + ylab = "Normalized response") } # Plot the data either as raw data or as error bars @@ -123,10 +131,10 @@ drplot <- function(drresults, data, } plot(0,type="n", - xlim=c(lld - 0.5, lhd + 2), - ylim= c(-0.1, hr + 0.2), - xlab=paste("Decadic Logarithm of the dose in ", unit), - ylab="Normalized response") + xlim = xlim, + ylim = ylim, + xlab = paste("Decadic Logarithm of the dose in ", unit), + ylab = "Normalized response") } if (!overlay) legend(lpos, i,lty = 1, col = color, inset=0.05) tmp$dosefactor <- factor(tmp$dose) # necessary because the old diff --git a/man/drfit.Rd b/man/drfit.Rd index 6e9ef47..9d87b4c 100644 --- a/man/drfit.Rd +++ b/man/drfit.Rd @@ -7,8 +7,8 @@ } \usage{ drfit(data, startlogED50 = NA, chooseone = TRUE, probit = TRUE, logit = FALSE, - weibull = FALSE, linlogit = FALSE, linlogitWrong = NA, allWrong = NA, - s0 = 0.5, b0 = 2, f0 = 0) + weibull = FALSE, linlogit = FALSE, conf = FALSE, linlogitWrong = NA, + allWrong = NA, s0 = 0.5, b0 = 2, f0 = 0) } \arguments{ \item{data}{ @@ -44,6 +44,10 @@ 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.} + \item{conf}{ + A boolean defining if a confidence interval of the log ED50 is to be + listed for alpha = 5 \% (two-sided). If FALSE (default), the standard + deviation is listed in the output of drfit.} \item{linlogitWrong}{ An optional vector containing the names of the substances for which the linlogit function produces a wrong fit.} @@ -73,8 +77,8 @@ 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 - rounded mean of the number of replicates at each dose level in the raw - data, \code{lld} is the decadic logarithm of the lowest dose and + 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. For the \dQuote{linlogit}, \dQuote{logit} and \dQuote{probit} models, the parameter \code{a} that is reported coincides with the logED50, i.e the diff --git a/man/drplot.Rd b/man/drplot.Rd index 4c9e298..0249476 100644 --- a/man/drplot.Rd +++ b/man/drplot.Rd @@ -6,8 +6,9 @@ either combined or separately, for one or more substances. } \usage{ - drplot(drresults, data, dtype, alpha, ctype, path, fileprefix, overlay, - postscript, pdf, png, bw, pointsize, colors, devoff, lpos) + drplot(drresults, data, dtype, alpha, ctype, path, + fileprefix, overlay, xlim, ylim, postscript, pdf, png, bw, + pointsize, colors, devoff, lpos) } \arguments{ \item{drresults}{ @@ -51,6 +52,12 @@ (postscript=FALSE) only works correctly for data plots. Dose-response models will all be put into the last graph in this case. } + \item{xlim}{ + The plot limits (min,max) on the dose axis. + } + \item{ylim}{ + The plot limits (min,max) on the response axis. + } \item{postscript}{ If TRUE, (a) postscript graph(s) will be created. Otherwise, and if the pdf and png arguments are also FALSE, graphics will be |