From 2cf977cb94fcf6e68b4a75cdf33a29a5308b9457 Mon Sep 17 00:00:00 2001 From: ranke Date: Wed, 3 May 2006 22:53:31 +0000 Subject: Changed the way confidence intervals are being calculated. Now the function confint.nls from the MASS package is being used, after stumbling across the relevant paragraph of the MASS book. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@76 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drfit.R | 97 ++++++++++++++++++++++++++++++++++++--------------------------- 1 file changed, 55 insertions(+), 42 deletions(-) (limited to 'R') diff --git a/R/drfit.R b/R/drfit.R index 1bbffa8..b33d754 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -1,27 +1,28 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, probit = TRUE, logit = FALSE, weibull = FALSE, - linlogit = FALSE, conf = FALSE, + linlogit = FALSE, level = 0.95, 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 - # ok set to "no fit" + require(MASS) + if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data + # with ok set to + # "no fit" substances <- levels(data$substance) - ri <- rix <- 0 # ri is the index over the result rows - # rix is used later to check if any - # model result was appended - rsubstance <- array() # the substance names in the results - rndl <- vector() # number of dose levels - rn <- vector() # mean number of replicates - # in each dose level - runit <- vector() # vector of units for each result row - rlhd <- rlld <- vector() # highest and lowest doses tested - mtype <- array() # the modeltypes - sigma <- array() # the standard deviation of the residuals + ri <- rix <- 0 # ri is the index over the result rows + # rix is used later to check if any + # model result was appended + rsubstance <- array() # the substance names in the results + rndl <- vector() # number of dose levels + rn <- vector() # mean number of replicates + # in each dose level + runit <- vector() # vector of units for each result row + rlhd <- rlld <- vector() # highest and lowest doses tested + mtype <- array() # the modeltypes + sigma <- array() # the standard deviation of the residuals logED50 <- vector() - stderrlogED50 <- vector() - conflogED50 <- vector() + logED50low <- logED50high <- vector() a <- b <- c <- vector() splitted <- split(data,data$substance) @@ -87,15 +88,21 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA - stderrlogED50[[ri]] <- NA - conflogED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[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) + logED50conf <- try(confint(m,"logED50",level=level)) + if (!inherits(logED50conf, "try-error")) { + logED50low[[ri]] <- logED50conf[[1]] + logED50high[[ri]] <- logED50conf[[2]] + } else { + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + } a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["b"]] c[[ri]] <- coef(m)[["f"]] @@ -125,14 +132,20 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA - stderrlogED50[[ri]] <- NA - conflogED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[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) + logED50conf <- try(confint(m,"logED50",level=level)) + if (!inherits(logED50conf, "try-error")) { + logED50low[[ri]] <- logED50conf[[1]] + logED50high[[ri]] <- logED50conf[[2]] + } else { + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + } a[[ri]] <- coef(m)[["logED50"]] b[[ri]] <- coef(m)[["scale"]] } @@ -163,14 +176,20 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA - stderrlogED50[[ri]] <- NA - conflogED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[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) + logED50conf <- try(confint(m,"logED50",level=level)) + if (!inherits(logED50conf, "try-error")) { + logED50low[[ri]] <- logED50conf[[1]] + logED50high[[ri]] <- logED50conf[[2]] + } else { + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + } } } } @@ -200,17 +219,15 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, } logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate c[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA if (logED50[[ri]] > rlhd[[ri]]) { 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) } } } @@ -247,21 +264,17 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, } sigma[[ri]] <- NA logED50[[ri]] <- NA - stderrlogED50[[ri]] <- NA - conflogED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } } - 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") - } + results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, + logED50, logED50low, logED50high, runit, sigma, a, b) + names(results) <- c("Substance","ndl","n","lld","lhd","mtype", + "logED50",names(logED50conf)[[1]],names(logED50conf)[[2]],"unit","sigma","a","b") if (linlogit) { results$c <- c -- cgit v1.2.1