aboutsummaryrefslogtreecommitdiff
path: root/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
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')
-rw-r--r--R/drdata.R39
-rw-r--r--R/drfit.R282
-rw-r--r--R/linlogitf.R254
3 files changed, 287 insertions, 288 deletions
diff --git a/R/drdata.R b/R/drdata.R
new file mode 100644
index 0000000..9da7d96
--- /dev/null
+++ b/R/drdata.R
@@ -0,0 +1,39 @@
+drdata <- function(substances, experimentator = "%", db = "cytotox",
+ celltype="IPC-81",enzymetype="AChE",
+ organism="Vibrio fischeri",endpoint="Luminescence",whereClause="1",
+ ok="'ok','no fit'")
+{
+ 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
+ } else {
+ responsetype <- "response"
+ testtype <- "organism"
+ type <- organism
+ }
+ }
+
+ 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)
+}
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)
}
diff --git a/R/linlogitf.R b/R/linlogitf.R
index 2af8b7e..2d69256 100644
--- a/R/linlogitf.R
+++ b/R/linlogitf.R
@@ -2,257 +2,3 @@ linlogitf <- function(x,k,f,mu,b)
{
k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))
}
-
-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)
-{
- 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 {
- 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
- }
- rownames(results) <- 1:ri
- return(results)
-}

Contact - Imprint