From 699448335d6fc5668b569577cb38a628c258eaec Mon Sep 17 00:00:00 2001 From: ranke Date: Wed, 21 Dec 2005 15:33:22 +0000 Subject: Changed the terminology from EC50 to ED50 reflecting that the package is meant for dose-response analysis in general and not only for concentration-response analysis. This is the first revision that has been prepared with the knowledge of the existence of the drc package for dose-response curve analysis. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@48 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drfit.R | 88 +++++++++++++++++++++++++++++++-------------------------------- 1 file changed, 44 insertions(+), 44 deletions(-) (limited to 'R') diff --git a/R/drfit.R b/R/drfit.R index d8631e6..e74667d 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -44,7 +44,7 @@ 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, +drfit <- function(data, startlogED50 = NA, chooseone=TRUE, probit = TRUE, logit = FALSE, weibull = FALSE, linlogit = FALSE, linlogitWrong = NA, allWrong = NA, s0 = 0.5, b0 = 2, f0 = 0) @@ -62,8 +62,8 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, rlhd <- rlld <- vector() # highest and lowest doses tested mtype <- array() # the modeltypes sigma <- array() # the standard deviation of the residuals - logEC50 <- vector() - stderrlogEC50 <- vector() + logED50 <- vector() + stderrlogED50 <- vector() a <- b <- c <- vector() splitted <- split(data,data$substance) @@ -89,9 +89,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, } rix <- ri if (!nodata) { - if (is.na(startlogEC50[i])){ + if (is.na(startlogED50[i])){ w <- 1/abs(tmp$response - 0.3) - startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w) + startlogED50[[i]] <- sum(w * log10(tmp$dose))/sum(w) } highestdose <- max(tmp$dose) lowestdose <- min(tmp$dose) @@ -104,9 +104,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, if (linlogit && length(subset(linlogitWrong,linlogitWrong == i))==0 && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ linlogitf(dose,1,f,logEC50,b), + m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b), data=tmp, - start=list(f=f0,logEC50=startlogEC50[[i]],b=b0))) + start=list(f=f0,logED50=startlogED50[[i]],b=b0))) if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 @@ -117,18 +117,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - logEC50[[ri]] <- coef(m)[["logEC50"]] - if (logEC50[[ri]] > rlhd[[ri]]) { + logED50[[ri]] <- coef(m)[["logED50"]] + if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA + logED50[[ri]] <- NA + stderrlogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } else { mtype[[ri]] <- "linlogit" - stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] - a[[ri]] <- coef(m)[["logEC50"]] + stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] + a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["b"]] c[[ri]] <- coef(m)[["f"]] } @@ -137,9 +137,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, if (probit && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,scale), + m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale), data=tmp, - start=list(logEC50=startlogEC50[[i]],scale=1))) + start=list(logED50=startlogED50[[i]],scale=1))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -151,18 +151,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - logEC50[[ri]] <- coef(m)[["logEC50"]] + logED50[[ri]] <- coef(m)[["logED50"]] c[[ri]] <- NA - if (logEC50[[ri]] > rlhd[[ri]]) { + if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA + logED50[[ri]] <- NA + stderrlogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "probit" - stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] - a[[ri]] <- coef(m)[["logEC50"]] + stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] + a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["scale"]] } } @@ -171,9 +171,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, if (logit && length(subset(allWrong,allWrong == i))==0) { - m <- try(nls(response ~ plogis(-log10(dose),-logEC50,scale), + m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale), data=tmp, - start=list(logEC50=startlogEC50[[i]],scale=1))) + start=list(logED50=startlogED50[[i]],scale=1))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -185,18 +185,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - logEC50[[ri]] <- a[[ri]] <- coef(m)[["logEC50"]] + logED50[[ri]] <- a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["scale"]] c[[ri]] <- NA - if (logEC50[[ri]] > rlhd[[ri]]) { + if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA + logED50[[ri]] <- NA + stderrlogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "logit" - stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"] } } } @@ -206,7 +206,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ pweibull(-log10(dose)+location,shape), data=tmp, - start=list(location=startlogEC50[[i]],shape=s0))) + start=list(location=startlogED50[[i]],shape=s0))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -223,17 +223,17 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, sqrdev <- function(logdose) { (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 } - logEC50[[ri]] <- nlm(sqrdev,startlogEC50[[i]])$estimate + logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate c[[ri]] <- NA - if (logEC50[[ri]] > rlhd[[ri]]) { + if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA + logED50[[ri]] <- NA + stderrlogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "weibull" - stderrlogEC50[[ri]] <- NA + stderrlogED50[[ri]] <- NA } } } @@ -263,15 +263,15 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, } } sigma[[ri]] <- NA - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA + logED50[[ri]] <- NA + stderrlogED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } } - 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") + results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b) + names(results) <- c("Substance","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b") if (linlogit) { results$c <- c } @@ -303,8 +303,8 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, hr <- max(data$response) dsubstances <- levels(data$substance) } else { - lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2 - lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2 + lld <- min(drresults[["logED50"]],na.rm=TRUE) - 2 + lhd <- max(drresults[["logED50"]],na.rm=TRUE) + 2 if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) { hr <- 1.8 } else { @@ -419,15 +419,15 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, if (nf > 0) { for (j in 1:nf) { - logEC50 <- fits[j,"logEC50"] + logED50 <- fits[j,"logED50"] mtype <- as.character(fits[j, "mtype"]) if (mtype == "probit") { scale <- fits[j,"b"] - plot(function(x) pnorm(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) + plot(function(x) pnorm(-x,-logED50,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) + plot(function(x) plogis(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) } if (mtype == "weibull") { location <- fits[j,"a"] @@ -435,7 +435,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95, plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color) } if (mtype == "linlogit") { - plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logEC50"],fits[j,"b"]), + 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) } -- cgit v1.2.1