diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/drdata.R | 39 | ||||
| -rw-r--r-- | R/drfit.R | 282 | ||||
| -rw-r--r-- | R/linlogitf.R | 254 |
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) +} @@ -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) -} |