From fdd212ca853af2b41736ba6e1cd539feb5f145de Mon Sep 17 00:00:00 2001
From: ranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>
Date: Mon, 2 Jan 2006 16:41:43 +0000
Subject: - Added a little documentation of the sample data in inst/doc. -
 Added the possibility to show standard deviation or confidence   interval for
 the mean of the control values, if present in   the dose-response data. Also,
 control data (with dose=0)   are not listed in the results data frame of
 drfit any more.

git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@54 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc
---
 R/drfit.R | 46 ++++++++++++++++++++++++++++++++++++----------
 1 file changed, 36 insertions(+), 10 deletions(-)

(limited to 'R')

diff --git a/R/drfit.R b/R/drfit.R
index 9184dc9..ac3ad3e 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -287,7 +287,8 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
     return(results)
 }
 
-drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
+drplot <- function(drresults, data, 
+        dtype = "std", alpha = 0.95, ctype = "none",
         path = "./", fileprefix = "drplot", overlay = FALSE,
         postscript = FALSE, png = FALSE, bw = TRUE,
         pointsize = 12,
@@ -301,16 +302,30 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
         unit <- "different units"
     }
 
-    # Get the plot limits on the x-axis (log of the dose)
+    # Get the plot limits on the x-axis (log of the dose) and y axis
     if(is.data.frame(data)) {
-        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))
+        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
+        sdc <- sd(zerodata$response)    # Standard deviation of control responses
+        controlconf <- sdc * qt((1 + alpha)/2, nc - 1) / sqrt(nc)
+        if (nc > 0) {
+            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"
+            }
         }
-        lhd <- log10(max(data$dose))
-        hr <- max(data$response)
-        dsubstances <- levels(data$substance)    
+        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
@@ -349,7 +364,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
 
     # Plot the data either as raw data or as error bars
     if(is.data.frame(data)) {
-        splitted <- split(data,data$substance)
+        splitted <- split(nonzerodata,nonzerodata$substance)
         # n is the index for the dose-response curves
         n <- 0
         if (bw) colors <- rep("black",length(dsubstances))
@@ -389,6 +404,17 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
                                                     # 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") {
-- 
cgit v1.2.1