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 --- ChangeLog | 5 +++-- DESCRIPTION | 2 +- R/drcfit.R | 34 ++++++++++++++++++---------------- R/drfit.R | 3 +++ R/drplot.R | 5 ++++- 5 files changed, 29 insertions(+), 20 deletions(-) diff --git a/ChangeLog b/ChangeLog index 79ee431..35a8a35 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-02-25 Johannes Ranke +2014-03-18 Johannes Ranke * DESCRIPTION, NAMESPACE: New verson 0.6.1, with the drc package as dependency. Import needed functionality to avoid masking standard functions. @@ -7,7 +7,8 @@ drfit(), but internally uses the methods supplied by the drc package * R/{drfit,drcfit}.R: Return the list of fitted models in an attribute of - the resulting dataframe + the resulting dataframe. Deal with the case of only NA values in the unit + column of the data. * R/drfit.R: Use the fitted nls() models directly to calculate EDx values diff --git a/DESCRIPTION b/DESCRIPTION index f706688..4bf9df4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: drfit Version: 0.6.1 -Date: 2014-02-25 +Date: 2014-03-18 Title: Dose-response data evaluation Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre"), email = "jranke@uni-bremen.de")) 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