From b83cc0e0ec978b5d047924a222a387152b086f4b Mon Sep 17 00:00:00 2001 From: ranke Date: Sat, 16 Jul 2005 09:12:22 +0000 Subject: Several changes in this revision: - Addition of the weibull fit - Renaming of lognorm to probit, logis to logit, and linearlogis to linlogit - Fix of a bug that prevented the logis fit to work in recent version - Different way of listing parameters in the output. - Documentation of changes git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@29 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drfit.R | 135 ++++++++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 87 insertions(+), 48 deletions(-) (limited to 'R') 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) } -- cgit v1.2.1