aboutsummaryrefslogtreecommitdiff
path: root/R/drplot.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/drplot.R')
-rw-r--r--R/drplot.R220
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()
+ }
+ }
+}

Contact - Imprint