From 53005a97d56b43a7901f91781145acf30f8bf250 Mon Sep 17 00:00:00 2001 From: ranke Date: Mon, 28 Aug 2006 14:07:59 +0000 Subject: Fix for the weibull fitting. Now there is a warning if the log ED 50 can not be found, and the log ED50 is not reported. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@85 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- DESCRIPTION | 4 ++-- R/drfit.R | 44 ++++++++++++++++++++++++++------------------ 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8f8b1e2..b78e243 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drfit -Version: 0.05-83 -Date: 2006-08-14 +Version: 0.05-85 +Date: 2006-08-28 Title: Dose-response data evaluation Author: Johannes Ranke Maintainer: Johannes Ranke diff --git a/R/drfit.R b/R/drfit.R index e993c61..1414499 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -203,32 +203,40 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, start=list(location=startlogED50[[i]],shape=ws0))) 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 - rndl[[ri]] <- ndl - rn[[ri]] <- n - runit[[ri]] <- unit - rlld[[ri]] <- log10(lowestdose) - rlhd[[ri]] <- log10(highestdose) a[[ri]] <- coef(m)[["location"]] b[[ri]] <- coef(m)[["shape"]] sqrdev <- function(logdose) { (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 } 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 - a[[ri]] <- NA - b[[ri]] <- NA + if (sqrdev(logED50[[ri]]) > 0.1) { + cat("\nCan't find ED50 for fitted weibull model of ",i, + "data\nwith startlogED50", startlogED50[[i]],"\n") + ri <- ri - 1 + length(a) <- length(b) <- ri + length(logED50) <- ri } else { - mtype[[ri]] <- "weibull" + c[[ri]] <- NA + fit <- TRUE + s <- summary(m) + sigma[[ri]] <- s$sigma + rsubstance[[ri]] <- i + rndl[[ri]] <- ndl + rn[[ri]] <- n + runit[[ri]] <- unit + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + if (logED50[[ri]] > rlhd[[ri]]) { + mtype[[ri]] <- "no fit" + logED50[[ri]] <- NA + a[[ri]] <- NA + b[[ri]] <- NA + } else { + mtype[[ri]] <- "weibull" + } } } } -- cgit v1.2.1