From b8a7ce784ce498130058a9cff999dc52955461fa Mon Sep 17 00:00:00 2001 From: ranke Date: Wed, 23 Feb 2005 18:10:47 +0000 Subject: I addded the checksubstance function, and I improved the drfit function, so it now also works with larger datasets, and if some substances have no data. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@11 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc --- R/drfit.R | 269 +++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 169 insertions(+), 100 deletions(-) diff --git a/R/drfit.R b/R/drfit.R index 62434bc..20f4426 100644 --- a/R/drfit.R +++ b/R/drfit.R @@ -21,7 +21,8 @@ drdata <- function(substances, experimentator = "%", db = "cytotox", return(data) } -drfit <- function(data, startlogEC50 = NA, lognorm = TRUE, logis = FALSE, +drfit <- function(data, startlogEC50 = NA, chooseone=TRUE, + lognorm = TRUE, logis = FALSE, linearlogis = FALSE, b0 = 2, f0 = 0) { substances <- levels(data$substance) @@ -40,6 +41,7 @@ drfit <- function(data, startlogEC50 = NA, lognorm = TRUE, logis = FALSE, rn <- vector() # number of dose-response curves rlhd <- rlld <- vector() # highest and lowest doses tested mtype <- array() # the modeltypes + sigma <- array() # the standard deviation of the residuals logEC50 <- vector() stderrlogEC50 <- vector() slope <- vector() @@ -47,128 +49,143 @@ drfit <- function(data, startlogEC50 = NA, lognorm = TRUE, logis = FALSE, f <- vector() splitted <- split(data,data$substance) - for (i in substances) - { + for (i in substances) { tmp <- splitted[[i]] - n <- round(length(tmp$response)/9) - if (is.na(startlogEC50[i])){ - w <- 1/abs(tmp$response - 0.3) - startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w) + if (length(tmp$response) == 0) { + nodata = TRUE + } else { + nodata = FALSE } - highestdose <- max(tmp$dose) - lowestdose <- min(tmp$dose) - lhd <- log10(highestdose) - lld <- log10(lowestdose) - responseathighestdose <- mean(subset(tmp,dose==highestdose)$response) - rix <- ri # rix is used late to check if any - # model result was appended - if (responseathighestdose < 0.5) { - if (lognorm) - { - m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope), - data=tmp, - start=list(logEC50=startlogEC50[[i]],slope=1))) - if (!inherits(m, "try-error")) - { + if (!nodata) { + n <- round(length(tmp$response)/9) + if (is.na(startlogEC50[i])){ + w <- 1/abs(tmp$response - 0.3) + startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w) + } + highestdose <- max(tmp$dose) + lowestdose <- min(tmp$dose) + lhd <- log10(highestdose) + lld <- log10(lowestdose) + responseathighestdose <- mean(subset(tmp,dose==highestdose)$response) + rix <- ri # ri is the index of result lines + # rix is used later to check if any + # model result was appended + if (responseathighestdose < 0.5) { + inactive <- FALSE + + if (lognorm) { + m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope), + data=tmp, + start=list(logEC50=startlogEC50[[i]],slope=1))) + if (!inherits(m, "try-error")) { + s <- summary(m) ri <- ri + 1 - rsubstance[[ri]] <- i - rn[[ri]] <- n - rlld[[ri]] <- log10(lowestdose) - rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "lognorm" - s <- summary(m) - logEC50[[ri]] <- coef(m)[["logEC50"]] - if (logEC50[[ri]] > rlhd[[ri]]) - { - logEC50[[ri]] <- NA - slope[[ri]] <- NA - stderrlogEC50[[ri]] <- NA - } else - { - slope[[ri]] <- coef(m)[["slope"]] - stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] - } + sigma[[ri]] <- s$sigma + rsubstance[[ri]] <- i + rn[[ri]] <- n + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + mtype[[ri]] <- "lognorm" + logEC50[[ri]] <- coef(m)[["logEC50"]] + b[[ri]] <- NA + f[[ri]] <- NA + if (logEC50[[ri]] > rlhd[[ri]]) { + logEC50[[ri]] <- NA + slope[[ri]] <- NA + stderrlogEC50[[ri]] <- NA + } else { + slope[[ri]] <- coef(m)[["slope"]] + stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + } } - } + } - 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))) - if (!inherits(m, "try-error")) - { - ri <- ri + 1 - rsubstance[[ri]] <- i - rn[[ri]] <- n - rlld[[ri]] <- log10(lowestdose) - rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "logis" - s <- summary(m) - logEC50[[ri]] <- coef(m)[["logEC50"]] - if (logEC50[[ri]] > rlhd[[ri]]) - { - logEC50[[ri]] <- NA + 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))) + if (!inherits(m, "try-error")) { + s <- summary(m) + ri <- ri + 1 + rsubstance[[ri]] <- i + rn[[ri]] <- n + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + mtype[[ri]] <- "logis" + sigma[[ri]] <- s$sigma + logEC50[[ri]] <- coef(m)[["logEC50"]] + b[[ri]] <- NA + f[[ri]] <- NA + if (logEC50[[ri]] > rlhd[[ri]]) { + logEC50[[ri]] <- NA slope[[ri]] <- NA stderrlogEC50[[ri]] <- NA - } else - { - slope[[ri]] <- coef(m)[["slope"]] + } else { + slope[[ri]] <- coef(m)[["slope"]] stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + } } } - } - if (linearlogis) - { - m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b), - data=tmp, - start=list(f=f0,logEC50=startlogEC50[[i]],b=b0))) - if (!inherits(m, "try-error")) - { - ri <- ri + 1 - rsubstance[[ri]] <- i - rn[[ri]] <- n - rlld[[ri]] <- log10(lowestdose) - rlhd[[ri]] <- log10(highestdose) - mtype[[ri]] <- "linearlogis" - s <- summary(m) -#print(s) - logEC50[[ri]] <- coef(m)[["logEC50"]] - if (logEC50[[ri]] > rlhd[[ri]]) - { - logEC50[[ri]] <- NA + if (linearlogis) { + m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b), + data=tmp, + start=list(f=f0,logEC50=startlogEC50[[i]],b=b0))) + if (!inherits(m, "try-error")) { + s <- summary(m) + ri <- ri + 1 + rsubstance[[ri]] <- i + rn[[ri]] <- n + rlld[[ri]] <- log10(lowestdose) + rlhd[[ri]] <- log10(highestdose) + mtype[[ri]] <- "linearlogis" + sigma[[ri]] <- s$sigma + logEC50[[ri]] <- coef(m)[["logEC50"]] + slope[[ri]] <- NA + if (logEC50[[ri]] > rlhd[[ri]]) { + logEC50[[ri]] <- NA stderrlogEC50[[ri]] <- NA b[[ri]] <- NA f[[ri]] <- NA - } else - { - stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] + } else { + stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"] b[[ri]] <- coef(m)[["b"]] f[[ri]] <- coef(m)[["f"]] + } } } + } else { + inactive <- TRUE } - } - if (ri == rix) # if no entry was appended for this substance - { + } + if (ri == rix) { # if no entry was appended for this substance ri <- ri + 1 - rsubstance[[ri]] <- i - rn[[ri]] <- n + rsubstance[[ri]] <- i + rn[[ri]] <- n + if (nodata) { + rlld[[ri]] <- rlhd[[i]] <- NA + mtype[[ri]] <- "no data" + } else { rlld[[ri]] <- log10(lowestdose) rlhd[[i]] <- log10(highestdose) - mtype[[ri]] <- "none" - logEC50[[ri]] <- NA - stderrlogEC50[[ri]] <- NA - slope[[ri]] <- NA - b[[ri]] <- NA - f[[ri]] <- NA + if (inactive) { + mtype[[ri]] <- "inactive" + } else { + mtype[[ri]] <- "no fit" + } + } + sigma[[ri]] <- NA + logEC50[[ri]] <- NA + stderrlogEC50[[ri]] <- NA + slope[[ri]] <- NA + b[[ri]] <- NA + f[[ri]] <- NA } } - results <- data.frame(rsubstance,rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, unit) - names(results) <- c("Substance","n", "lld","lhd","mtype","logEC50","std","unit") + results <- data.frame(rsubstance,rn, rlhd, mtype, logEC50, stderrlogEC50, unit, sigma) + names(results) <- c("Substance","n", "lhd","mtype","logEC50","std","unit","sigma") if (lognorm || logis) { results$slope <- slope } @@ -375,7 +392,8 @@ drplot <- function(drresults, data = FALSE, dtype = "std", alpha = 0.95, } } -checkplate <- function(plate,db="cytotox") { +checkplate <- function(plate,db="cytotox") +{ library(RODBC) channel <- odbcConnect(db,uid=db,pwd=db,case="tolower") @@ -444,3 +462,54 @@ checkplate <- function(plate,db="cytotox") { title(main=paste("Plate ",plate," - ",levels(platedata$experimentator)," - ",levels(platedata$celltype))) } } + +checksubstance <- function(substance,db="cytotox",experimentator="%",celltype="%",whereClause="1",ok="%") +{ + library(RODBC) + channel <- odbcConnect(db,uid=db,pwd=db,case="tolower") + query <- paste("SELECT experimentator,substance,celltype,plate,conc,unit,viability,ok FROM cytotox WHERE substance LIKE '", + substance,"' AND experimentator LIKE '", + experimentator,"' AND celltype LIKE '", + celltype,"' AND ", + whereClause," AND ok LIKE '",ok,"'",sep="") + + data <- sqlQuery(channel,query) + odbcClose(channel) + + data$experimentator <- factor(data$experimentator) + data$substance <- factor(data$substance) + substances <- levels(data$substance) + data$celltype <- factor(data$celltype) + data$plate <- factor(data$plate) + plates <- levels(data$plate) + concentrations <- split(data$conc,data$conc) + concentrations <- as.numeric(names(concentrations)) + data$unit <- factor(data$unit) + data$ok <- factor(data$ok) + + if (length(plates)>6) { + palette(rainbow(length(plates))) + } + + plot(log10(data$conc),data$viability, + xlim=c(-2.5, 4.5), + ylim= c(-0.1, 2), + xlab=paste("Decadic Logarithm of the concentration in ",levels(data$unit)), + ylab="Viability") + + platelist <- split(data,data$plate) + + for (i in 1:length(platelist)) { + points(log10(platelist[[i]]$conc),platelist[[i]]$viability,col=i); + } + + legend(3.5,1.7,plates, pch=1, col=1:length(plates)) + title(main=paste(substance," - ",levels(data$experimentator)," - ",levels(data$celltype))) + + cat("Substanz ",substance,"\n", + "\tExperimentator(s):",levels(data$experimentator),"\n", + "\tCell type(s):\t",levels(data$celltype),"\n", + "\tSubstance(s):\t",levels(data$substance),"\n", + "\tPlate(s):\t",plates,"\n\n") + +} -- cgit v1.2.1