rewritten drplot function, documentation not up to date

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

2 files changed, 110 insertions, 152 deletions

diff --git a/R/drfit.R b/R/drfit.R
index fde2243..b328171 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -31,20 +31,17 @@ drdata <- function(substances, experimentator = "%", db = "cytotox",
     return(data)

 }

     

+linearlogisf <- function(x,k,f,mu,b)

+{

+    k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))

+}

+

 drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,

         lognorm = TRUE, logis = FALSE,

         linearlogis = FALSE, b0 = 2, f0 = 0)

 {

     substances <- levels(data$substance)

     unit <- levels(as.factor(data$unit))

-    logisf <- function(x,x0,b,k=1)

-    {

-        k / (1 + (x/x0)^b)

-    }

-    linearlogisf <- function(x,k,f,mu,b)

-    {

-        k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))

-    }

 

     ri <- rix <- 0                        # ri is the index over the result rows

                                           # rix is used later to check if any

@@ -114,8 +111,6 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                 }

 

                 if (logis) {

-                # Instead of plogis(), the function logisf() defined above

-                # could be used for regression against dose, not log10(dose)

                     m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope),

                             data=tmp,

                             start=list(logEC50=startlogEC50[[i]],slope=1)))

@@ -215,30 +210,33 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
     return(results)

 }

 

-drplot <- function(drresults, data = FALSE, dtype = "std", alpha = 0.95,

+drplot <- function(drresults, data, dtype = "std", alpha = 0.95,

         path = "./", fileprefix = "drplot", overlay = FALSE,

-        postscript = FALSE, 

-        color = TRUE,

-        datacolors = 1:8, fitcolors = "default")

+        postscript = FALSE, png = FALSE, bw = TRUE,

+        colors = 1:8)

 {

-    # Prepare plots

-    devices <- 1

-    if (postscript && !overlay) psdevices <- vector()

-    if (!postscript && !overlay) x11devices <- vector()

-

-    unit <- levels(as.factor(drresults$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)

     if(is.data.frame(data))

     {

-        lld <- log10(min(data$dose))

+        if (min(data$dose == 0)) {

+            cat("At least one of the dose levels is 0 - this is not a valid dose.")

+        } else {

+            lld <- log10(min(data$dose))

+        }

         lhd <- log10(max(data$dose))

         hr <- max(data$response)

         dsubstances <- levels(data$substance)    

     } else {

         lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2

         lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2

-        if (length(subset(drresults,mtype=="linearlogis")$substance) != 0) {

+        if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) {

             hr <- 1.8 

         } else {

             hr <- 1.0

@@ -248,43 +246,55 @@ drplot <- function(drresults, data = FALSE, dtype = "std", alpha = 0.95,
     # Prepare overlay plot if requested

     if (overlay)

     {

-        devices <- devices + 1

         if (postscript) {

             filename = paste(path,fileprefix,".eps",sep="")

             postscript(file=filename,

                     paper="special",width=7,height=7,horizontal=FALSE,pointsize=12) 

             cat("Created File: ",filename,"\n")

-        } else {

+        } 

+        if (png) {

+            filename = paste(path,fileprefix,".png",sep="")

+            png(filename=filename,

+                width=500, height=500,pointsize=12)

+            cat("Created File: ",filename,"\n")

+        }

+        if (!postscript && !png) {

             x11(width=7,height=7,pointsize=12)

         }

             

         plot(0,type="n",

-            xlim=c(lld - 0.5, lhd + 2),

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

-    {

+    if(is.data.frame(data)) {

         splitted <- split(data,data$substance)

         n <- 0

-        for (i in dsubstances)

-        {

+        if (bw) colors <- rep("black",length(dsubstances))

+        for (i in dsubstances) {

+            n <- n + 1

+            tmp <- splitted[[i]]

+            color <- colors[[n]]

             # Prepare the single graphs if an overlay is not requested

             if (!overlay)

             {

-                devices <- devices + 1

                 if (postscript) {

                     filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".eps",sep="")

                     postscript(file=filename,

                             paper="special",width=7,height=7,horizontal=FALSE,pointsize=12)

-                    psdevices[[i]] <- devices

                     cat("Created File: ",filename,"\n")

-                } else {

+                } 

+                if (png) {

+                    filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", i)),".png",sep="")

+                    png(filename=filename,

+                        width=500, height=500,pointsize=12)

+                    cat("Created File: ",filename,"\n")

+                }

+                if (!postscript && !png) {

                     x11(width=7,height=7,pointsize=12)

-                    x11devices[[i]] <- devices

                 }

                     

                 plot(0,type="n",

@@ -293,122 +303,72 @@ drplot <- function(drresults, data = FALSE, dtype = "std", alpha = 0.95,
                     xlab=paste("Decadic Logarithm of the dose in ", unit),    

                     ylab="Normalized response")

             }

-            if (color == FALSE) datacolors <- rep("black",length(dsubstances))

-            n <- n + 1

-            if (!overlay) legend(lhd - 1, hr + 0.1, i,lty = 1, col = datacolors[[n]])

-            tmp <- splitted[[i]]

+            if (!overlay) legend(lhd - 1, hr + 0.1, i,lty = 1, col = color)

             tmp$dosefactor <- factor(tmp$dose)  # necessary because the old

                                                 # factor has all levels, not 

                                                 # only the ones tested with

                                                 # this substance

-            if (dtype == "raw") {

-                points(log10(tmp$dose),tmp$response,col=datacolors[[n]])

-            } 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=datacolors[[n]])

-                points(x,means,col=datacolors[[n]])

-                smidge <- 0.05

-                segments(x - smidge,bottoms,x + smidge,bottoms,col=datacolors[[n]])

-                segments(x - smidge,tops,x + smidge,tops,col=datacolors[[n]])

+            # 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 fitted dose response models from drresults

-    fits <- subset(drresults,!is.na(logEC50))

-    nf <-  length(fits$Substance)  # number of fits to plot

-    if (fitcolors[[1]] == "default")

-    {

-        if (nf <= 8)

-        {

-            defaultfitcolors <- palette()

-        } else

-        {

-            defaultfitcolors <- rainbow(nf)

-        }

-    }

-    legendcolors <- vector()

-    for (i in 1:nf)

-    {

-        s <- as.character(fits[i,"Substance"]) # The substance name to display

-        if (!overlay && !is.data.frame(data))

-        {

-            devices <- devices + 1

-            if (postscript) {

-                filename = paste(path,fileprefix,sub(" ","_",gsub("([\(\) ])", "", s)),".eps",sep="")

-                postscript(file=filename,

-                        paper="special",width=7,height=7,horizontal=FALSE,pointsize=12) 

-                psdevices[[s]] <- devices

-                cat("Created File: ",filename,"\n")

-            } else {

-                x11(width=7,height=7,pointsize=12)

-                x11devices[[s]] <- devices

-            }

-                

-            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 (postscript && !overlay) {

-            dev.set(psdevices[[s]]) }

-        if (!postscript && !overlay) {

-            dev.set(x11devices[[s]]) }

-

-        if (color == FALSE) {

-            fitcolor <- "black" 

-        } else {

-            if (fitcolors[[1]] == "default")

-            {

-                fitcolor <- defaultfitcolors[[i]]

-            } else {

-                fitcolor <- fitcolors[[i]]

+            # 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)

+                {

+                    logEC50 <- fits[j,"logEC50"]

+                    mtype <- as.character(fits[j, "mtype"])

+                    if (mtype == "lognorm") {

+                        slope <- fits[j,"slope"]

+                        plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)

+                    }

+                    if (mtype == "logis") {

+                        slope <- fits[j,"slope"]

+                        plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)

+                    }

+                    if (mtype == "linearlogis") {

+                        plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]),

+                            lld - 0.5, lhd + 2,

+                            add=TRUE,col=color)

+                    }

+                }

             }

-        }

-        if (!overlay) legend(lhd - 1, hr + 0.1, s,lty = 1, col = fitcolor)

-        legendcolors[[i]] <- fitcolor

-        logEC50 <- fits[i,"logEC50"]

-        mtype <- as.character(fits[i, "mtype"])

-        if (mtype == "lognorm")

-        {

-            slope <- fits[i,"slope"]

-            plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=fitcolor)

-        }

-        if (mtype == "logis")

-        {

-            slope <- fits[i,"slope"]

-            plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=fitcolor)

-        }

-    }

-    if (overlay) {

-        legend(lhd - 1, hr + 0.1, as.vector(fits$Substance), lty = 1, col = legendcolors)

-    }

-    if (devices > 1 && postscript)

-    {

-        for (i in 2:devices) {

-            dev.off(i)

+        if (!overlay && (postscript || png)) dev.off()

         }

     }

+    if (overlay) legend(lhd - 1, hr + 0.1, dsubstances,lty = 1, col = colors)

+    if (overlay && (postscript || png)) dev.off()

 }

 

 checkplate <- function(plate,db="cytotox")

diff --git a/man/drplot.Rd b/man/drplot.Rd
index 3708d4b..418de15 100644
--- a/man/drplot.Rd
+++ b/man/drplot.Rd
@@ -14,20 +14,18 @@
     A data frame as returned from \code{\link{drfit}}.
     }
   \item{data}{ 
-    A data frame as returned from \code{\link{drdata}}.  If data is to be
-    plotted, the data frame has to contain at least a factor called
-    "substance", a vector called "unit" containing the unit used for the dose,
-    a column "response" with the response values of the test system normalized
-    between 0 and 1, a column "dose" with the numeric dose values and a factor
-    called "dosefactor" containing the dose as a factor. If plotting of the data is 
-    not required, data can be set to FALSE.
+    A data frame as returned from \code{\link{drdata}}. The data frame has to
+    contain at least a factor called "substance", a vector called "unit"
+    containing the unit used for the dose, a column "response" with the
+    response values of the test system normalized between 0 and 1, a column
+    "dose" with the numeric dose values and a factor called "dosefactor"
+    containing the dose as a factor. 
     }
   \item{dtype}{
-    A string describing if the raw data should be plotted ("raw"), or 
-    an error bar should be constructed from the standard deviations of the
-    responses at each dose level ("std", default value) or from the confidence
-    intervals ("conf"). Of course, dtype only makes a difference, if a valid data
-    object has been referenced.
+    A string describing if the raw data should be plotted ("raw"), or an error
+    bar should be constructed from the standard deviations of the responses at
+    each dose level ("std", default value) or from the confidence intervals
+    ("conf"). If you don't want to see the data, set it to "none".
     }
   \item{alpha}{
     The confidence level, defaulting to 0.95, only used if dtype "conf" has been