Some bugfixes and the addition of the linearlogisWrong argument

to the drfit function.



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

5 files changed, 118 insertions, 111 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 15ffa95..69f3d1d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: drfit
-Version: 0.03-13
-Date: 2005-02-24
+Version: 0.03-16
+Date: 2005-02-25
 Title: Dose-response data evaluation
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
diff --git a/R/drfit.R b/R/drfit.R
index b328171..9516180 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -38,7 +38,8 @@ linearlogisf <- function(x,k,f,mu,b)
 

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

         lognorm = TRUE, logis = FALSE,

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

+        linearlogis = FALSE, linearlogisWrong = NA, 

+        b0 = 2, f0 = 0)

 {

     substances <- levels(data$substance)

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

@@ -67,6 +68,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
         } else {

             nodata = FALSE

         }

+        rix <- ri

         if (!nodata) {

             if (is.na(startlogEC50[i])){

                 w <- 1/abs(tmp$response - 0.3)

@@ -77,11 +79,11 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
             lhd <- log10(highestdose)

             lld <- log10(lowestdose)

             responseathighestdose <- mean(subset(tmp,dose==highestdose)$response)

-            rix <- ri

             if (responseathighestdose < 0.5) {

                 inactive <- FALSE

 

-                if (linearlogis) {

+                if (linearlogis && 

+                    length(subset(linearlogisWrong,linearlogisWrong == i))==0) {

                     m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b),

                             data=tmp,

                             start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))

@@ -223,8 +225,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
     }

 

     # Get the plot limits on the x-axis (log of the dose)

-    if(is.data.frame(data))

-    {

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

@@ -277,94 +278,98 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
         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)

-            {

-                if (postscript) {

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

-                    postscript(file=filename,

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

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

-                } 

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

-                }

-                    

-                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(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

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

+            if (length(tmp$response) != 0) {

+                color <- colors[[n]]

+                # Prepare the single graphs if an overlay is not requested

+                if (!overlay)

                 {

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

+                    if (postscript) {

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

+                        postscript(file=filename,

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

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

+                    } 

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

+                    }

+                        

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

                 }

-            }

-

-            # 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 (!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

+                # 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 (mtype == "logis") {

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

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

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

                     }

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

+                }

+

+                # 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 && (postscript || png)) dev.off()

+            } else {

+                cat("No data for ",i,"\n")

             }

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

         }

     }

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

diff --git a/man/drfit.Rd b/man/drfit.Rd
index 30f3fd4..ddcf80b 100644
--- a/man/drfit.Rd
+++ b/man/drfit.Rd
@@ -7,7 +7,7 @@
 }
 \usage{
   drfit(data, startlogEC50 = NA, chooseone = TRUE, lognorm = TRUE, logis = FALSE, 
-    linearlogis = FALSE, b0 = 2, f0 = 0)
+    linearlogis = FALSE, linearlogisWrong = NA, b0 = 2, f0 = 0)
 }
 \arguments{
   \item{data}{
@@ -30,8 +30,12 @@
     A boolean defining if cumulative densitiy curves of logistic distributions
     are fitted to the data. Default is FALSE.} 
   \item{linearlogis}{
-    A boolean defining if the linear-logistic function as defined by van Ewijk and Hoekstra
-    1993 is fitted to the data. Default is FALSE.}
+    A boolean defining if the linear-logistic function
+    \code{\link{linearlogisf}} as defined by van Ewijk and Hoekstra 1993 is
+    fitted to the data. Default is FALSE.}
+  \item{linearlogisWrong}{
+    An optional vector containing the names of the substances for which the
+    linearlogis function produces a wrong fit.}
   \item{chooseone}{
     If TRUE (default), the models are tried in the order linearlogis, logis and lognorm,
     and the first model that produces a valid fit is used. Usually this will be the one 
diff --git a/man/drplot.Rd b/man/drplot.Rd
index 418de15..784b65b 100644
--- a/man/drplot.Rd
+++ b/man/drplot.Rd
@@ -7,7 +7,7 @@
 }
 \usage{
   drplot(drresults, data, dtype, alpha, path, fileprefix, overlay,
-    postscript, color, datacolors, fitcolors)
+    postscript, png, bw, colors)
 }
 \arguments{
   \item{drresults}{
@@ -46,23 +46,21 @@
     will all be put into the last graph in this case.
     }
   \item{postscript}{
-    If TRUE, (a) postscript graph(s) will be created. Otherwise, graphics will be
+    If TRUE, (a) postscript graph(s) will be created. Otherwise, and if 
+    the png argument is also FALSE, graphics will be
     displayed with a screen graphics device.
     }
-  \item{color}{
-    If TRUE, a sensible selection of colors will be attempted. If false, everything
-    will be drawn in black
+  \item{png}{
+    If TRUE, (a) png graph(s) will be created. Otherwise, and if the
+    postscript argument is also FALSE, graphics will be displayed with a
+    screen graphics device.
     }
-  \item{datacolors}{
-    This is a vector of colors, defaulting to 1:8, used for plotting the data.
+  \item{bw}{
+    A boolean deciding if the plots will be black and white or not. Default 
+    is TRUE.
     }
-  \item{fitcolors}{
-    Here you can specify a palette for the colors of the dose-response fits. The 
-    default value is "default", which produces the default palette, if the 
-    number of fits to be plotted is 8 or less. Otherwise, rainbow colors 
-    will be plotted. Unless there is more than one fit per substance to be plotted,
-    or the number of fits is larger than 8, the fitcolors will match the
-    datacolors.
+  \item{colors}{
+    This is a vector of colors, defaulting to 1:8, used for plotting the data.
     }
 }
 \value{
@@ -73,9 +71,9 @@
     
 } 
 \note{
-  Turn off the colors if you don't like them and don't want to fiddle with
-  them. Treatment of legends is quite bad. Be sure all devices are closed
-  (e.g. by calling \code{dev.off()}) before calling \code{drplot} again.
+  Treatment of legends is not really well solved. Be sure all devices are
+  closed (e.g. by calling \code{dev.off()}) before calling \code{drplot}
+  again after a failure.
 }
 \examples{
 data(antifoul)
diff --git a/man/linearlogisf.Rd b/man/linearlogisf.Rd
index d49dc9a..1ebe8be 100644
--- a/man/linearlogisf.Rd
+++ b/man/linearlogisf.Rd
@@ -23,9 +23,9 @@
 \value{
   The response at dose x.
 } 
-\examples{
-  
-}
+\references{
+  van Ewijk, P. H. and Hoekstra, J. A. (1993) \emph{Ecotox Environ Safety}
+  \bold{25} 25-32}
 \author{
   Johannes Ranke 
   \email{jranke@uni-bremen.de}