Split of the drfit.R file to have a separate file for each

function.



git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@60 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc

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()
+        }
+    }
+}