From 51d591316daab73ca0f7ec3cf0b0c5925005d29b Mon Sep 17 00:00:00 2001 From: ranke Date: Fri, 25 Feb 2005 21:43:31 +0000 Subject: 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 --- R/drfit.R | 177 ++++++++++++++++++++++++++++++++------------------------------ 1 file changed, 91 insertions(+), 86 deletions(-) (limited to 'R') 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) -- cgit v1.2.1