aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--DESCRIPTION4
-rw-r--r--R/drfit.R135
-rw-r--r--man/drfit.Rd84
-rw-r--r--man/drplot.Rd5
-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
diff --git a/R/drfit.R b/R/drfit.R
index 38f8e75..603cf94 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -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}{

Contact - Imprint