R/drcfit.R

Defines functions plotfit2pdf plot.drcfit print.drcfit drcfit

Documented in drcfit plot.drcfit plotfit2pdf print.drcfit

### fit different models to each dose-response curve and choose the best fit
drcfit <- function(itemselect,
                   information.criterion = c("AICc", "BIC", "AIC"),
                   deltaAICminfromnullmodel = 2,
                   postfitfilter = TRUE,
                   preventsfitsoutofrange = TRUE,
                   enablesfequal0inGP = TRUE,
                   enablesfequal0inLGP = TRUE,
                   progressbar = TRUE,
                   parallel = c("no", "snow", "multicore"), ncpus) {
  # Checks
  if (!inherits(itemselect, "itemselect"))
    stop("Use only with 'itemselect' objects, created with the function itemselect.")
  
  parallel <- match.arg(parallel, c("no", "snow", "multicore"))
  if (parallel == "multicore" && .Platform$OS.type == "windows") {
    parallel <- "snow"
    warning(strwrap(prefix = "\n", initial = "\n",
                    "As the multicore option is not supported on Windows it was replaced by snow."))
  }
  if ((parallel == "snow" || parallel == "multicore") && missing(ncpus))
    stop("You have to specify the number of available processors to parallelize the fitting.")
  if (parallel != "no")
    progressbar <- FALSE
  
  if (progressbar)
    cat("The fitting may be long if the number of selected items is high.\n")
  
  # definition of necessary data
  selectindex <- itemselect$selectindex
  adjpvalue <- itemselect$adjpvalue
  dose <- itemselect$omicdata$dose
  doseranks <- as.numeric(as.factor(itemselect$omicdata$dose))
  data <- itemselect$omicdata$data
  data.mean <- itemselect$omicdata$data.mean
  
  containsNA <- itemselect$omicdata$containsNA
  
  # calculations for starting values and other uses
  dosemin <- min(dose)
  dosemax <- max(dose)
  dosemed <- stats::median(dose[dose != 0]) # lintr: assigned but not used
  doseu <- as.numeric(colnames(data.mean)) # sorted unique doses  # lintr: assigned but not used
  
  # number of points per dose-response curve
  npts <- length(dose)
  ndoses <- length(unique(dose))
  lessthan5doses <- ndoses < 5  # lintr: assigned but not used
  
  nselect <- length(selectindex)
  
  # Information criterion definition
  AICdigits <- 2 # number of digits for rounding the AIC values
  information.criterion <- match.arg(information.criterion, c("AICc", "BIC", "AIC"))
  
  # kcrit gives the argument k to pass to function AIC()
  # depending of the number of parameters of the model
  # (1 to 5, corresponding to the index of the vector)
  if (information.criterion == "AIC") {
    kcrit <- rep(2, 5)  # lintr: assigned but not used
  } else if (information.criterion == "AICc") {
    nparwithsigma <- 1:5 + 1
    kcrit <- 2 * npts / (npts - nparwithsigma - 1)  # lintr: assigned but not used
  } else { # BIC last choice
    lnpts <- log(npts)
    kcrit <- rep(lnpts, 5)  # lintr: assigned but not used
  }
  
  # progress bar
  if (progressbar)
    pb <- utils::txtProgressBar(min = 0, max = length(selectindex), style = 3)
  
  # function to fit all the models an choose the best on one item
  ################################################################
  fitoneitem <- function(i) {
    keeplin <- TRUE
    keepExpo <- TRUE
    keepHill <- TRUE
    keepLGauss <- TRUE
    keepGauss <- TRUE
    
    equalcdG <- FALSE # use to define the value of c equal to d in the Gauss4p model if needed
    equalcdLG <- FALSE # use to define the value of c equal to d in the LGauss4p model if needed
    fequal0LG <- FALSE # use to define the value of f at 0 in the LGauss5p model if needed
    fequal0G <- FALSE # use to define the value of f at 0 in the Gauss5p model if needed
    
    signal <- data[selectindex[i], ]
    signalm <- as.vector(data.mean[selectindex[i], ]) # means per dose
    signalmin <- min(signal, na.rm = TRUE)
    signalmax <- max(signal, na.rm = TRUE)
    
    # preparation of data for modelling with nls
    dset <- data.frame(signal = signal, dose = dose, doseranks = doseranks)
    
    if (containsNA)
      if (!all(stats::complete.cases(dset))) {
        if (!all(stats::complete.cases(signalm))) {
          doseu <- doseu[!is.na(signalm)]
          signalm <- signalm[!is.na(signalm)]
        }
        
        # removing lines with NA values for the signal
        dset <- dset[stats::complete.cases(dset$signal), ]
        npts <- nrow(dset)
        ndoses <- length(unique(dset$dose))
        lessthan5doses <- ndoses < 5
        
        # kcrit gives the argument k to pass to function AIC()
        # dependeing of the number of parameters of the model
        # (1 to 5, corresponding to the index of the vector)
        if (information.criterion == "AIC") {
          kcrit <- rep(2, 5)
        } else if (information.criterion == "AICc") {
          nparwithsigma <- 1:5 + 1
          kcrit <- 2 * npts / (npts - nparwithsigma - 1)
        } else { # BIC last choice
          lnpts <- log(npts)
          kcrit <- rep(lnpts, 5)
        }
      }
    
    # for choice of the linear trend (decreasing or increasing)
    modlin <- stats::lm(signal ~ doseranks, data = dset)
    increaseranks <- stats::coef(modlin)[2] >= 0
    increaseminmax <- dset$dose[which.min(dset$signal)] < dset$dose[which.max(dset$signal)]
    
    # for choice of the quadratic trend (Ushape or Umbrella shape)
    modquad <- stats::lm(signal ~ doseranks + I(doseranks^2), data = dset)
    Ushape <- stats::coef(modquad)[3] >= 0
    
    ################ Expo fit (npar = 3) ###############################
    if (keepExpo) {
      # fit of the exponential model with two starting values for abs(e)
      # 0.1*max(dose) or max(dose)
      startExpo3p.1 <- startvalExp3pnls.1(xm = doseu, ym = signalm,
                                          increase = increaseranks,
                                          Ushape = Ushape)
      startExpo3p.2 <- startvalExp3pnls.2(xm = doseu, ym = signalm,
                                          increase = increaseranks,
                                          Ushape = Ushape)
      if ((increaseranks && !Ushape) || (!increaseranks && Ushape)) { # e < 0
        # Fit of the 3 par model
        Expo3p.1 <- suppressWarnings(try(stats::nls(formExp3p, start = startExpo3p.1, data = dset,
                                                    lower = c(-Inf, -Inf, -Inf),
                                                    upper = c(Inf, Inf, 0),
                                                    algorithm = "port"), silent = TRUE))
        
        Expo3p.2 <- suppressWarnings(try(stats::nls(formExp3p, start = startExpo3p.2, data = dset,
                                                    lower = c(-Inf, -Inf, -Inf),
                                                    upper = c(Inf, Inf, 0),
                                                    algorithm = "port"), silent = TRUE))
      } else { # e > 0
        # Fit of the 3 par model
        Expo3p.1 <- suppressWarnings(try(stats::nls(formExp3p, start = startExpo3p.1, data = dset,
                                                    lower = c(-Inf, -Inf, 0),
                                                    algorithm = "port"), silent = TRUE))
        Expo3p.2 <- suppressWarnings(try(stats::nls(formExp3p, start = startExpo3p.2, data = dset,
                                                    lower = c(-Inf, -Inf, 0),
                                                    algorithm = "port"), silent = TRUE))
      }
      #### convergence of both models
      if ((!inherits(Expo3p.1, "try-error")) && (!inherits(Expo3p.2, "try-error"))) {
        AICExpo3p.1 <- round(stats::AIC(Expo3p.1, k = kcrit[3]), digits = AICdigits)
        AICExpo3p.2 <- round(stats::AIC(Expo3p.2, k = kcrit[3]), digits = AICdigits)
        if (AICExpo3p.1 < AICExpo3p.2) {
          Expo <- Expo3p.1
          AICExpoi <- AICExpo3p.1
        } else {
          Expo <- Expo3p.2
          AICExpoi <- AICExpo3p.2
        }
        
      } else if (inherits(Expo3p.1, "try-error") && inherits(Expo3p.2, "try-error")) {
        #### no convergence of both models
        # keepExpo <- FALSE
        AICExpoi <- Inf
        Expo <- Expo3p.1 # we could have given Expo3p.2
        
      } else if ((!inherits(Expo3p.2, "try-error")) && inherits(Expo3p.1, "try-error")) {
        #### convergence only of Expo3p.2
        Expo <- Expo3p.2
        AICExpoi <- round(stats::AIC(Expo3p.2, k = kcrit[3]), digits = AICdigits)
        
      } else if ((!inherits(Expo3p.1, "try-error")) && inherits(Expo3p.2, "try-error")) {
        #### convergence only of Expo3p.1
        Expo <- Expo3p.1
        AICExpoi <- round(stats::AIC(Expo3p.1, k = kcrit[3]), digits = AICdigits)
      }
      
    } else {
      (AICExpoi <- Inf)
    }
    
    ################## Hill fit (npar = 4) ##########################
    if (keepHill) {
      startHill <- startvalHillnls2(x = dose, y = signal, xm = doseu, ym = signalm,
                                    increase = increaseminmax)
      Hill <- suppressWarnings(try(stats::nls(formHill, start = startHill, data = dset,
                                              lower = c(0, -Inf, -Inf, 0), algorithm = "port"), silent = TRUE))
      if (!inherits(Hill, "try-error")) {
        AICHilli <- round(stats::AIC(Hill, k = kcrit[4]), digits = AICdigits)
      } else {
        # keepHill <- FALSE
        AICHilli <- Inf
      }
    } else {
      (AICHilli <- Inf)
    }
    
    ################# LGauss fit ####################
    if (keepLGauss) {
      if (!lessthan5doses) {
        startLGauss5p <- startvalLGauss5pnls(xm = doseu, ym = signalm,
                                             Ushape = Ushape)
        LGauss5p <- suppressWarnings(try(stats::nls(formLGauss5p, start = startLGauss5p, data = dset,
                                                    lower = c(0, -Inf, -Inf, 0, -Inf), algorithm = "port"), silent = TRUE))
        LGauss5psucces <- !inherits(LGauss5p, "try-error")
        # state a failure if the fitted model is out of the range of the signal
        if (LGauss5psucces && preventsfitsoutofrange)
          LGauss5psucces <- !fLGauss5poutofrange(LGauss5p, signalmin, signalmax)
      }
      startLGauss4p <- startvalLGauss4pnls(xm = doseu, ym = signalm,
                                           Ushape = Ushape)
      LGauss4p <- suppressWarnings(try(stats::nls(formLGauss4p, start = startLGauss4p, data = dset,
                                                  lower = c(0, -Inf, 0, -Inf), algorithm = "port"), silent = TRUE))
      LGauss4psucces <- !inherits(LGauss4p, "try-error")
      # state a failure if the fitted model is out of the range of the signal
      if (LGauss4psucces && preventsfitsoutofrange)
        LGauss4psucces <- !fLGauss4poutofrange(LGauss4p, signalmin, signalmax)
      
      if (lessthan5doses) {
        if (LGauss4psucces) {
          equalcdLG <- TRUE
          LGauss <- LGauss4p
          AICLGaussi <- round(stats::AIC(LGauss4p, k = kcrit[4]), digits = AICdigits)
        } else {
          (AICLGaussi <- Inf)
        }
        
      } else { # if (lessthan5doses)
        #### convergence of both models
        if ((LGauss4psucces) && (LGauss5psucces)) {
          AICLGauss4p <- round(stats::AIC(LGauss4p, k = kcrit[4]), digits = AICdigits)
          AICLGauss5p <- round(stats::AIC(LGauss5p, k = kcrit[5]), digits = AICdigits)
          if (AICLGauss5p < AICLGauss4p) {
            LGauss <- LGauss5p
            AICLGaussi <- AICLGauss5p
          } else {
            LGauss <- LGauss4p
            equalcdLG <- TRUE
            AICLGaussi <- AICLGauss4p
          }
          
        } else if ((!LGauss4psucces) && (!LGauss5psucces)) {
          #### no convergence of both models
          # keepLGauss <- FALSE
          AICLGaussi <- Inf
          LGauss <- LGauss5p # we could have given LGauss4p
          
        } else if ((LGauss4psucces) && (!LGauss5psucces)) {
          #### convergence only of LGauss4p
          equalcdLG <- TRUE
          LGauss <- LGauss4p
          AICLGaussi <- round(stats::AIC(LGauss4p, k = kcrit[4]), digits = AICdigits)
          
        } else if ((LGauss5psucces) && (!LGauss4psucces)) {
          #### convergence only of LGauss5p
          LGauss <- LGauss5p
          AICLGaussi <- round(stats::AIC(LGauss5p, k = kcrit[5]), digits = AICdigits)
          
        } else {
          (AICLGaussi <- Inf)
        }
        
        # If LGauss5p chosen, try with f = 0
        if (enablesfequal0inLGP && (is.finite(AICLGaussi)) && (!equalcdLG)) {
          parLG5p <- stats::coef(LGauss)
          startLprobit <- list(b = parLG5p["b"], c = parLG5p["c"], d = parLG5p["d"], e = parLG5p["e"])
          Lprobit <- suppressWarnings(try(stats::nls(formLprobit, start = startLprobit, data = dset,
                                                     lower = c(0, -Inf, -Inf, 0), algorithm = "port"), silent = TRUE))
          if (!inherits(Lprobit, "try-error")) {
            AICwithfat0 <- round(stats::AIC(Lprobit, k = kcrit[4]), digits = AICdigits)
            if (AICwithfat0 <= AICLGaussi) {
              AICLGaussi <- AICwithfat0
              LGauss <- Lprobit
              fequal0LG <- TRUE
            }
          }
        }
      } # END of if (lessthan5doses)
    } # END of if (keepLGauss)
    
    
    
    
    ################### Gauss fit ########################
    if (keepGauss) {
      if (!lessthan5doses) {
        startGauss5p <- startvalGauss5pnls(xm = doseu, ym = signalm,
                                           Ushape = Ushape)
        Gauss5p <- suppressWarnings(try(stats::nls(formGauss5p, start = startGauss5p, data = dset,
                                                   lower = c(0, -Inf, -Inf, 0, -Inf), algorithm = "port"), silent = TRUE))
        Gauss5psucces <- !inherits(Gauss5p, "try-error")
        # state a failure if the fitted model is out of the range of the signal
        if (Gauss5psucces && preventsfitsoutofrange)
          Gauss5psucces <- !fGauss5poutofrange(Gauss5p, signalmin, signalmax)
      }
      startGauss4p <- startvalGauss4pnls(xm = doseu, ym = signalm,
                                         Ushape = Ushape)
      Gauss4p <- suppressWarnings(try(stats::nls(formGauss4p, start = startGauss4p, data = dset,
                                                 lower = c(0, -Inf, 0, -Inf), algorithm = "port"), silent = TRUE))
      Gauss4psucces <- !inherits(Gauss4p, "try-error")
      # state a failure if the fitted model is out of the range of the signal
      if (Gauss4psucces && preventsfitsoutofrange)
        Gauss4psucces <- !fGauss4poutofrange(Gauss4p, signalmin, signalmax)
      
      if (lessthan5doses) {
        if (Gauss4psucces) {
          equalcdG <- TRUE
          Gauss <- Gauss4p
          AICGaussi <- round(stats::AIC(Gauss4p, k = kcrit[4]), digits = AICdigits)
        } else {
          (AICGaussi <- Inf)
        }
      } else { # if (lessthan5doses)
        #### convergence of both models
        if ((Gauss4psucces) && (Gauss5psucces)) {
          AICGauss4p <- round(stats::AIC(Gauss4p, k = kcrit[4]), digits = AICdigits)
          AICGauss5p <- round(stats::AIC(Gauss5p, k = kcrit[5]), digits = AICdigits)
          if (AICGauss5p < AICGauss4p) {
            Gauss <- Gauss5p
            AICGaussi <- AICGauss5p
          } else {
            Gauss <- Gauss4p
            equalcdG <- TRUE
            AICGaussi <- AICGauss4p
          }
          
        } else if ((!Gauss4psucces) && (!Gauss5psucces)) {
          #### no convergence of both models
          # keepGauss <- FALSE
          AICGaussi <- Inf
          Gauss <- Gauss5p # we could have given Gauss4p
          
        } else if ((Gauss4psucces) && (!Gauss5psucces)) {
          #### convergence only of Gauss4p
          equalcdG <- TRUE
          Gauss <- Gauss4p
          AICGaussi <- round(stats::AIC(Gauss4p, k = kcrit[4]), digits = AICdigits)
          
        } else if ((Gauss5psucces) && (!Gauss4psucces)) {
          #### convergence only of Gauss5p
          Gauss <- Gauss5p
          AICGaussi <- round(stats::AIC(Gauss5p, k = kcrit[5]), digits = AICdigits)
        } else {
          (AICGaussi <- Inf)
        }
        
        # If Gauss5p chosen, try with f = 0
        if (enablesfequal0inGP && (is.finite(AICGaussi)) && (!equalcdG)) {
          parG5p <- stats::coef(Gauss)
          startprobit <- list(b = parG5p["b"], c = parG5p["c"], d = parG5p["d"], e = parG5p["e"])
          probit <- suppressWarnings(try(stats::nls(formprobit, start = startprobit, data = dset,
                                                    lower = c(0, -Inf, -Inf, 0), algorithm = "port"), silent = TRUE))
          if (!inherits(probit, "try-error")) {
            AICwithfat0 <- round(stats::AIC(probit, k = kcrit[4]), digits = AICdigits)
            if (AICwithfat0 <= AICGaussi) {
              AICGaussi <- AICwithfat0
              Gauss <- probit
              fequal0G <- TRUE
            }
          }
        }
      } # END of if (lessthan5doses)
    } # END of if (keepGauss)
    
    
    ######### Fit of the linear model ############################
    if (keeplin) {
      lin <- stats::lm(signal ~ dose, data = dset)
      AIClini <- round(stats::AIC(lin, k = kcrit[2]), digits = AICdigits)
    } else {
      (AIClii <- Inf)  # lintr: assigned but not used
    }
    
    
    ######## Fit of the null model (constant) ###########################
    constmodel <- stats::lm(signal ~ 1, data = dset)
    AICconsti <-  round(stats::AIC(constmodel, k = kcrit[1]), digits = AICdigits)
    
    ######### Choice of the best fit #####################################
    ######################################################################
    AICvec <- c(AICGaussi, AICLGaussi, AICHilli, AICExpoi, AIClini)
    # the nb. of the model, 1 to 5, is used in the following with the last
    # being the null model : nb. 6
    indmodeli <- which.min(AICvec)
    AICmin <- AICvec[indmodeli]
    if (AICmin > AICconsti - deltaAICminfromnullmodel) { # we keep the null model
      fit <- constmodel
      indmodeli <- 6 # constant model
      nbpari <- 1
      b.i <- NA
      c.i <- mean(dset$signal)
      d.i <- NA
      e.i <- NA
      f.i <- NA
      SDres.i <- stats::sigma(constmodel)
    } else {
      if (indmodeli == 1) {
        fit <- Gauss
        par <- stats::coef(fit)
        b.i <- par["b"]
        c.i <- ifelse(equalcdG, par["d"], par["c"])
        d.i <- par["d"]
        e.i <- par["e"]
        f.i <- ifelse(fequal0G, 0, par["f"])
        SDres.i <- stats::sigma(fit)
        if (enablesfequal0inGP) {
          nbpari <- ifelse(equalcdG | fequal0G, 4, 5)
        } else {
          nbpari <- ifelse(equalcdG, 4, 5)
        }
        
      } else if (indmodeli == 2) {
        fit <- LGauss
        par <- stats::coef(fit)
        b.i <- par["b"]
        c.i <- ifelse(equalcdLG, par["d"], par["c"])
        d.i <- par["d"]
        e.i <- par["e"]
        f.i <- ifelse(fequal0LG, 0, par["f"])
        SDres.i <- stats::sigma(fit)
        if (enablesfequal0inLGP) {
          nbpari <- ifelse(equalcdLG | fequal0LG, 4, 5)
        } else {
          nbpari <- ifelse(equalcdLG, 4, 5)
        }
        
      } else if (indmodeli == 3) {
        fit <- Hill
        par <- stats::coef(fit)
        b.i <- par["b"]
        c.i <- par["c"]
        d.i <- par["d"]
        e.i <- par["e"]
        f.i <- NA
        SDres.i <- stats::sigma(fit)
        nbpari <- 4
        
      } else if (indmodeli == 4) {
        fit <- Expo
        par <- stats::coef(fit)
        b.i <- par["b"]
        c.i <- NA
        d.i <- par["d"]
        e.i <- par["e"]
        f.i <- NA
        SDres.i <- stats::sigma(fit)
        nbpari <- 3
        
      } else if (indmodeli == 5) {
        fit <- lin
        par <- stats::coef(fit)
        b.i <- par[2]
        c.i <- NA
        d.i <- par[1]
        e.i <- NA
        f.i <- NA
        SDres.i <- stats::sigma(fit)
        nbpari <- 2
      }
    }
    
    # diagnostics on residuals (quadratic trend on residuals)
    # answer to the question: correct mean function ?
    dset$resi <- stats::residuals(fit)
    modquad.resi <- stats::lm(resi ~ doseranks + I(doseranks^2), data = dset)
    mod0.resi <- stats::lm(resi ~ 1, data = dset)
    resimeantrendPi <- stats::anova(modquad.resi, mod0.resi)[[6]][2]
    
    # diagnostics on absolute value of residuals
    # answer to the question of homoscedasticity ?
    # (quadratic trend on abs(residuals))
    dset$absresi <- abs(dset$resi)
    modquad.absresi <- stats::lm(absresi ~ doseranks + I(doseranks^2), data = dset)
    mod0.absresi <- stats::lm(absresi ~ 1, data = dset)
    resivartrendPi <- stats::anova(modquad.absresi, mod0.absresi)[[6]][2]
    
    if (progressbar) {
      utils::setTxtProgressBar(pb, i)
    }
    
    return(c(indmodeli, nbpari, b.i, c.i, d.i, e.i, f.i, SDres.i,
             AIClini, AICExpoi, AICHilli, AICLGaussi,
             AICGaussi, resimeantrendPi, resivartrendPi, AICconsti))
    
  } ##################################### END of fitoneitem
  
  # Loop on items
  # parallel or sequential computation
  if (parallel != "no") {
    if (parallel == "snow") type <- "PSOCK" else if (parallel == "multicore") type <- "FORK"
    clus <- parallel::makeCluster(ncpus, type = type)
    res <- parallel::parSapply(clus, 1:nselect, fitoneitem)
    parallel::stopCluster(clus)
  } else {
    res <- sapply(1:nselect, fitoneitem)
  }
  
  # close progress bar
  if (progressbar)
    close(pb)
  
  dres <- as.data.frame(t(res))
  colnames(dres) <- c("model", "nbpar", "b", "c", "d", "e", "f", "SDres",
                      "InfoCrit.L", "InfoCrit.E", "InfoCrit.H", "InfoCrit.lGP", "InfoCrit.GP",
                      "resimeantrendP", "resivartrendP", "InfoCrit.nullmodel")
  
  dres <- cbind(data.frame(id = row.names(data)[selectindex],
                           irow = selectindex,
                           adjpvalue = adjpvalue),
                dres)
  
  # correction of dres$resimeantrendP and dres$resivartrendP
  # using Benjamini Hochberg method
  # not done to be more cautious and alert about potential heteroscedasticity pb
  # and eliminate bad fits
  # dres$resimeantrendadjP <- p.adjust(dres$resimeantrendP, method = "BH")
  # dres$resivartrendadjP <- p.adjust(dres$resivartrendP, method = "BH")
  
  
  # removing of null models (const, model no 6) and
  # fits eliminated by the quadratic trend test on residuals
  if (postfitfilter) {
    lines.success <- (dres$model != 6) &
      ((dres$resimeantrendP > 0.05) | is.na(dres$resimeantrendP))
    # is.na(resimeantrendP because anova of two models with very close RSS
    # may return NA for pvalue)
  } else {
    lines.success <- (dres$model != 6)
  }
  
  
  dres.failure <- dres[!lines.success, ]
  dfail <- dres.failure[, c("id", "irow", "adjpvalue")]
  dfail$cause <- character(length = nrow(dfail))
  dfail$cause[dres.failure$model == 6] <- "constant.model"
  dfail$cause[dres.failure$model != 6] <- "trend.in.residuals"
  
  dres <- dres[lines.success, ]
  
  # update of nselect
  nselect <- nrow(dres)
  
  dc <- dres[, c("id", "irow", "adjpvalue", "model", "nbpar", "b", "c", "d", "e", "f", "SDres")]
  dresitests <- dres[, c("resimeantrendP", "resivartrendP")]
  modelnames <- c("Gauss-probit", "log-Gauss-probit", "Hill", "exponential", "linear")
  dc$model <- modelnames[dc$model]
  
  # Calculation of the theoretical value at the control : y0
  # of the theoretical value at the maximal dose : yatdosemax
  # of the theoretical signal range on the range of tested concentration : yrange
  # of the maximal change of y from its value at the minimal dose : maxychange
  # of the x-value that corresponds to the extremum for U and bell curves : xextrem
  ##################################################################################
  y0 <- numeric(length = nselect)
  yatdosemax <- numeric(length = nselect)
  yrange <- numeric(length = nselect)
  maxychange <- numeric(length = nselect)
  xextrem <- numeric(length = nselect)
  xextrem[1:nselect] <- NA # will remain at NA for monotonic curves
  yextrem <- numeric(length = nselect)
  yextrem[1:nselect] <- NA # will remain at NA for monotonic curves
  
  # calculation of y0, maxychange and yrange for linear curves
  indlin <- which(dc$model == "linear")
  vb <- dc$b[indlin]
  vd <- dc$d[indlin]
  ydosemin <- flin(dosemin, vb, vd)
  ydosemax <- flin(dosemax, vb, vd)
  yatdosemax[indlin] <- ydosemax
  yrange[indlin] <- abs(ydosemin - ydosemax)
  y0[indlin] <- vd
  maxychange[indlin] <- yrange[indlin]
  
  # calculation of y0 and yrange for exponential curves
  indExpo <- which(dc$model == "exponential")
  vb <- dc$b[indExpo]
  vd <- dc$d[indExpo]
  ve <- dc$e[indExpo]
  ydosemin <- fExpo(dosemin, vb, vd, ve)
  ydosemax <- fExpo(dosemax, vb, vd, ve)
  yatdosemax[indExpo] <- ydosemax
  yrange[indExpo] <- abs(ydosemin - ydosemax)
  y0[indExpo] <- vd
  maxychange[indExpo] <- yrange[indExpo]
  
  # calculation of y0 and yrange for Hill curves
  indHill <- which(dc$model == "Hill")
  vb <- dc$b[indHill]
  vc <- dc$c[indHill]
  vd <- dc$d[indHill]
  ve <- dc$e[indHill]
  ydosemin <- fHill(dosemin, vb, vc, vd, ve)
  ydosemax <- fHill(dosemax, vb, vc, vd, ve)
  yatdosemax[indHill] <- ydosemax
  yrange[indHill] <- abs(ydosemin - ydosemax)
  y0[indHill] <- vd
  maxychange[indHill] <- yrange[indHill]
  
  # calculation of y0, xextrem and yrange for Gauss-probit curves
  # when f != 0
  indGP <- which(dc$model == "Gauss-probit" & dc$f != 0)
  vb <- dc$b[indGP]
  vc <- dc$c[indGP]
  vd <- dc$d[indGP]
  ve <- dc$e[indGP]
  vf <- dc$f[indGP]
  xextr <- xextrem[indGP] <- ve + (vc - vd) * vb / (vf * sqrt(2 * pi))
  yextr <- yextrem[indGP] <- fGauss5p(xextr, vb, vc, vd, ve, vf)
  ydosemin <- fGauss5p(dosemin, vb, vc, vd, ve, vf)
  ydosemax <- fGauss5p(dosemax, vb, vc, vd, ve, vf)
  yatdosemax[indGP] <- ydosemax
  yrange[indGP] <- pmax(abs(ydosemin - yextr), abs(yextr - ydosemax))
  maxychange[indGP] <- pmax(abs(ydosemin - yextr), abs(ydosemax - ydosemin))
  y0[indGP] <- fGauss5p(0, vb, vc, vd, ve, vf)
  
  # when f == 0
  indGPf0 <- which((dc$model == "Gauss-probit" & dc$f == 0))
  vb <- dc$b[indGPf0]
  vc <- dc$c[indGPf0]
  vd <- dc$d[indGPf0]
  ve <- dc$e[indGPf0]
  vf <- dc$f[indGPf0]
  ydosemin <- fGauss5p(dosemin, vb, vc, vd, ve, vf)
  ydosemax <- fGauss5p(dosemax, vb, vc, vd, ve, vf)
  yatdosemax[indGPf0] <- ydosemax
  yrange[indGPf0] <- abs(ydosemin - ydosemax)
  y0[indGPf0] <- fGauss5p(0, vb, vc, vd, ve, vf)
  maxychange[indGPf0] <- yrange[indGPf0]
  
  # calculation of y0, xextrem and yrange for log-Gauss-probit curves
  # when f != 0
  indlGP <- which(dc$model == "log-Gauss-probit" & dc$f != 0)
  vb <- dc$b[indlGP]
  vc <- dc$c[indlGP]
  vd <- dc$d[indlGP]
  ve <- dc$e[indlGP]
  vf <- dc$f[indlGP]
  xextr <- xextrem[indlGP] <- exp(log(ve) + (vc - vd) * vb / (vf * sqrt(2 * pi)))
  yextr <- yextrem[indlGP] <- fLGauss5p(xextr, vb, vc, vd, ve, vf)
  ydosemin <- fLGauss5p(dosemin, vb, vc, vd, ve, vf)
  ydosemax <- fLGauss5p(dosemax, vb, vc, vd, ve, vf)
  yatdosemax[indlGP] <- ydosemax
  yrange[indlGP] <- pmax(abs(ydosemin - yextr), abs(yextr - ydosemax))
  maxychange[indlGP] <- pmax(abs(ydosemin - yextr), abs(ydosemax - ydosemin))
  y0[indlGP] <- vd
  
  # when f == 0
  indlGPf0 <- which(dc$model == "log-Gauss-probit" & dc$f == 0)
  vb <- dc$b[indlGPf0]
  vc <- dc$c[indlGPf0]
  vd <- dc$d[indlGPf0]
  ve <- dc$e[indlGPf0]
  vf <- dc$f[indlGPf0]
  ydosemin <- fLGauss5p(dosemin, vb, vc, vd, ve, vf)
  ydosemax <- fLGauss5p(dosemax, vb, vc, vd, ve, vf)
  yatdosemax[indlGPf0] <- ydosemax
  yrange[indlGPf0] <- abs(ydosemin - ydosemax)
  y0[indlGPf0] <- vd
  maxychange[indlGPf0] <- yrange[indlGPf0]
  
  # definition of trend and typology
  typology <- character(length = nselect)
  trend <- character(length = nselect)
  
  if (nselect != 0) {
    for (i in 1:nselect) {
      di <- dc[i, ]
      if (di$model == "exponential" && di$e > 0 && di$b > 0) {
        typology[i] <- "E.inc.convex"
        trend[i] <- "inc"
      } else if (di$model == "exponential" && di$e <= 0 && di$b > 0) {
        typology[i] <- "E.dec.convex"
        trend[i] <- "dec"
      } else if (di$model == "exponential" && di$e <= 0 && di$b <= 0) {
        typology[i] <- "E.inc.concave"
        trend[i] <- "inc"
      } else if (di$model == "exponential" && di$e > 0 && di$b <= 0) {
        typology[i] <- "E.dec.concave"
        trend[i] <- "dec"
      } else if (di$model == "Hill" && di$c > di$d) {
        typology[i] <- "H.inc"
        trend[i] <- "inc"
      } else if (di$model == "Hill" && di$c <= di$d) {
        typology[i] <- "H.dec"
        trend[i] <- "dec"
      } else if (di$model == "log-Gauss-probit" && di$f < 0) {
        typology[i] <- "lGP.U"
        trend[i] <- "U"
      } else if (di$model == "log-Gauss-probit" && di$f >= 0) {
        typology[i] <- "lGP.bell"
        trend[i] <- "bell"
      } else if (di$model == "Gauss-probit" && di$f < 0) {
        typology[i] <- "GP.U"
        trend[i] <- "U"
      } else if (di$model == "Gauss-probit" && di$f >= 0) {
        typology[i] <- "GP.bell"
        trend[i] <- "bell"
      } else if (di$model == "linear" && di$b > 0) {
        typology[i] <- "L.inc"
        trend[i] <- "inc"
      } else if (di$model == "linear" && di$b <= 0) {
        typology[i] <- "L.dec"
        trend[i] <- "dec"
      }
      
      if (enablesfequal0inLGP) {
        if (di$model == "log-Gauss-probit" && di$f == 0) {
          if (di$c > di$d) {
            typology[i] <- "lGP.inc"
            trend[i] <- "inc"
          } else
            if (di$c <= di$d) {
              typology[i] <- "lGP.dec"
              trend[i] <- "dec"
            }
        }
      }
      if (enablesfequal0inGP) {
        if (di$model == "Gauss-probit" && di$f == 0) {
          if (di$c > di$d) {
            typology[i] <- "GP.inc"
            trend[i] <- "inc"
          } else
            if (di$c <= di$d) {
              typology[i] <- "GP.dec"
              trend[i] <- "dec"
            }
        }
      }
      
    } # END of the for
  } else {
    warning(strwrap(prefix = "\n", initial = "\n",
                    "THERE IS NO SUCCESSFUL FIT."))
  }
  
  dc$typology <- typology
  
  # correction of the trend and typology for Gauss-probit curves with xextrem == 0
  indnullxextr <- which((dc$model == "Gauss-probit") & (xextrem == 0))
  trend[indnullxextr] <- ifelse(dc$f[indnullxextr] > 0, "dec", "inc")
  typology[indnullxextr] <- ifelse(dc$f[indnullxextr] > 0, "GP.dec", "GP.inc")
  
  dc$trend <- factor(trend)
  dc$y0 <- y0
  dc$yatdosemax <- yatdosemax
  dc$yrange <- yrange
  dc$maxychange <- maxychange
  dc$xextrem <- xextrem
  dc$yextrem <- yextrem
  
  # number of null models
  n.failure <- length(itemselect$selectindex) - nrow(dc)
  
  dc$model <- factor(dc$model, # to specify the order
                     levels = c("Hill", "linear", "exponential", "Gauss-probit", "log-Gauss-probit"))
  dc$typology <- factor(dc$typology)
  dInfoCrit <- dres[, c("InfoCrit.L", "InfoCrit.E", "InfoCrit.H", "InfoCrit.lGP", "InfoCrit.GP", "InfoCrit.nullmodel")]
  
  reslist <- list(fitres = dc, omicdata = itemselect$omicdata,
                  information.criterion = information.criterion, information.criterion.val = dInfoCrit,
                  n.failure = n.failure, unfitres = dfail,
                  residualtests = dresitests)
  
  return(structure(reslist, class = "drcfit"))
}

print.drcfit <- function(x, ...) {
  if (!inherits(x, "drcfit"))
    stop("Use only with 'drcfit' objects.")
  
  cat("Results of the fitting using the", x$information.criterion, "to select the best fit model\n")
  ttrend <- table(x$fitres$trend)
  tfit <- table(x$fitres$model)
  nsucces <- nrow(x$fitres)
  nfirstselect <- x$n.failure + nsucces
  if (x$n.failure > 0)
    cat(strwrap(paste(x$n.failure, "dose-response curves out of", nfirstselect, "previously selected were removed 
                       because no model could be fitted reliably.")), fill = TRUE)
  ncaseheterosced <- length(which(x$residualtests$resivartrendP < 0.05))
  ntot <- nrow(x$residualtests)
  pc.heterosced <- round(ncaseheterosced / ntot * 100)
  if (pc.heterosced > 50)
    cat(strwrap(paste0(pc.heterosced, "% of the fitted dose-response curves show a significant heteroscedasticity. 
                       (non constant variance).")), fill = TRUE)
  cat("Distribution of the chosen models among the", nsucces, "fitted dose-response curves:\n")
  print(tfit)
  cat("Distribution of the trends (curve shapes) among the", nsucces, "fitted dose-response curves:\n")
  print(ttrend)
  # ttypology <- table(x$fitres$typology)
  # cat("Distribution of the typology of the ",nsucces," fitted dose-response curves :\n")
  # print(ttypology)
}

plot.drcfit <- function(x, items,
                        plot.type = c("dose_fitted", "dose_residuals", "fitted_residuals"),
                        dose_log_transfo = TRUE,
                        BMDoutput, BMDtype = c("zSD", "xfold"), ...) {
  
  plot.type <- match.arg(plot.type, c("dose_fitted", "dose_residuals", "fitted_residuals"))
  if (!inherits(x, "drcfit"))
    stop("Use only with 'drcfit' objects.")
  
  if (missing(items)) {
    items <- 20
  }
  if (!(is.numeric(items) || is.character(items)))
    stop("Wrong argument 'items'. It must be a number inferior or equal to 20 or
    a character vector indicating the identifiers of the items who want to plot.")
  if (is.numeric(items)) {
    inditems <- seq_len(min(nrow(x$fitres), items))
    subd <- x$fitres[inditems, ]
  } else if (is.character(items)) {
    inditems <- match(items, x$fitres$id)
    if (anyNA(inditems))
      stop("At least one of the chosen items was not selected as responding. You should use targetplot() in that case.")
    subd <- x$fitres[inditems, ]
  }
  subd$id <- factor(subd$id, levels = subd$id)
  g <- plotfitsubset(subd,
                     dose = x$omicdata$dose,
                     data = x$omicdata$data,
                     data.mean = x$omicdata$data.mean,
                     npts = 500,
                     plot.type = plot.type,
                     dose_log_transfo = dose_log_transfo) + theme_classic()
  
  addBMD <- FALSE
  addCI <- FALSE
  ## optional add of BMD values on fits
  if (!(missing(BMDoutput)) && (plot.type == "dose_fitted")) {
    BMDtype <- match.arg(BMDtype, c("zSD", "xfold"))
    addBMD <- TRUE
    if (inherits(BMDoutput, "bmdcalc") || inherits(BMDoutput, "bmdboot")) {
      subbmdres <- BMDoutput$res[inditems, ]
      if (inherits(BMDoutput, "bmdboot"))
        addCI <- TRUE
      else
        addCI <- FALSE
      
      if (any(subd$id != subbmdres$id) || any(subd$yrange != subbmdres$yrange)) {
        warning(strwrap(prefix = "\n", initial = "\n",
                        "To add BMD values on the plot you must 
                        first apply bmdcalc() (and if you want also BMD confidence intervals
                        bmdboot()) on the R object of class drcfit that is specified as first
                        argument of the current plot function, and then 
                        give in the argument BMDoutput the R object given in output
                    of bmdcalc() or bmdboot()."))
        addBMD <- FALSE
      } else {
        if (BMDtype == "zSD") {
          zvalue <- BMDoutput$z
          subbmdres$lowhline <- subbmdres$y0 - zvalue * subbmdres$SDres
          subbmdres$uphline <- subbmdres$y0 + zvalue * subbmdres$SDres
          subbmdres$BMD <- subbmdres$BMD.zSD
          if (inherits(BMDoutput, "bmdboot")) {
            subbmdres$BMDlower <- subbmdres$BMD.zSD.lower
            subbmdres$BMDupper <- subbmdres$BMD.zSD.upper
          }
        } else { # so BMDxfold
          xvalue <- BMDoutput$x
          subbmdres$lowhline <- subbmdres$y0 * (1 + xvalue / 100)
          subbmdres$uphline <- subbmdres$y0 * (1 - xvalue / 100)
          subbmdres$BMD <- subbmdres$BMD.xfold
          if (inherits(BMDoutput, "bmdboot")) {
            addCI <- TRUE
            subbmdres$BMDlower <- subbmdres$BMD.xfold.lower
            subbmdres$BMDupper <- subbmdres$BMD.xfold.upper
          }
        }
      }
    } else {
      warning(strwrap(prefix = "\n", initial = "\n",
                      "To add BMD values on the plot you must 
                        first apply bmdcalc() (and if you want also BMD confidence intervals
                        bmdboot()) on the R object of class drcfit that is specified as first
                        argument of the current plot function, and then 
                        give in the argument BMDoutput the R object given in output
                    of bmdcalc() or bmdboot()."))
      addBMD <- FALSE
    }
  } # END if !missing(BMDoutput)
  
  if (addBMD) {
    subbmdres$id <- factor(subbmdres$id, levels = subd$id)
    g <- g + geom_vline(data = subbmdres,
                        aes(xintercept = .data$BMD),
                        linetype = 1, colour = "red") +
      # geom_ribbon(data = subbmdres,
      #             aes(ymin = .data$lowhline,
      #                  ymax = .data$uphline),
      #             fill = "red", alpha = 0.1)
      geom_hline(data = subbmdres, aes(yintercept = .data$uphline),
                 linetype = 3, colour = "red") +
      geom_hline(data = subbmdres, aes(yintercept = .data$lowhline),
                 linetype = 3, colour = "red")
    if (addCI) {
      g <- g + geom_vline(data = subbmdres, aes(xintercept = .data$BMDlower),
                          linetype = 2, colour = "red") +
        geom_vline(data = subbmdres, aes(xintercept = .data$BMDupper),
                   linetype = 2, colour = "red")
    }
  }
  print(g)
  
}

plotfit2pdf <- function(x, items,
                        plot.type = c("dose_fitted", "dose_residuals", "fitted_residuals"),
                        dose_log_transfo = TRUE,
                        BMDoutput, BMDtype = c("zSD", "xfold"),
                        nrowperpage = 6, ncolperpage = 4,
                        path2figs = getwd()) {
  plot.type <- match.arg(plot.type, c("dose_fitted", "dose_residuals", "fitted_residuals"))
  if (!inherits(x, "drcfit"))
    stop("Use only with 'drcfit' objects.")
  
  # a ggplot alternative
  if (missing(items)) {
    items <- x$fitres$id
  }
  if (!(is.numeric(items) || is.character(items)))
    stop("Wrong argument 'items'. It must be a number inferior or equal to 20 or
        a character vector indicating the identifiers of the items who want to plot.")
  if (is.numeric(items)) {
    inditems <- seq_len(min(nrow(x$fitres), items))
    subd <- x$fitres[seq_len(min(nrow(x$fitres), items)), ]
  } else if (is.character(items)) {
    inditems <- match(items, x$fitres$id)
    if (anyNA(inditems))
      stop("At least one of the chosen items was not selected as responding. You should use targetplot() in that case.")
    subd <- x$fitres[inditems, ]
  }
  
  file2plot <- paste0(path2figs, "/drcfitplot.pdf")
  message(strwrap(prefix = "\n", initial = "\n",
                  paste0("Figures are stored in ", normalizePath(path2figs), ".")))
  grDevices::pdf(file2plot, width = 7, height = 10, onefile = TRUE) # w and h in inches
  
  nplotsperpage <- nrowperpage * ncolperpage
  npage <- ceiling(nrow(subd) / nplotsperpage)
  
  addBMD <- FALSE
  addCI <- FALSE
  if (!(missing(BMDoutput)) && (plot.type == "dose_fitted")) {
    BMDtype <- match.arg(BMDtype, c("zSD", "xfold"))
    addBMD <- TRUE
    if (inherits(BMDoutput, "bmdcalc") || inherits(BMDoutput, "bmdboot")) {
      subbmdres <- BMDoutput$res[inditems, ]
      if (inherits(BMDoutput, "bmdboot"))
        addCI <- TRUE
      else
        addCI <- FALSE
      
      if (any(subd$id != subbmdres$id) || any(subd$yrange != subbmdres$yrange)) {
        warning(strwrap(prefix = "\n", initial = "\n",
                        "To add BMD values on the plot you must 
                        first apply bmdcalc() (and if you want also BMD confidence intervals
                        bmdboot()) on the R object of class drcfit that is specified as first
                        argument of the current plot function, and then 
                        give in the argument BMDoutput the R object given in output
                    of bmdcalc() or bmdboot()."))
        addBMD <- FALSE
      } else {
        if (BMDtype == "zSD") {
          zvalue <- BMDoutput$z
          subbmdres$lowhline <- subbmdres$y0 - zvalue * subbmdres$SDres
          subbmdres$uphline <- subbmdres$y0 + zvalue * subbmdres$SDres
          subbmdres$BMD <- subbmdres$BMD.zSD
          if (inherits(BMDoutput, "bmdboot")) {
            subbmdres$BMDlower <- subbmdres$BMD.zSD.lower
            subbmdres$BMDupper <- subbmdres$BMD.zSD.upper
          }
        } else { # so BMDxfold
          xvalue <- BMDoutput$x
          subbmdres$lowhline <- subbmdres$y0 * (1 + xvalue / 100)
          subbmdres$uphline <- subbmdres$y0 * (1 - xvalue / 100)
          subbmdres$BMD <- subbmdres$BMD.xfold
          if (inherits(BMDoutput, "bmdboot")) {
            addCI <- TRUE
            subbmdres$BMDlower <- subbmdres$BMD.xfold.lower
            subbmdres$BMDupper <- subbmdres$BMD.xfold.upper
          }
        }
      }
    } else {
      warning(strwrap(prefix = "\n", initial = "\n",
                      "To add BMD values on the plot you must 
                        first apply bmdcalc() (and if you want also BMD confidence intervals
                        bmdboot()) on the R object of class drcfit that is specified as first
                        argument of the current plot function, and then 
                        give in the argument BMDoutput the R object given in output
                    of bmdcalc() or bmdboot()."))
      addBMD <- FALSE
    }
  }# END if !missing(BMDoutput)
  
  
  for (i in 1:npage) {
    if (i == npage) indmax <- nrow(subd) else indmax <- i * nplotsperpage
    ind2plot <- seq((i - 1) * nplotsperpage + 1, indmax, 1)
    g <- plotfitsubset(subd[ind2plot, ],
                       dose = x$omicdata$dose,
                       data = x$omicdata$data,
                       data.mean = x$omicdata$data.mean,
                       npts = 500,
                       plot.type = plot.type,
                       dose_log_transfo = dose_log_transfo,
                       nr = nrowperpage,
                       nc = ncolperpage) + theme_classic()
    if (addBMD) {
      data4BMDadd <- subbmdres[ind2plot, ]
      # to keep the ordering of items by p-value and not by alphabetic order
      data4BMDadd$id <- factor(data4BMDadd$id, levels = data4BMDadd$id)
      g <- g + geom_vline(data = data4BMDadd,
                          aes(xintercept = .data$BMD),
                          linetype = 1, colour = "red") +
        # geom_ribbon(data = subbmdres[ind2plot, ],
        #             aes(ymin = .data$lowhline,
        #                  ymax = .data$uphline),
        #             fill = "red", alpha = 0.1)
        geom_hline(data = data4BMDadd, aes(yintercept = .data$uphline),
                   linetype = 3, colour = "red") +
        geom_hline(data = data4BMDadd, aes(yintercept = .data$lowhline),
                   linetype = 3, colour = "red")
      if (addCI) {
        g <- g + geom_vline(data = data4BMDadd, aes(xintercept = .data$BMDlower),
                            linetype = 2, colour = "red") +
          geom_vline(data = data4BMDadd, aes(xintercept = .data$BMDupper),
                     linetype = 2, colour = "red")
      }
    }
    print(g)
  }
  grDevices::dev.off()
}
aursiber/DRomics documentation built on Sept. 12, 2024, 8:52 p.m.