diff options
| -rw-r--r-- | ChangeLog | 8 | ||||
| -rw-r--r-- | DESCRIPTION | 24 | ||||
| -rw-r--r-- | NAMESPACE | 3 | ||||
| -rw-r--r-- | R/checkexperiment.R | 28 | ||||
| -rw-r--r-- | R/checksubstance.R | 52 | ||||
| -rw-r--r-- | R/drcfit.R | 68 | ||||
| -rw-r--r-- | R/drdata.R | 2 | ||||
| -rw-r--r-- | R/drfit.R | 47 | ||||
| -rw-r--r-- | R/drplot.R | 36 | ||||
| -rw-r--r-- | data/XY.rda | bin | 1602 -> 630 bytes | |||
| -rw-r--r-- | inst/doc/drfit-Rnews.pdf | bin | 229223 -> 169560 bytes |
11 files changed, 140 insertions, 128 deletions
@@ -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 <jranke@uni-bremen.de> +Repository: CRAN +Date/Publication: 2016-09-02 21:28:14 @@ -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 = " ") @@ -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"] } } } @@ -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 '", @@ -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 } } } @@ -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 Binary files differindex 4b0d426..a4f6242 100644 --- a/data/XY.rda +++ b/data/XY.rda diff --git a/inst/doc/drfit-Rnews.pdf b/inst/doc/drfit-Rnews.pdf Binary files differindex 65ad25f..0eb0de4 100644 --- a/inst/doc/drfit-Rnews.pdf +++ b/inst/doc/drfit-Rnews.pdf |