diff options
author | ranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc> | 2014-03-18 11:21:20 +0000 |
---|---|---|
committer | ranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc> | 2014-03-18 11:21:20 +0000 |
commit | ec7e9cce71d3b46a84be7268da1de4b0953d54c6 (patch) | |
tree | ccb734b11acaf2d3cbfc0ac5a117360e25c714bf /R | |
parent | 189aad5e39df513075fc7c7a64e646b1a354bc4c (diff) |
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
Diffstat (limited to 'R')
-rw-r--r-- | R/drcfit.R | 34 | ||||
-rw-r--r-- | R/drfit.R | 3 | ||||
-rw-r--r-- | R/drplot.R | 5 |
3 files changed, 25 insertions, 17 deletions
@@ -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]] } } } @@ -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 @@ -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) { |