if(getRversion() >= '2.15.1') utils::globalVariables(c("ok", "dose")) drfit <- function(data, startlogED50 = NA, chooseone=TRUE, probit = TRUE, logit = FALSE, weibull = FALSE, linlogit = FALSE, level = 0.95, linlogitWrong = NA, allWrong = NA, ps0 = 1, ls0 = 0.5, ws0 = 0.5, b0 = 2, f0 = 0, showED50 = FALSE, EDx = NULL, EDx.tolerance = 1e-4) { 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() logED50low <- logED50high <- vector() a <- b <- c <- vector() models <- list() # a list containing the dose-response models 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))) message("\n",i,": Fitting data...\n") } else { unit <- "" message("\n",i,": No data\n") } if (length(unit) == 0) { unit <- NA } if (length(unit) > 1) { message("More than one unit for substance ",i,", halting\n\n") 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 <- length(tmp$response) 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, algorithm="port", start=list(f=f0,logED50=startlogED50[[i]],b=b0)), silent = TRUE) if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 models[[ri]] <- m 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 logED50low[[ri]] <- NA logED50high[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } else { mtype[[ri]] <- "linlogit" logED50conf <- try(confint(m,"logED50",level=level), silent = TRUE) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] } else { logED50low[[ri]] <- NA logED50high[[ri]] <- NA } 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, algorithm="port", start=list(logED50=startlogED50[[i]],scale=ps0)), silent = TRUE) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 models[[ri]] <- m 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 logED50low[[ri]] <- NA logED50high[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "probit" logED50conf <- try(confint(m,"logED50",level=level), silent = TRUE) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] } else { logED50low[[ri]] <- NA logED50high[[ri]] <- NA } 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, algorithm="port", start=list(logED50=startlogED50[[i]],scale=ls0)), silent = TRUE) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { fit <- TRUE ri <- ri + 1 models[[ri]] <- m 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 logED50low[[ri]] <- NA logED50high[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "logit" logED50conf <- try(confint(m,"logED50",level=level), silent = TRUE) if (!inherits(logED50conf, "try-error")) { logED50low[[ri]] <- logED50conf[[1]] logED50high[[ri]] <- logED50conf[[2]] } else { logED50low[[ri]] <- NA logED50high[[ri]] <- NA } } } } } if (weibull && length(subset(allWrong,allWrong == i))==0) { m <- try(nls(response ~ pweibull(-log10(dose)+location,shape), data=tmp, algorithm="port", start=list(location=startlogED50[[i]],shape=ws0)), silent = TRUE) if (chooseone==FALSE || fit==FALSE) { if (!inherits(m, "try-error")) { ri <- ri + 1 models[[ri]] <- m 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 if (sqrdev(logED50[[ri]]) > 0.1) { message("\nCan't find ED50 for fitted weibull model of ",i, " data\nwith startlogED50 ", startlogED50[[i]],"\n") ri <- ri - 1 length(a) <- length(b) <- ri length(logED50) <- ri } else { c[[ri]] <- NA fit <- TRUE 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) logED50low[[ri]] <- NA logED50high[[ri]] <- NA if (logED50[[ri]] > rlhd[[ri]]) { mtype[[ri]] <- "no fit" logED50[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA } else { mtype[[ri]] <- "weibull" } } } } } } } 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 logED50low[[ri]] <- NA logED50high[[ri]] <- NA a[[ri]] <- NA b[[ri]] <- NA c[[ri]] <- NA } } results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, logED50low, logED50high, runit, sigma, a, b) lower_level_percent = paste(100 * (1 - level)/2, "%", sep = "") upper_level_percent = paste(100 * (1 + level)/2, "%", sep = "") names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50", lower_level_percent, upper_level_percent, "unit","sigma","a","b") if (linlogit) { results$c <- c } if (showED50) { results[c("ED50", paste("ED50", c(lower_level_percent, upper_level_percent)))] <- 10^results[7:9] } if (!is.null(EDx)) { for (row.i in 1:ri) { if (mtype[[row.i]] %in% c("probit", "logit", "weibull", "linlogit")) { for (ED in EDx) { of <- function(x) { abs(predict(models[[row.i]], data.frame(dose = 10^x)) - (1 - (ED/100))) } # Search over interval starting an order of magnitude below # the lowest dose up to one order of magnitude above the # highest dose o = optimize(of, results[row.i, c("lld", "lhd")] + c(-1, 1)) # Only keep results within the tolerance if ((o$objective) < EDx.tolerance) { logdose.ED = o$minimum results[row.i, paste0("EDx", ED)] <- 10^logdose.ED } } } } } rownames(results) <- 1:ri attr(results, "models") <- models return(results) } # vim: set ts=4 sw=4 expandtab: