diff options
| author | ranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc> | 2005-02-23 18:10:47 +0000 |
|---|---|---|
| committer | ranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc> | 2005-02-23 18:10:47 +0000 |
| commit | b8a7ce784ce498130058a9cff999dc52955461fa (patch) | |
| tree | 4c39d0f422fb28a0a509fd4f615c9bf25f7715db /R | |
| parent | 407b6818886836727824864935b012ee177163a6 (diff) | |
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
Diffstat (limited to 'R')
| -rw-r--r-- | R/drfit.R | 269 |
1 files changed, 169 insertions, 100 deletions
@@ -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")
+
+}
|