diff options
| -rw-r--r-- | DESCRIPTION | 4 | ||||
| -rw-r--r-- | R/drfit.R | 135 | ||||
| -rw-r--r-- | man/drfit.Rd | 84 | ||||
| -rw-r--r-- | man/drplot.Rd | 5 | ||||
| -rwxr-xr-x[-rw-r--r--] | man/linlogitf.Rd (renamed from man/linearlogisf.Rd) | 6 |
5 files changed, 149 insertions, 85 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 3eee40b..6b1b85b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,10 @@ Package: drfit -Version: 0.03-28 +Version: 0.03-29 Date: 2005-07-16 Title: Dose-response data evaluation Author: Johannes Ranke <jranke@uni-bremen.de> Maintainer: Johannes Ranke <jranke@uni-bremen.de> -Depends: R (>=2.1.0),stats,RODBC +Depends: R (>= 2.1.0),stats,RODBC Description: drfit provides basic functions for fitting dose-response curves to dose-response data, calculating some (eco)toxicological parameters and plotting the results. Functions that are fitted are the cumulative densitiy @@ -31,14 +31,14 @@ drdata <- function(substances, experimentator = "%", db = "cytotox", return(data) } -linearlogisf <- function(x,k,f,mu,b) +linlogitf <- function(x,k,f,mu,b) { k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b))) } drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, - lognorm = TRUE, logis = FALSE, - linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA, + probit = TRUE, logit = FALSE, weibull = FALSE, + linlogit = FALSE, linlogitWrong = NA, allWrong = NA, b0 = 2, f0 = 0) { if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data where ok @@ -56,9 +56,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, sigma <- array() # the standard deviation of the residuals logEC50 <- vector() stderrlogEC50 <- vector() - slope <- vector() - b <- vector() - f <- vector() + a <- b <- c <- vector() splitted <- split(data,data$substance) for (i in substances) { @@ -95,10 +93,10 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, if (responseathighestdose < 0.5) { inactive <- FALSE - if (linearlogis && - length(subset(linearlogisWrong,linearlogisWrong == i))==0 && + if (linlogit && + length(subset(linlogitWrong,linlogitWrong == i))==0 && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b), + m <- try(nls(response ~ linlogitf(dose,1,f,logEC50,b), data=tmp, start=list(f=f0,logEC50=startlogEC50[[i]],b=b0))) if (!inherits(m, "try-error")) { @@ -111,27 +109,28 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "linearlogis" + mtype[[ri]] <- "linlogit" logEC50[[ri]] <- coef(m)[["logEC50"]] - slope[[ri]] <- NA if (logEC50[[ri]] > rlhd[[ri]]) { logEC50[[ri]] <- NA stderrlogEC50[[ri]] <- NA + a[[ri]] <- NA b[[ri]] <- NA - f[[ri]] <- NA + c[[ri]] <- NA } else { stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + a[[ri]] <- coef(m)[["logEC50"]] b[[ri]] <- coef(m)[["b"]] - f[[ri]] <- coef(m)[["f"]] + c[[ri]] <- coef(m)[["f"]] } } } - if (logis && + if (weibull && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope), + m <- try(nls(response ~ pweibull(-log10(dose)+location,shape), data=tmp, - start=list(logEC50=startlogEC50[[i]],slope=1))) + start=list(location=-0.16,shape=0.4659))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -140,29 +139,67 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, sigma[[ri]] <- s$sigma rsubstance[[ri]] <- i rn[[ri]] <- n + runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "logis" + mtype[[ri]] <- "weibull" + a[[ri]] <- coef(m)[["location"]] + b[[ri]] <- coef(m)[["shape"]] + sqrdev <- function(logdose) { + (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 + } + logEC50[[ri]] <- nlm(sqrdev,startlogEC50[[i]])$estimate + c[[ri]] <- NA + if (logEC50[[ri]] > rlhd[[ri]]) { + logEC50[[ri]] <- NA + stderrlogEC50[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + } else { + stderrlogEC50[[ri]] <- NA + } + } + } + } + + if (logit && + length(subset(allWrong,allWrong == i))==0) { + m <- try(nls(response ~ plogis(-log10(dose),-logEC50,scale), + data=tmp, + start=list(logEC50=startlogEC50[[i]],scale=1))) + if (chooseone==FALSE || fit==FALSE) { + if (!inherits(m, "try-error")) { + fit <- TRUE + ri <- ri + 1 + s <- summary(m) + sigma[[ri]] <- s$sigma + rsubstance[[ri]] <- i + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + mtype[[ri]] <- "logit" logEC50[[ri]] <- coef(m)[["logEC50"]] - b[[ri]] <- NA - f[[ri]] <- NA + c[[ri]] <- NA if (logEC50[[ri]] > rlhd[[ri]]) { logEC50[[ri]] <- NA - slope[[ri]] <- NA stderrlogEC50[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA } else { - slope[[ri]] <- coef(m)[["slope"]] stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + a[[ri]] <- coef(m)[["logEC50"]] + b[[ri]] <- coef(m)[["scale"]] } } } } - if (lognorm && + if (probit && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope), + m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,scale), data=tmp, - start=list(logEC50=startlogEC50[[i]],slope=1))) + start=list(logEC50=startlogEC50[[i]],scale=1))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -174,17 +211,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "lognorm" + mtype[[ri]] <- "probit" logEC50[[ri]] <- coef(m)[["logEC50"]] - b[[ri]] <- NA - f[[ri]] <- NA + c[[ri]] <- NA if (logEC50[[ri]] > rlhd[[ri]]) { logEC50[[ri]] <- NA - slope[[ri]] <- NA stderrlogEC50[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA } else { - slope[[ri]] <- coef(m)[["slope"]] stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + a[[ri]] <- coef(m)[["logEC50"]] + b[[ri]] <- coef(m)[["scale"]] } } } @@ -215,19 +253,15 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, sigma[[ri]] <- NA logEC50[[ri]] <- NA stderrlogEC50[[ri]] <- NA - slope[[ri]] <- NA + a[[ri]] <- NA b[[ri]] <- NA - f[[ri]] <- NA + c[[ri]] <- NA } } - results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma) - names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma") - if (lognorm || logis) { - results$slope <- slope - } - if (linearlogis) { - results$b <- b - results$f <- f + results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma, a, b) + names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma","a","b") + if (linlogit) { + results$c <- c } return(results) } @@ -257,7 +291,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, } else { lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2 lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2 - if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) { + if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) { hr <- 1.8 } else { hr <- 1.0 @@ -372,16 +406,21 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, { logEC50 <- fits[j,"logEC50"] mtype <- as.character(fits[j, "mtype"]) - if (mtype == "lognorm") { - slope <- fits[j,"slope"] - plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color) + if (mtype == "probit") { + scale <- fits[j,"b"] + plot(function(x) pnorm(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) + } + if (mtype == "logit") { + scale <- fits[j,"b"] + plot(function(x) plogis(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) } - if (mtype == "logis") { - slope <- fits[j,"slope"] - plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color) + if (mtype == "weibull") { + location <- fits[j,"a"] + shape <- fits[j,"b"] + plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color) } - if (mtype == "linearlogis") { - plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]), + if (mtype == "linlogit") { + plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logEC50"],fits[j,"b"]), lld - 0.5, lhd + 2, add=TRUE,col=color) } diff --git a/man/drfit.Rd b/man/drfit.Rd index 2ec4a73..0d9ee6f 100644 --- a/man/drfit.Rd +++ b/man/drfit.Rd @@ -6,58 +6,84 @@ biometric results for (eco)toxicity evaluation } \usage{ - drfit(data, startlogEC50 = NA, chooseone = TRUE, lognorm = TRUE, logis = FALSE, - linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA, b0 = 2, f0 = 0) + drfit(data, startlogEC50 = NA, chooseone = TRUE, probit = TRUE, logit = FALSE, + weibull = FALSE, linlogit = FALSE, linlogitWrong = NA, allWrong = NA, + b0 = 2, f0 = 0) } \arguments{ \item{data}{ - A data frame as returned from \code{\link{drdata}}. The data frame has to + A data frame containing dose-response data. The data frame has to contain at least a factor called "substance", a vector called "unit" containing the unit used for the dose, a column "response" with the response values of the test system normalized between 0 and 1 and a column - "dose" with the numeric dose values. For later use of the - \code{\link{drplot}} function, a factor called "dosefactor" also has to be - present, containing the dose as a factor. If there is a column called "ok" - and it is set to "no fit", then the corresponding data point will be excluded - from the fitting procedure. + "dose" with the numeric dose values. 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 "ok" and it is set to "no fit" in a specific + line, then the corresponding data point will be excluded from the fitting + procedure, although it will be plotted. } \item{startlogEC50}{ - Especially for the linearlogis model, a suitable log10 of the EC50 has to be given + Especially for the linlogit model, a suitable log10 of the EC50 has to be given by the user, since it is not correctly estimated for data showing hormesis with the default estimation method.} - \item{lognorm}{ + \item{probit}{ A boolean defining if cumulative density curves of normal distributions - are fitted to the data. Default ist TRUE.} - \item{logis}{ - A boolean defining if cumulative densitiy curves of logistic distributions - are fitted to the data. Default is FALSE.} - \item{linearlogis}{ + \code{\link{pnorm}} are fitted against the decadic logarithm of the dose. + Default ist TRUE.} + \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.} + \item{weibull}{ + A boolean defining if the cumulative density curves of weibull distributions + (\code{\link{pweibull}} with additionall location parameter and scale=1) + are fitted to the decadic logarithm of the dose. Default is FALSE.} + \item{linlogit}{ A boolean defining if the linear-logistic function - \code{\link{linearlogisf}} as defined by van Ewijk and Hoekstra 1993 is + \code{\link{linlogitf}} as defined by van Ewijk and Hoekstra 1993 is fitted to the data. Default is FALSE.} - \item{linearlogisWrong}{ + \item{linlogitWrong}{ An optional vector containing the names of the substances for which the - linearlogis function produces a wrong fit.} + linlogit function produces a wrong fit.} \item{allWrong}{ An optional vector containing the names of the substances for which all - functions produces a wrong fit.} + functions produce a wrong fit.} \item{chooseone}{ - If TRUE (default), the models are tried in the order linearlogis, logis and lognorm, - and the first model that produces a valid fit is used. Usually this will be the one - with the lowest residual standard deviation. If FALSE, all models that are set to TRUE - and that can be fitted will be reported.} + If TRUE (default), the models are tried in the order linlogit, weibull, + logit, probit, 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{b0,f0}{ If the linearlogistic model is fitted, b0 and f0 give the possibility to adapt the starting values for the parameters b and f.} } \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 "none". - Every successful fit is reported in one line. Parameters of the fitted curves are only - reported if the fitted EC50 is not higher than the highest dose.} - -} + 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 "none". + Every successful fit is reported in one line. Parameters of the fitted + curves are only reported if the fitted EC50 is not higher than the + highest dose. + \code{n} is the number of dose-response curves in the raw data (repetitions + in each point), \code{lld} is the decadic logarithm of the lowest dose and + \code{lhd} is the decadic logarithm of the highest dose. + For the "linlogit", "logit" and "probit" models, the parameter + \code{a} that is reported coincides with the logEC50, i.e the logEC50 is + one of the model parameters that is being fitted, and therefore + a standard deviation \code{std} is reported for the logEC50. In the + case of the "weibull" model, \code{a} is a location parameter. + Parameter \code{b} in the case of the "linlogit" fit is the variable + b from the \code{\link{linlogitf}} function. In the case of "probit" fit + it is the standard deviation of the fitted normal distribution, in the case + of the "logit" fit it is the \code{scale} parameter in the \code{\link{plogis}} + function, and in the "weibull" fit it is the \code{shape} parameter of the + fitted \code{\link{pweibull}} function. Only the "linlogit" fit produces a + third parameter \code{c} which is the variable f from the + \code{\link{linlogitf}} function.} + } \examples{ data(antifoul) r <- drfit(antifoul) diff --git a/man/drplot.Rd b/man/drplot.Rd index ecd6895..e5201a2 100644 --- a/man/drplot.Rd +++ b/man/drplot.Rd @@ -80,9 +80,8 @@ } \note{ - Treatment of legends is not really well solved. Be sure all devices are - closed (e.g. by calling \code{dev.off()}) before calling \code{drplot} - again after a failure. + Be sure all devices are closed (e.g. by calling \code{dev.off()}) before + calling \code{drplot} again after a failure. } \examples{ data(antifoul) diff --git a/man/linearlogisf.Rd b/man/linlogitf.Rd index 1ebe8be..2eca8c4 100644..100755 --- a/man/linearlogisf.Rd +++ b/man/linlogitf.Rd @@ -1,12 +1,12 @@ -\name{linearlogisf} -\alias{linearlogisf} +\name{linlogitf} +\alias{linlogitf} \title{Linear-logistic function} \description{ Helper function describing a special type of dose-response curves, showing a stimulus at subtoxic doses. } \usage{ - linearlogisf(x,k,f,mu,b) + linlogitf(x,k,f,mu,b) } \arguments{ \item{x}{ |