From baf8e4341b58c09532b4ca37ebf80c9068934db4 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 24 Mar 2017 15:28:44 +0100 Subject: Import version 0.6.7 from CRAN I had deleted the git commits between 0.6.4 and 0.6.7 while removing drfit from github, while the latest changes were not pushed to jrwb.de yet. --- ChangeLog | 8 ++++++ DESCRIPTION | 24 ++++++++++------- NAMESPACE | 3 +-- R/checkexperiment.R | 28 +++++++++---------- R/checksubstance.R | 52 ++++++++++++++++++------------------ R/drcfit.R | 68 ++++++++++++++++++++++++++--------------------- R/drdata.R | 2 +- R/drfit.R | 47 ++++++++++++++------------------ R/drplot.R | 36 ++++++++++++------------- data/XY.rda | Bin 1602 -> 630 bytes inst/doc/drfit-Rnews.pdf | Bin 229223 -> 169560 bytes 11 files changed, 140 insertions(+), 128 deletions(-) diff --git a/ChangeLog b/ChangeLog index dce0368..170eab9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2016-09-02 Johannes Ranke + + * Make drcfit compatible with drc >= 3.0.1 (the row names of the matrix + returned by the ED() function were changed). + * Switch back to using requireNamespace(RODBC) in order to address + a NOTE issued by R CMD check + * Use single quotes to quote other packages in the DESCRIPTION file + 2015-12-15 Johannes Ranke * Use require(RODBC) instead of requireNamespace(RODBC) as the former fails diff --git a/DESCRIPTION b/DESCRIPTION index 173e9c6..4a0f79b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,10 @@ Package: drfit -Version: 0.7.1 -Date: 2017-03-24 +Version: 0.6.7 +Date: 2016-09-02 Title: Dose-Response Data Evaluation Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre"), email = "jranke@uni-bremen.de")) -Imports: graphics, grDevices, MASS, drc, reshape2, qcc +Imports: graphics, grDevices, MASS, drc Suggests: RODBC Description: A somewhat outdated package of basic and easy-to-use functions for fitting dose-response curves to continuous dose-response data, calculating some @@ -12,16 +12,22 @@ Description: A somewhat outdated package of basic and easy-to-use functions for the more powerful and actively developed 'drc' package. Functions that are fitted are the cumulative density function of the lognormal distribution (probit fit), of the logistic distribution (logit fit), of the weibull - distribution (weibull fit) and a linear-logistic model ("linlogit" fit), + distribution (weibull fit) and a linear-logistic model ('linlogit' fit), derived from the latter, which is used to describe data showing stimulation at low doses (hormesis). In addition, functions checking, plotting and retrieving - dose-response data retrieved from a database accessed via RODBC are included. - As an alternative to the original fitting methods, the algorithms from the drc + dose-response data retrieved from a database accessed via 'RODBC' are included. + As an alternative to the original fitting methods, the algorithms from the 'drc' package can be used. Encoding: latin1 License: GPL (>= 2) LazyLoad: yes LazyData: yes -URL: http://www.r-project.org, - http://www.uft.uni-bremen.de/chemie/ranke?page=drfit, - http://cgit.jrwb.de/drfit/ +URL: http://www.r-project.org, + http://www.uft.uni-bremen.de/chemie/ranke?page=drfit, + http://cgit.jrwb.de/drfit/ +NeedsCompilation: no +Packaged: 2016-09-02 14:24:52 UTC; jranke +Author: Johannes Ranke [aut, cre] +Maintainer: Johannes Ranke +Repository: CRAN +Date/Publication: 2016-09-02 21:28:14 diff --git a/NAMESPACE b/NAMESPACE index 120d052..cadcc62 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,5 +8,4 @@ import( ) importFrom("drc", drm, ED, LN.2, LL.2, BC.4, W1.2) -importFrom("reshape2", melt, acast) -importFrom("qcc", qcc) +importFrom("utils", "packageVersion") diff --git a/R/checkexperiment.R b/R/checkexperiment.R index b69b81f..8c2f472 100644 --- a/R/checkexperiment.R +++ b/R/checkexperiment.R @@ -25,21 +25,21 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") " WHERE ", exptype, " = ", id) commentdata <- RODBC::sqlQuery(channel,commentquery) comment <- as.character(commentdata[[1]]) - + expquery <- paste("SELECT experimentator,substance, ", testtype, ",conc,unit,", responsename, ",performed,ok", - " FROM ",db," WHERE ",exptype,"=", id, + " FROM ",db," WHERE ",exptype,"=", id, sep = "") if (db == "ecotox") { - expquery <- paste(expquery, " AND type LIKE '", + expquery <- paste(expquery, " AND type LIKE '", endpoint, "'", sep = "") } expdata <- RODBC::sqlQuery(channel,expquery) if (db %in% c("cytotox","enzymes")) { - controlquery <- paste("SELECT type,response FROM controls + controlquery <- paste("SELECT type,response FROM controls WHERE plate=",id) controldata <- RODBC::sqlQuery(channel,controlquery) } @@ -63,7 +63,7 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") numberOfBlinds <- NA meanOfBlinds <- NA stdOfBlinds <- NA - + } numberOfControls <- length(controls$response) if (numberOfControls > 0) { @@ -78,7 +78,7 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") if (length(expdata$experimentator) < 1) { stop("There is no response data for ",exptype," ", id," in database ",db,"\n") - } + } exptypestring <- paste(toupper(substring(exptype,1,1)), substring(exptype,2),sep="") expdata$experimentator <- factor(expdata$experimentator) @@ -87,7 +87,7 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") expdata$substance <- factor(expdata$substance) expdata$unit <- factor(expdata$unit) expdata$ok <- factor(expdata$ok) - + cat("\n",exptypestring,id,"from database",db,":\n\n", "\tExperimentator(s):\t",levels(expdata$experimentator),"\n", "\tType(s):\t\t",levels(expdata$type),"\n", @@ -100,7 +100,7 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") "\tblind\t",numberOfBlinds,"\t",meanOfBlinds,"\t",stdOfBlinds,"\n", "\tcontrol\t",numberOfControls,"\t",meanOfControls,"\t", stdOfControls,"\t\t",percentstdOfcontrols,"\n") - + if (db == "ecotox") { boxplot(controls$response, @@ -116,22 +116,22 @@ checkexperiment <- function(id, db = "ecotox", endpoint = "%") boxwex=0.4, main=paste("Plate ",id)) } - + drdata <- expdata[c(2,4,6)] drdata$substance <- factor(drdata$substance) substances <- levels(drdata$substance) - + lld <- log10(min(subset(drdata,conc!=0)$conc)) lhd <- log10(max(drdata$conc)) plot(1,type="n", - xlim=c(lld - 0.5, lhd + 2), - ylim= c(-0.1, 2), + xlim=c(lld - 0.5, lhd + 2), + ylim= c(-0.1, 2), xlab=paste("decadic logarithm of the concentration in ",levels(expdata$unit)), ylab=responsename) - + drdatalist <- split(drdata,drdata$substance) - + for (i in 1:length(drdatalist)) { points(log10(drdatalist[[i]]$conc),drdatalist[[i]][[responsename]],col=i); } diff --git a/R/checksubstance.R b/R/checksubstance.R index 06cc635..3e07f92 100644 --- a/R/checksubstance.R +++ b/R/checksubstance.R @@ -1,7 +1,7 @@ checksubstance <- function(substance, db = "cytotox", experimentator = "%", - celltype = "%", enzymetype = "%", organism = "%", + celltype = "%", enzymetype = "%", organism = "%", endpoint = "%", - whereClause = "1", ok= "%") + whereClause = "1", ok= "%") { databases <- data.frame( responsename=c("viability","activity","response"), @@ -23,13 +23,13 @@ checksubstance <- function(substance, db = "cytotox", experimentator = "%", if (db == "cytotox") { type <- celltype - } + } if (db == "enzymes") { type <- enzymetype - } + } if (db == "ecotox") { type <- organism - } + } query <- paste("SELECT experimentator,substance,", testtype, ",", exptype, ",conc,unit,",responsename,",ok", @@ -41,7 +41,7 @@ checksubstance <- function(substance, db = "cytotox", experimentator = "%", sep = "") if (db == "ecotox") { - query <- paste(query, " AND type LIKE '", + query <- paste(query, " AND type LIKE '", endpoint, "'", sep = "") } @@ -51,39 +51,39 @@ checksubstance <- function(substance, db = "cytotox", experimentator = "%", if (length(data$experimentator) < 1) { stop(paste("\nNo response data for",substance,"in database", db,"found with these parameters\n")) - } - - data$experimentator <- factor(data$experimentator) + } + + data$experimentator <- factor(data$experimentator) data$substance <- factor(data$substance) - substances <- levels(data$substance) - data$type <- factor(data[[testtype]]) - data[[exptype]] <- factor(data[[exptype]]) + substances <- levels(data$substance) + data$type <- factor(data[[testtype]]) + data[[exptype]] <- factor(data[[exptype]]) experiments <- levels(data[[exptype]]) concentrations <- split(data$conc,data$conc) concentrations <- as.numeric(names(concentrations)) - data$unit <- factor(data$unit) + data$unit <- factor(data$unit) data$ok <- factor(data$ok) if (length(experiments)>6) { - palette(rainbow(length(experiments))) + palette(rainbow(length(experiments))) } - + plot(log10(data$conc),data[[responsename]], - xlim=c(-2.5, 4.5), - ylim= c(-0.1, 2), - xlab=paste("decadic logarithm of the concentration in ",levels(data$unit)), - ylab=responsename) - + xlim=c(-2.5, 4.5), + ylim= c(-0.1, 2), + xlab=paste("decadic logarithm of the concentration in ",levels(data$unit)), + ylab=responsename) + explist <- split(data,data[[exptype]]) - - for (i in 1:length(explist)) { - points(log10(explist[[i]]$conc),explist[[i]][[responsename]],col=i); - } - + + for (i in 1:length(explist)) { + points(log10(explist[[i]]$conc),explist[[i]][[responsename]],col=i); + } + legend("topleft", experiments, pch=1, col=1:length(experiments), inset=0.05) title(main=paste(substance," - ",levels(data$experimentator)," - ",levels(data$type))) - exptypename <- paste(toupper(substring(exptype,1,1)), + exptypename <- paste(toupper(substring(exptype,1,1)), substring(exptype,2), sep = "") experimentators <- paste(levels(data$experimentator), collapse = " ") types <- paste(levels(data$type), collapse = " ") diff --git a/R/drcfit.R b/R/drcfit.R index 0d27c1c..008c53c 100644 --- a/R/drcfit.R +++ b/R/drcfit.R @@ -14,7 +14,7 @@ drcfit <- function(data, chooseone=TRUE, # model result was appended rsubstance <- array() # the substance names in the results rndl <- vector() # number of dose levels - rn <- vector() # mean number of replicates + 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 @@ -27,6 +27,16 @@ drcfit <- function(data, chooseone=TRUE, models <- list() # a list containing the dose-response models splitted <- split(data,data$substance) + + # The indexing of the results of the ED50 function changed with drc 3.0.1 + if (packageVersion("drc") > 3) { + ED50_row_index = "e:1:50" + EDx_row_index_prefix = "e:1" + } else { + ED50_row_index = "1:50" + EDx_row_index_prefix = "1" + } + for (i in substances) { tmp <- splitted[[i]] fit <- FALSE @@ -69,8 +79,7 @@ drcfit <- function(data, chooseone=TRUE, active <- FALSE if (linlogit) { - m <- try(drm(response ~ dose, data = tmp, fct = BC.4(fixed = c(NA, 1, NA, NA))), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = BC.4(fixed = c(NA, 1, NA, NA)))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -91,26 +100,25 @@ drcfit <- function(data, chooseone=TRUE, a[[ri]] <- coef(m)[[2]] b[[ri]] <- coef(m)[[1]] c[[ri]] <- coef(m)[[3]] - ED50 <- try(ED(m, 50, interval = "delta", + ED50 <- try(ED(m, 50, interval = "delta", lower = lowestdose / 10, upper = highestdose * 10, - display = FALSE), - silent = TRUE) + 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"]) + logED50[[ri]] <- log10(ED50[ED50_row_index, "Estimate"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" } - } + } } } } if (probit) { - m <- try(drm(response ~ dose, data = tmp, - fct = LN.2()), silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, + fct = LN.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -136,16 +144,15 @@ drcfit <- function(data, chooseone=TRUE, } 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"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) } } } } if (logit) { - m <- try(drm(response ~ dose, data = tmp, fct = LL.2()), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = LL.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -173,8 +180,8 @@ drcfit <- function(data, chooseone=TRUE, } else { mtype[[ri]] <- "logit" ED50 <- ED(m, 50, interval = "delta", display = FALSE) - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) } } } @@ -182,8 +189,7 @@ drcfit <- function(data, chooseone=TRUE, } if (weibull) { - m <- try(drm(response ~ dose, data = tmp, fct = W1.2()), - silent = TRUE) + m <- try(drm(response ~ dose, data = tmp, fct = W1.2())) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -199,7 +205,7 @@ drcfit <- function(data, chooseone=TRUE, rlhd[[ri]] <- log10(highestdose) c[[ri]] <- NA ED50 <- ED(m, 50, interval = "delta", display = FALSE) - logED50[[ri]] <- log10(ED50["1:50", "Estimate"]) + logED50[[ri]] <- log10(ED50[ED50_row_index, "Estimate"]) if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA @@ -209,8 +215,8 @@ drcfit <- function(data, chooseone=TRUE, b[[ri]] <- NA } else { mtype[[ri]] <- "weibull" - logED50low[[ri]] <- log10(ED50["1:50", "Lower"]) - logED50high[[ri]] <- log10(ED50["1:50", "Upper"]) + logED50low[[ri]] <- log10(ED50[ED50_row_index, "Lower"]) + logED50high[[ri]] <- log10(ED50[ED50_row_index, "Upper"]) a[[ri]] <- logED50[[ri]] b[[ri]] <- coef(m)[[1]] } @@ -256,7 +262,7 @@ drcfit <- function(data, chooseone=TRUE, } } - results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, + results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, logED50low, logED50high, runit, sigma, a, b) lower_level_percent = paste(100 * (1 - level)/2, "%", sep = "") upper_level_percent = paste(100 * (1 + level)/2, "%", sep = "") @@ -279,14 +285,14 @@ drcfit <- function(data, chooseone=TRUE, mtype <- as.character(results[row.i, "mtype"]) if (mtype %in% c("probit", "logit", "weibull", "linlogit")) { for (EDi in EDx) { - EDx.drc = try(ED(m, EDi, interval = "delta", display = FALSE, level = level), - silent = TRUE) + EDx.drc = try(ED(m, EDi, interval = "delta", display = FALSE, level = level)) if (!inherits(EDx.drc, "try-error")) { - results[row.i, paste0("EDx", EDi)] <- EDx.drc[paste0("1:", EDi), "Estimate"] - results[row.i, paste0("EDx", EDi, " ", lower_level_percent)] <- EDx.drc[paste0("1:", EDi), - "Lower"] - results[row.i, paste0("EDx", EDi, " ", upper_level_percent)] <- EDx.drc[paste0("1:", EDi), - "Upper"] + results[row.i, paste0("EDx", EDi)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Estimate"] + results[row.i, paste0("EDx", EDi, " ", lower_level_percent)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Lower"] + results[row.i, paste0("EDx", EDi, " ", upper_level_percent)] <- + EDx.drc[paste(EDx_row_index_prefix, EDi, sep = ":"), "Upper"] } } } diff --git a/R/drdata.R b/R/drdata.R index dc741a6..0bf9597 100644 --- a/R/drdata.R +++ b/R/drdata.R @@ -21,7 +21,7 @@ drdata <- function(substances, experimentator = "%", db = "cytotox", type <- organism } } - + query <- paste("SELECT conc,",responsetype,",unit,experimentator,substance,",testtype, ",ok FROM ", db, " WHERE substance IN ('", slist,"') AND experimentator LIKE '", diff --git a/R/drfit.R b/R/drfit.R index 65ab8e4..7fa53bf 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -18,7 +18,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, # model result was appended rsubstance <- array() # the substance names in the results rndl <- vector() # number of dose levels - rn <- vector() # mean number of replicates + 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 @@ -75,13 +75,12 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, active <- TRUE } else { active <- FALSE - if (linlogit && + if (linlogit && length(subset(linlogitWrong,linlogitWrong == i))==0 && length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b), - data=tmp, algorithm="port", - start=list(f=f0,logED50=startlogED50[[i]],b=b0)), - silent = TRUE) + data=tmp, algorithm="port", + start=list(f=f0,logED50=startlogED50[[i]],b=b0))) if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 @@ -105,8 +104,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, c[[ri]] <- NA } else { mtype[[ri]] <- "linlogit" - logED50conf <- try(confint(m,"logED50",level=level), - silent = TRUE) + logED50conf <- try(confint(m,"logED50",level=level)) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] @@ -124,9 +122,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (probit && length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale), - data=tmp, algorithm="port", - start=list(logED50=startlogED50[[i]],scale=ps0)), - silent = TRUE) + data=tmp, algorithm="port", + start=list(logED50=startlogED50[[i]],scale=ps0))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -151,8 +148,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, b[[ri]] <- NA } else { mtype[[ri]] <- "probit" - logED50conf <- try(confint(m,"logED50",level=level), - silent = TRUE) + logED50conf <- try(confint(m,"logED50",level=level)) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] @@ -170,9 +166,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (logit && length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale), - data=tmp, algorithm="port", - start=list(logED50=startlogED50[[i]],scale=ls0)), - silent = TRUE) + data=tmp, algorithm="port", + start=list(logED50=startlogED50[[i]],scale=ls0))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE @@ -198,8 +193,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, b[[ri]] <- NA } else { mtype[[ri]] <- "logit" - logED50conf <- try(confint(m,"logED50",level=level), - silent = TRUE) + logED50conf <- try(confint(m,"logED50",level=level)) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] @@ -215,9 +209,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (weibull && length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ pweibull(-log10(dose)+location,shape), - data=tmp, algorithm="port", - start=list(location=startlogED50[[i]],shape=ws0)), - silent = TRUE) + data=tmp, algorithm="port", + start=list(location=startlogED50[[i]],shape=ws0))) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { ri <- ri + 1 @@ -225,7 +218,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, a[[ri]] <- coef(m)[["location"]] b[[ri]] <- coef(m)[["shape"]] sqrdev <- function(logdose) { - (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 + (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 } logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate if (sqrdev(logED50[[ri]]) > 0.1) { @@ -298,7 +291,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, c[[ri]] <- NA } } - results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, + results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, logED50low, logED50high, runit, sigma, a, b) lower_level_percent = paste(100 * (1 - level)/2, "%", sep = "") upper_level_percent = paste(100 * (1 + level)/2, "%", sep = "") @@ -318,18 +311,18 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE, if (mtype[[row.i]] %in% c("probit", "logit", "weibull", "linlogit")) { for (ED in EDx) { of <- function(x) { - abs(predict(models[[row.i]], data.frame(dose = 10^x)) - + abs(predict(models[[row.i]], data.frame(dose = 10^x)) - (1 - (ED/100))) - } - # Search over interval starting an order of magnitude below + } + # Search over interval starting an order of magnitude below # the lowest dose up to one order of magnitude above the # highest dose - o = optimize(of, + o = optimize(of, results[row.i, c("lld", "lhd")] + c(-1, 1)) # Only keep results within the tolerance if ((o$objective) < EDx.tolerance) { logdose.ED = o$minimum - results[row.i, paste0("EDx", ED)] <- 10^logdose.ED + results[row.i, paste0("EDx", ED)] <- 10^logdose.ED } } } diff --git a/R/drplot.R b/R/drplot.R index 1e6b2cd..11c4d42 100644 --- a/R/drplot.R +++ b/R/drplot.R @@ -1,12 +1,12 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("dose", "Substance", "mtype")) -drplot <- function(drresults, data, +drplot <- function(drresults, data, dtype = "std", alpha = 0.95, ctype = "none", path = "./", fileprefix = "drplot", overlay = FALSE, xlim = c("auto","auto"), ylim = c("auto","auto"), xlab = paste("Decadic Logarithm of the dose in ", unit), ylab = "Normalized response", axes = TRUE, frame.plot = TRUE, - postscript = FALSE, pdf = FALSE, png = FALSE, + postscript = FALSE, pdf = FALSE, png = FALSE, bw = TRUE, pointsize = 12, colors = 1:8, ltys = 1:8, pchs = "auto", @@ -24,7 +24,7 @@ drplot <- function(drresults, data, if(is.data.frame(data)) { # Get rid of pseudo substance names of controls nonzerodata <- subset(data,dose!=0) - nonzerodata$substance <- factor(nonzerodata$substance) + nonzerodata$substance <- factor(nonzerodata$substance) zerodata <- subset(data,dose==0) nc <- length(zerodata$dose) # Number of control points if (nc > 0) { @@ -44,12 +44,12 @@ drplot <- function(drresults, data, hr <- max(nonzerodata$response) if (ctype == "std") hr <- max(hr,1 + sdc) if (ctype == "conf") hr <- max(hr,1 + controlconf) - dsubstances <- levels(nonzerodata$substance) + dsubstances <- levels(nonzerodata$substance) } else { lld <- min(drresults[["logED50"]],na.rm=TRUE) - 2 lhd <- max(drresults[["logED50"]],na.rm=TRUE) + 2 if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) { - hr <- 1.8 + hr <- 1.8 } else { hr <- 1.0 } @@ -70,20 +70,20 @@ drplot <- function(drresults, data, postscript(file=filename, paper="special",width=7,height=7,horizontal=FALSE, pointsize=pointsize) message("Created File: ",filename,"\n") - } + } if (pdf) { filename = paste(path,fileprefix,".pdf",sep="") pdf(file=filename, paper="special",width=7,height=7,horizontal=FALSE, pointsize=pointsize) message("Created File: ",filename,"\n") - } + } if (png) { filename = paste(path,fileprefix,".png",sep="") png(filename=filename, width=500, height=500, pointsize=pointsize) message("Created File: ",filename,"\n") } - + plot(0,type="n", xlim = xlim, ylim = ylim, @@ -97,7 +97,7 @@ drplot <- function(drresults, data, if (!postscript && !png && !pdf && length(dsubstances) > 1) { op <- par(ask=TRUE) on.exit(par(op)) - } + } } # nl is the overall number of fits to draw by different line types nl <- 0 @@ -123,20 +123,20 @@ drplot <- function(drresults, data, postscript(file=filename, paper="special",width=7,height=7,horizontal=FALSE,pointsize=pointsize) message("Created File: ",filename,"\n") - } + } if (pdf) { filename = paste(path,fileprefix,sub(" ","_",i),".pdf",sep="") pdf(file=filename, paper="special",width=7,height=7,horizontal=FALSE,pointsize=pointsize) message("Created File: ",filename,"\n") - } + } if (png) { filename = paste(path,fileprefix,sub(" ","_",i),".png",sep="") png(filename=filename, width=500, height=500, pointsize=pointsize) message("Created File: ",filename,"\n") } - + plot(0,type="n", xlim = xlim, ylim = ylim, @@ -147,7 +147,7 @@ drplot <- function(drresults, data, } if (!overlay) legend(lpos, i, lty = 1, col = color, pch = pch, inset=0.05) tmp$dosefactor <- factor(tmp$dose) # necessary because the old - # factor has all levels, not + # factor has all levels, not # only the ones tested with # this substance @@ -157,7 +157,7 @@ drplot <- function(drresults, data, abline(h = 1 + sdc, lty = 2) } if (ctype == "conf") { - abline(h = 1 - controlconf, lty = 2) + abline(h = 1 - controlconf, lty = 2) abline(h = 1 + controlconf, lty = 2) } @@ -179,7 +179,7 @@ drplot <- function(drresults, data, } if (dtype == "conf") { - confidencedeltas <- qt((1 + alpha)/2, lengths - 1) * sqrt(vars) + confidencedeltas <- qt((1 + alpha)/2, lengths - 1) * sqrt(vars) tops <- means + confidencedeltas bottoms <- means - confidencedeltas } @@ -207,13 +207,13 @@ drplot <- function(drresults, data, if (drresults[[fit.row, "mtype"]] %in% c("probit", "logit", "weibull", - "linlogit")) + "linlogit")) { m <- attr(drresults, "models")[[as.numeric(fit.row)]] of <- function(x) { predict(m, data.frame(dose = 10^x)) - } - plot(of, lld - 0.5, lhd + 2, + } + plot(of, lld - 0.5, lhd + 2, add = TRUE, col = color, lty = lty) } } diff --git a/data/XY.rda b/data/XY.rda index 4b0d426..a4f6242 100644 Binary files a/data/XY.rda and b/data/XY.rda differ diff --git a/inst/doc/drfit-Rnews.pdf b/inst/doc/drfit-Rnews.pdf index 65ad25f..0eb0de4 100644 Binary files a/inst/doc/drfit-Rnews.pdf and b/inst/doc/drfit-Rnews.pdf differ -- cgit v1.2.1