diff options
Diffstat (limited to 'R/drplot.R')
| -rw-r--r-- | R/drplot.R | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/R/drplot.R b/R/drplot.R new file mode 100644 index 0000000..baec69d --- /dev/null +++ b/R/drplot.R @@ -0,0 +1,220 @@ +drplot <- function(drresults, data, + dtype = "std", alpha = 0.95, ctype = "none", + path = "./", fileprefix = "drplot", overlay = FALSE, + postscript = FALSE, pdf = FALSE, png = FALSE, + bw = TRUE, + pointsize = 12, + colors = 1:8, devoff=T, lpos="topright") +{ + # Check if all data have the same unit + unitlevels <- levels(as.factor(drresults$unit)) + if (length(unitlevels) == 1) { + unit <- unitlevels + } else { + unit <- "different units" + } + + # Get the plot limits on the x-axis (log of the dose) and y axis + if(is.data.frame(data)) { + nonzerodata <- subset(data,dose!=0) + nonzerodata$substance <- factor(nonzerodata$substance) # Get rid of pseudo substance names of controls + zerodata <- subset(data,dose==0) + nc <- length(zerodata$dose) # Number of control points + if (nc > 0) { + sdc <- sd(zerodata$response) # Standard deviation of control responses + controlconf <- sdc * qt((1 + alpha)/2, nc - 1) / sqrt(nc) + cat("There are ",nc,"data points with dose 0 (control values)\n") + cat("with a standard deviation of",sdc,"\n") + cat("and a confidence interval of",controlconf,"\n") + if (nc < 3) { + cat("\nThere are less than 3 control points, therefore their scatter\n") + cat("will not be displayed\n") + ctype = "none" + } + } else { + if (ctype != "none") { + stop("There are no controls in the dataset, and therefore ", + "their scatter cannot be displayed\n") + } + } + lld <- log10(min(nonzerodata$dose)) + lhd <- log10(max(nonzerodata$dose)) + hr <- max(nonzerodata$response) + if (ctype == "std") hr <- max(hr,1 + sdc) + if (ctype == "conf") hr <- max(hr,1 + controlconf) + 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 + } else { + hr <- 1.0 + } + } + + # Prepare overlay plot if requested + if (overlay) + { + if (postscript) { + filename = paste(path,fileprefix,".eps",sep="") + postscript(file=filename, + paper="special",width=7,height=7,horizontal=FALSE, pointsize=pointsize) + cat("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) + cat("Created File: ",filename,"\n") + } + if (png) { + filename = paste(path,fileprefix,".png",sep="") + png(filename=filename, + width=500, height=500, pointsize=pointsize) + cat("Created File: ",filename,"\n") + } + if (!postscript && !png && !pdf) { + get(getOption("device"))(width=7,height=7) + } + + plot(0,type="n", + xlim=c(lld - 0.5, lhd + 1), + ylim= c(-0.1, hr + 0.2), + xlab=paste("Decadic Logarithm of the dose in ", unit), + ylab="Normalized response") + } + + # Plot the data either as raw data or as error bars + if(is.data.frame(data)) { + splitted <- split(nonzerodata,nonzerodata$substance) + # n is the index for the dose-response curves + n <- 0 + if (bw) colors <- rep("black",length(dsubstances)) + # Loop over the substances in the data + for (i in dsubstances) { + n <- n + 1 + tmp <- splitted[[i]] + if (length(tmp$response) != 0) { + color <- colors[[n]] + # Prepare the single graphs if an overlay is not requested + if (!overlay) + { + if (postscript) { + filename = paste(path,fileprefix,sub(" ","_",i),".eps",sep="") + postscript(file=filename, + paper="special",width=7,height=7,horizontal=FALSE,pointsize=pointsize) + cat("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) + cat("Created File: ",filename,"\n") + } + if (png) { + filename = paste(path,fileprefix,sub(" ","_",i),".png",sep="") + png(filename=filename, + width=500, height=500, pointsize=pointsize) + cat("Created File: ",filename,"\n") + } + if (!postscript && !png && !pdf) { + get(getOption("device"))(width=7,height=7) + } + + plot(0,type="n", + xlim=c(lld - 0.5, lhd + 2), + ylim= c(-0.1, hr + 0.2), + xlab=paste("Decadic Logarithm of the dose in ", unit), + ylab="Normalized response") + } + if (!overlay) legend(lpos, i,lty = 1, col = color, inset=0.05) + tmp$dosefactor <- factor(tmp$dose) # necessary because the old + # factor has all levels, not + # only the ones tested with + # this substance + + # Plot the control lines, if requested + if (ctype == "std") { + abline(h = 1 - sdc, lty = 2) + abline(h = 1 + sdc, lty = 2) + } + if (ctype == "conf") { + abline(h = 1 - controlconf, lty = 2) + abline(h = 1 + controlconf, lty = 2) + } + + # Plot the data, if requested + if (dtype != "none") { + if (dtype == "raw") { + points(log10(tmp$dose),tmp$response,col=color) + } else { + splitresponses <- split(tmp$response,tmp$dosefactor) + means <- sapply(splitresponses,mean) + lengths <- sapply(splitresponses,length) + vars <- sapply(splitresponses,var) + standarddeviations <- sqrt(vars) + } + if (dtype == "std") + { + tops <- means + standarddeviations + bottoms <- means - standarddeviations + } + if (dtype == "conf") + { + confidencedeltas <- qt((1 + alpha)/2, lengths - 1) * sqrt(vars) + tops <- means + confidencedeltas + bottoms <- means - confidencedeltas + } + if (dtype != "raw") + { + x <- log10(as.numeric(levels(tmp$dosefactor))) + segments(x,bottoms,x,tops,col=color) + points(x,means,col=color) + smidge <- 0.05 + segments(x - smidge,bottoms,x + smidge,bottoms,col=color) + segments(x - smidge,tops,x + smidge,tops,col=color) + } + } + + # Plot the fits, if there are any + fits <- subset(drresults,Substance == i) + nf <- length(fits$Substance) # number of fits to plot + if (nf > 0) { + for (j in 1:nf) + { + logED50 <- fits[j,"logED50"] + mtype <- as.character(fits[j, "mtype"]) + if (mtype == "probit") { + scale <- fits[j,"b"] + plot(function(x) pnorm(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) + } + if (mtype == "logit") { + scale <- fits[j,"b"] + plot(function(x) plogis(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color) + } + if (mtype == "weibull") { + location <- fits[j,"a"] + shape <- fits[j,"b"] + plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color) + } + if (mtype == "linlogit") { + plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logED50"],fits[j,"b"]), + lld - 0.5, lhd + 2, + add=TRUE,col=color) + } + } + } + if (!overlay && (postscript || png || pdf)) dev.off() + } else { + cat("No data for ",i,"\n") + } + } + } + if (overlay) legend(lpos, dsubstances,lty = 1, col = colors, inset=0.05) + if (overlay && (postscript || png || pdf)) { + if (devoff) { + dev.off() + } + } +} |