Corrected the split of the drift.R file.

Added example data IM1xIPC81.




git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@61 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc

1 files changed, 248 insertions, 34 deletions

diff --git a/R/drfit.R b/R/drfit.R
index 9da7d96..ff03d90 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,39 +1,253 @@
-drdata <- function(substances, experimentator = "%", db = "cytotox",
-    celltype="IPC-81",enzymetype="AChE",
-    organism="Vibrio fischeri",endpoint="Luminescence",whereClause="1",
-    ok="'ok','no fit'")
+drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
+        probit = TRUE, logit = FALSE, weibull = FALSE,
+        linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
+        s0 = 0.5, b0 = 2, f0 = 0)
 {
-    library(RODBC) 
-    channel <- odbcConnect(db,uid="cytotox",pwd="cytotox",case="tolower")
-    slist <- paste(substances,collapse="','")
-    if (db == "cytotox") {
-        responsetype <- "viability"
-        testtype <- "celltype"
-        type <- celltype
-    } else {
-        if (db == "enzymes") {
-            responsetype <- "activity"
-            testtype <- "enzyme"
-            type <- enzymetype
+    if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data with 
+                                                            # ok set to "no fit"
+    substances <- levels(data$substance)
+
+    ri <- rix <- 0                        # ri is the index over the result rows
+                                          # rix is used later to check if any
+                                          # model result was appended
+    rsubstance <- array()                 # the substance names in the results
+    rndl <- vector()                      # number of dose levels
+    rn <- vector()                        # mean number of replicates 
+                                          # in each dose level
+    runit <- vector()                     # vector of units for each result row
+    rlhd <- rlld <- vector()              # highest and lowest doses tested
+    mtype <- array()                      # the modeltypes
+    sigma <- array()                      # the standard deviation of the residuals
+    logED50 <- vector()
+    stderrlogED50 <- vector()
+    a <- b <- c <- vector()
+
+    splitted <- split(data,data$substance)
+    for (i in substances) {
+        tmp <- splitted[[i]]
+        fit <- FALSE
+        if (length(tmp) != 0) {
+            unit <- levels(as.factor(as.vector(tmp$unit)))
+            cat("\n",i,": Fitting data...\n",sep="")
         } else {
-            responsetype <- "response"
-            testtype <- "organism"
-            type <- organism
+            unit <- ""
+            cat("\n",i,": No data\n",sep="")
         }
+        if (length(unit) > 1) {
+            cat("More than one unit for substance ",i,", halting\n\n",sep="")
+            break
+        }
+        if (length(tmp$response) == 0) {
+            nodata = TRUE
+        } else {
+            nodata = FALSE
+        }
+        rix <- ri
+        if (nodata) {
+            n <- ndl <- 0
+        } else {
+            ndl <- length(levels(factor(tmp$dose)))
+            n <- round(length(tmp$response)/ndl)
+            if (is.na(startlogED50[i])){
+                w <- 1/abs(tmp$response - 0.3)
+                startlogED50[[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)
+            responseatlowestdose <- mean(subset(tmp,dose==lowestdose)$response)
+            if (responseathighestdose < 0.5) {
+                inactive <- FALSE
+                if (responseatlowestdose < 0.5) {
+                    active <- TRUE
+                } else {
+                    active <- FALSE
+                    if (linlogit && 
+                        length(subset(linlogitWrong,linlogitWrong == i))==0 &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b),
+                                data=tmp,
+                                start=list(f=f0,logED50=startlogED50[[i]],b=b0)))
+                        if (!inherits(m, "try-error")) {
+                            fit <- TRUE
+                            ri <- ri + 1
+                            s <- summary(m)
+                            sigma[[ri]] <- s$sigma
+                            rsubstance[[ri]] <- i
+                            rndl[[ri]] <- ndl
+                            rn[[ri]] <- n
+                            runit[[ri]] <- unit
+                            rlld[[ri]] <- log10(lowestdose)
+                            rlhd[[ri]] <- log10(highestdose)
+                            logED50[[ri]] <- coef(m)[["logED50"]]
+                            if (logED50[[ri]] > rlhd[[ri]]) {
+                                mtype[[ri]] <- "no fit"
+                                logED50[[ri]] <- NA
+                                stderrlogED50[[ri]] <- NA
+                                a[[ri]] <- NA
+                                b[[ri]] <- NA
+                                c[[ri]] <- NA
+                            } else {
+                                mtype[[ri]] <- "linlogit"
+                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                a[[ri]] <- coef(m)[["logED50"]]
+                                b[[ri]] <- coef(m)[["b"]]
+                                c[[ri]] <- coef(m)[["f"]]
+                            }
+                        }
+                    }
+
+                    if (probit &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale),
+                                    data=tmp,
+                                    start=list(logED50=startlogED50[[i]],scale=1)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- coef(m)[["logED50"]]
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "probit"
+                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                    a[[ri]] <- coef(m)[["logED50"]]
+                                    b[[ri]] <- coef(m)[["scale"]]
+                                }
+                            }
+                        }
+                    }
+
+                    if (logit &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale),
+                                data=tmp,
+                                start=list(logED50=startlogED50[[i]],scale=1)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- a[[ri]] <- coef(m)[["logED50"]]
+                                b[[ri]] <- coef(m)[["scale"]]
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "logit"
+                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                }
+                            }
+                        }
+                    }
+
+                    if (weibull &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
+                                data=tmp,
+                                start=list(location=startlogED50[[i]],shape=s0)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                a[[ri]] <- coef(m)[["location"]]
+                                b[[ri]] <- coef(m)[["shape"]]
+                                sqrdev <- function(logdose) {
+                                    (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 
+                                }
+                                logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "weibull"
+                                    stderrlogED50[[ri]] <- NA
+                                }
+                            }
+                        }
+                    }
+
+                }
+
+            } else {
+                inactive <- TRUE
+            }
+        }
+        if (ri == rix) {          # if no entry was appended for this substance
+            ri <- ri + 1
+            rsubstance[[ri]] <- i
+            rndl[[ri]] <- ndl
+            rn[[ri]] <- n
+            if (nodata) {
+                rlld[[ri]] <- rlhd[[i]] <- NA
+                mtype[[ri]] <- "no data"
+                runit[[ri]] <- NA
+            } else {
+                rlld[[ri]] <- log10(lowestdose)
+                rlhd[[i]] <- log10(highestdose)
+                runit[[ri]] <- unit
+                if (inactive) {
+                    mtype[[ri]] <- "inactive"
+                } else {
+                    if (active) {
+                        mtype[[ri]] <- "active"
+                    } else {
+                        mtype[[ri]] <- "no fit"
+                    }
+                }
+            }
+            sigma[[ri]] <- NA
+            logED50[[ri]] <- NA
+            stderrlogED50[[ri]] <- NA
+            a[[ri]] <- NA
+            b[[ri]] <- NA
+            c[[ri]] <- NA
+        }
+    }
+    results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b)
+    names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b")
+    if (linlogit) {
+        results$c <- c
     }
-        
-    query <- paste("SELECT conc,",responsetype,",unit,experimentator,substance,",testtype,
-        ",ok FROM ", db, " WHERE substance IN ('",
-        slist,"') AND experimentator LIKE '",
-        experimentator,"' AND ",testtype," LIKE '",
-        type,"' AND ",
-        whereClause," AND ok in (",
-        ok,")",sep="")
-    if (db == "ecotox") query <- paste(query," AND type LIKE '",endpoint,"'",sep="")
-    data <- sqlQuery(channel,query)
-    odbcClose(channel)
-    names(data)[[1]] <- "dose"
-    names(data)[[2]] <- "response"
-    data$substance <- factor(data$substance,levels=substances)
-    return(data)
+    rownames(results) <- 1:ri
+    return(results)
 }