From abd8cb44371374aea0324a81287c42b809c304f5 Mon Sep 17 00:00:00 2001 From: ranke Date: Fri, 25 Feb 2005 14:57:22 +0000 Subject: 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 --- R/drfit.R | 240 ++++++++++++++++++++++++++------------------------------------ 1 file changed, 100 insertions(+), 140 deletions(-) (limited to 'R') 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") -- cgit v1.2.1