diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/drfit.R | 240 |
1 files changed, 100 insertions, 140 deletions
@@ -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")
|