aboutsummaryrefslogtreecommitdiff
path: root/R/drfit.R
diff options
context:
space:
mode:
authorranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2006-04-01 04:41:28 +0000
committerranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2006-04-01 04:41:28 +0000
commit27a255ea7e95c1924f34a5d3aa0bcd39ad902b98 (patch)
treeb3e8ff7813d5abf9a62487a32dc430f66c096594 /R/drfit.R
parentac626cca54de0981ee15430d035eedae41b7ed27 (diff)
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
Diffstat (limited to 'R/drfit.R')
-rw-r--r--R/drfit.R282
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)
}

Contact - Imprint