From ec7e9cce71d3b46a84be7268da1de4b0953d54c6 Mon Sep 17 00:00:00 2001 From: ranke Date: Tue, 18 Mar 2014 11:21:20 +0000 Subject: Improvements on the EDx calculations in drcfit. However, in cases that there is no hormesis in the data, drfit often fails to fit the linlogit model, and drcfit often fails to find the EDx values. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@100 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drcfit.R | 34 ++++++++++++++++++---------------- R/drfit.R | 3 +++ R/drplot.R | 5 ++++- 3 files changed, 25 insertions(+), 17 deletions(-) (limited to 'R') diff --git a/R/drcfit.R b/R/drcfit.R index ff9e6bb..8c52f6d 100644 --- a/R/drcfit.R +++ b/R/drcfit.R @@ -37,6 +37,9 @@ drcfit <- function(data, chooseone=TRUE, unit <- "" message("\n",i,": No data\n") } + if (length(unit) == 0) { + unit <- NA + } if (length(unit) > 1) { message("More than one unit for substance ",i,", halting\n\n") break @@ -81,23 +84,24 @@ drcfit <- function(data, chooseone=TRUE, runit[[ri]] <- unit rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) - ED50 <- try(ED(m, 50, interval = "delta", display = FALSE)) + logED50[[ri]] <- NA + logED50low[[ri]] <- NA + logED50high[[ri]] <- NA + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] + c[[ri]] <- coef(m)[[3]] + ED50 <- try(ED(m, 50, interval = "delta", + lower = lowestdose / 10, + upper = highestdose * 10, + display = FALSE)) if (!inherits(ED50, "try-error")) { logED50[[ri]] <- log10(ED50["1:50", "Estimate"]) logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" - logED50[[ri]] <- NA - logED50low[[ri]] <- NA - logED50high[[ri]] <- NA } - } else { - mtype[[ri]] <- "no ED50" - } - a[[ri]] <- coef(m)[[2]] - b[[ri]] <- coef(m)[[1]] - c[[ri]] <- coef(m)[[3]] + } } } } @@ -119,21 +123,19 @@ drcfit <- function(data, chooseone=TRUE, rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) logED50[[ri]] <- log10(coef(m)[[2]]) + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] c[[ri]] <- NA if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA logED50low[[ri]] <- NA logED50high[[ri]] <- NA - a[[ri]] <- NA - b[[ri]] <- NA } else { mtype[[ri]] <- "probit" ED50 <- ED(m, 50, interval = "delta", display = FALSE) logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) - a[[ri]] <- coef(m)[[2]] - b[[ri]] <- coef(m)[[1]] } } } @@ -155,6 +157,8 @@ drcfit <- function(data, chooseone=TRUE, rlld[[ri]] <- log10(lowestdose) rlhd[[ri]] <- log10(highestdose) logED50[[ri]] <- log10(coef(m)[[2]]) + a[[ri]] <- coef(m)[[2]] + b[[ri]] <- coef(m)[[1]] c[[ri]] <- NA if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" @@ -168,8 +172,6 @@ drcfit <- function(data, chooseone=TRUE, ED50 <- ED(m, 50, interval = "delta", display = FALSE) logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) - a[[ri]] <- coef(m)[[2]] - b[[ri]] <- coef(m)[[1]] } } } diff --git a/R/drfit.R b/R/drfit.R index ad7d16c..7fa53bf 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -41,6 +41,9 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, unit <- "" message("\n",i,": No data\n") } + if (length(unit) == 0) { + unit <- NA + } if (length(unit) > 1) { message("More than one unit for substance ",i,", halting\n\n") break diff --git a/R/drplot.R b/R/drplot.R index ded7aca..11c4d42 100644 --- a/R/drplot.R +++ b/R/drplot.R @@ -204,7 +204,10 @@ drplot <- function(drresults, data, fit.row = fit.rows[j] if (overlay) nl <- nl + 1 else nl = j lty <- ltys[nl] - if (drresults[[fit.row, "mtype"]] %in% c("probit", "logit", "weibull", "linlogit")) + if (drresults[[fit.row, "mtype"]] %in% c("probit", + "logit", + "weibull", + "linlogit")) { m <- attr(drresults, "models")[[as.numeric(fit.row)]] of <- function(x) { -- cgit v1.2.1