path: root/R/drfit.R

drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
        probit = TRUE, logit = FALSE, weibull = FALSE,
        linlogit = FALSE, conf = 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()
    conflogED50 <- 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 <- 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,
                                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
                                conflogED50[[ri]] <- NA
                                a[[ri]] <- NA
                                b[[ri]] <- NA
                                c[[ri]] <- NA
                            } else {
                                mtype[[ri]] <- "linlogit"
                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
                                conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 3)
                                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
                                    conflogED50[[ri]] <- NA
                                    a[[ri]] <- NA
                                    b[[ri]] <- NA
                                } else {
                                    mtype[[ri]] <- "probit"
                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                    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
                                    conflogED50[[ri]] <- NA
                                    a[[ri]] <- NA
                                    b[[ri]] <- NA
                                } else {
                                    mtype[[ri]] <- "logit"
                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                }
                            }
                        }
                    }

                    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
                                    conflogED50[[ri]] <- NA
                                    a[[ri]] <- NA
                                    b[[ri]] <- NA
                                } else {
                                    mtype[[ri]] <- "weibull"
                                    stderrlogED50[[ri]] <- NA
                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                }
                            }
                        }
                    }

                }

            } 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
            conflogED50[[ri]] <- NA
            a[[ri]] <- NA
            b[[ri]] <- NA
            c[[ri]] <- NA
        }
    }
    if (conf)
    {
        results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, conflogED50, runit, sigma, a, b)
        names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","conf","unit","sigma","a","b")
    } else {
        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)
}