R/dose-response-analysis.R

In QurvE: Robust and User-Friendly Analysis of Growth and Fluorescence Curves

#' Perform a dose-response analysis on response vs. concentration data
#'
#' \code{growth.drFit} serves to determine dose-response curves on every condition in a
#' dataset. The response parameter can be chosen from every physiological parameter in a
#' \code{gcTable} table which is obtained via \code{\link{growth.gcFit}}. \code{\link{growth.drFit}}
#' calls the functions \code{\link{growth.drFitSpline}} and \code{\link{growth.drBootSpline}}, or \code{\link{growth.drFitModel}} to
#' generate a table with estimates for EC50 and respecting statistics.
#'
#' @param gcTable A dataframe containing the data for the dose-response curve estimation. Such table of class \code{gcTable} can be obtained by running \code{\link{growth.gcFit}}.
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#' @param dr.method (Character) Define the method used to perform a dose-responde analysis: smooth spline fit (\code{'spline'}) or model fitting (\code{'model'}).
#' @param dr.model (Character) Provide a list of models from the R package 'drc' to include in the dose-response analysis (if \code{dr.method = 'model'}). If more than one model is provided, the best-fitting model will be chosen based on the Akaike Information Criterion.
#' @param dr.have.atleast (Numeric) Minimum number of different values for the response parameter one should have for estimating a dose response curve. Note: All fit procedures require at least six unique values. Default: \code{6}.
#' @param dr.parameter (Character or numeric) The response parameter in the output table to be used for creating a dose response curve. See \code{\link{growth.drFit}} for further details. Default: \code{'mu.linfit'}, which represents the maximum slope of the linear regression. Typical options include: \code{'mu.linfit'}, \code{'lambda.linfit'}, \code{'dY.linfit'}, \code{'mu.spline'}, \code{'dY.spline'}, \code{'mu.model'}, and \code{'A.model'}.
#' @param smooth.dr (Numeric) Smoothing parameter used in the spline fit by smooth.spline during dose response curve estimation. Usually (not necessesary) in (0; 1]. See \code{\link{smooth.spline}} for further details. Default: \code{NULL}.
#' @param log.x.dr (Logical) Indicates whether \code{ln(x+1)} should be applied to the concentration data of the dose response curves. Default: \code{FALSE}.
#' @param log.y.dr (Logical) Indicates whether \code{ln(y+1)} should be applied to the response data of the dose response curves. Default: \code{FALSE}.
#' @param nboot.dr (Numeric) Defines the number of bootstrap samples for EC50 estimation. Use \code{nboot.dr = 0} to disable bootstrapping. Default: \code{0}.
#'
#' @family growth fitting functions
#'
#' @details
#' Common response parameters used in dose-response analysis:<br><br><b>Linear fit:</b><br>- mu.linfit: Growth rate<br>- lambda.linfit: Lag time<br>- dY.linfit: Density increase<br>- A.linfit: Maximum measurement<br><br><b>Spline fit:</b><br>- mu.spline: Growth rate<br>- lambda.spline: Lag time<br>- A.spline: Maximum measurement<br>- dY.spline: Density increase<br>- integral.spline: Integral<br><br><b>Parametric fit:</b><br>- mu.model: Growth rate<br>- lambda.model: Lag time<br>- A.model: Maximum measurement<br>- integral.model: Integral'
#'
#' @return An object of class \code{drFit}.
#' \item{raw.data}{Data that passed to the function as \code{gcTable}.}
#' \item{drTable}{Dataframe containing condition identifiers, fit options, and results of the dose-response analysis.}
#' \item{drBootSplines}{List of all \code{drBootSpline} objects generated by the call of \code{\link{growth.drBootSpline}} for each distinct experiment.}
#' \item{drFittedSplines}{List of all \code{drFitSpline} objects generated by the call of \code{\link{growth.drFitSpline}} for each distinct experiment.}
#' \item{control}{Object of class \code{grofit.control} containing list of options passed to the function as \code{control}.}
#'
#'
#' @references Matthias Kahm, Guido Hasenbrink, Hella Lichtenberg-Frate, Jost Ludwig, Maik Kschischo (2010). _grofit: Fitting Biological Growth Curves with R_. Journal of Statistical Software, 33(7), 1-21. DOI: 10.18637/jss.v033.i07
#'
#' @export
#' @examples
#' \donttest{
#' # Create random growth data set
#' rnd.data1 <- rdm.data(d = 35, mu = 0.8, A = 5, label = 'Test1')
#' rnd.data2 <- rdm.data(d = 35, mu = 0.6, A = 4.5, label = 'Test2')
#'
#' rnd.data <- list()
#' rnd.data[['time']] <- rbind(rnd.data1$time, rnd.data2$time)
#' rnd.data[['data']] <- rbind(rnd.data1$data, rnd.data2$data)
#'
#' # Run growth curve analysis workflow
#' gcFit <- growth.gcFit(time = rnd.data$time,
#'                        data = rnd.data$data,
#'                        parallelize = FALSE,
#'                        control = growth.control(fit.opt = 's',
#'                                                 suppress.messages = TRUE))
#'
#' # Perform dose-response analysis
#' drFit <- growth.drFit(gcTable = gcFit$gcTable,
#'              control = growth.control(dr.parameter = 'mu.spline'))
#'
#' # Inspect results
#' summary(drFit)
#' plot(drFit)
#' }
growth.drFit <- function(
    gcTable, control = growth.control(
      dr.method = "model", dr.model = c(
        "gammadr", "multi2", "LL.2", "LL.3", "LL.4",
        "LL.5", "W1.2", "W1.3", "W1.4", "W2.2",
        "W2.3", "W2.4", "LL.3u", "LL2.2", "LL2.3",
        "LL2.3u", "LL2.4", "LL2.5", "AR.2", "AR.3",
        "MM.2"
      ),
      dr.have.atleast = 6, dr.parameter = "mu.linear",
      nboot.dr = 0, smooth.dr = NULL, log.x.dr = FALSE,
      log.y.dr = FALSE
    )
)
{
  if (methods::is(control) !=
      "grofit.control" && methods::is(control) !=
      "fl.control")
    stop(
      "control must be of class grofit.control or fl.control!"
    )
  EC50.table <- NULL
  if (is.character(control$dr.parameter))
  {
    dr.parameter <- match(control$dr.parameter, colnames(gcTable))
  }
  FitData <- gcTable[!is.na(gcTable[, 1]),
  ]
  table.tests <- table(
    (FitData[, 1])[which(
      (FitData[, 4] == TRUE) & (is.na(FitData[, dr.parameter]) ==
                                  FALSE)
    )]
  )
  distinct <- names(table.tests)
  EC50 <- list()
  EC50.boot <- list()
  validdata <- cbind(
    as.character(distinct),
    table.tests
  )
  colnames(validdata) <- c("TestID", "Number")
  rownames(validdata) <- rep("     ", times = dim(validdata)[1])
  if (control$suppress.messages == FALSE)
  {
    cat("\n")
    cat(
      "=== EC 50 Estimation ==============================\n"
    )
    cat(
      "---------------------------------------------------\n"
    )
    cat("--> Checking data ...\n")
    cat(
      paste(
        "--> Number of distinct tests found:",
        as.character(length(distinct))
      ),
      "\n"
    )
    cat("--> Valid datasets per test: \n")
    print(validdata, quote = FALSE)
  }
  if (TRUE %in% (table.tests < control$dr.have.atleast))
  {
    if (control$suppress.messages == FALSE){
      cat(
        paste(
          "Warning: following tests have not enough ( <",
          as.character(control$dr.have.atleast - 1),
          ") datasets:\n"
        )
      )
      cat(
        distinct[(table.tests < control$dr.have.atleast)]
      )
      cat("These tests will not be regarded\n")
      distinct <- distinct[table.tests >= control$dr.have.atleast]
    }
  }
  if ((length(distinct)) ==
      0)
  {
    warning(
      paste(
        "There are no tests having enough ( >",
        as.character(control$dr.have.atleast - 1),
        ") datasets!\n"
      )
    )
  } else
  {
    skip <- c()

    for (i in 1:length(distinct))
    {
      conc <- factor(
        (FitData[, 3])[which(FitData[, 1] == distinct[i])]
      )
      test <- (as.numeric(FitData[, dr.parameter]))[FitData[,
                                                            1] == distinct[i]]
      conc <- as.factor(as.character(conc[!is.na(test)]))
      test <- test[!is.na(test)]

      if (length(levels(conc)) <
          4)
      {
        message(
          paste0(
            distinct[i], " does not have enough unique concentrations. A condition must have at least 4 different concentrations to be considered for dose-response analysis."
          )
        )
        skip <- c(skip, i)
        next
      }

      names(test) <- rep(
        names(FitData)[dr.parameter],
        length(test)
      )
      drID <- distinct[i]
      if (control$dr.method == "spline")
      {
        EC50[[i]] <- growth.drFitSpline(conc, test, drID, control)
        if (control$nboot.dr > 0)
        {
          EC50.boot[[i]] <- growth.drBootSpline(conc, test, drID, control)
        } else
        {
          EC50.boot[[i]] <- list(
            raw.time = conc, raw.data = test,
            drID = drID, boot.x = NA, boot.y = NA,
            boot.drSpline = NA, ec50.boot = NA,
            bootFlag = FALSE, control = control
          )
          class(EC50.boot[[i]]) <- "drBootSpline"
        }
      } else
      {
        EC50[[i]] <- growth.drFitModel(conc, test, drID, control)
      }
      description <- data.frame(
        Test = distinct[i], log.x = control$log.x.dr,
        log.y = control$log.y.dr, Samples = control$nboot.dr
      )
      if (control$dr.method == "spline")
        out.row <- cbind(
          description, summary.drFitSpline(EC50[[i]]),
          summary.drBootSpline(EC50.boot[[i]])
        ) else out.row <- cbind(description, summary.drFitModel(EC50[[i]]))

      if (!is.na(out.row[1, 5]) &&
          out.row[1, 5] != "NA")
        EC50.table <- rbind(
          as.data.frame(EC50.table),
          out.row
        ) else
        {
          return(NA)
        }

      class(EC50.table) <- c("drTable", "list")
    }
  }
  if (exists("skip") &&
      !is.null(skip))
  {
    distinct <- distinct[-skip]
    EC50 <- EC50[-skip]
    EC50.boot <- EC50.boot[-skip]
  }
  if (control$dr.method == "spline")
  {
    if (length(EC50) ==
        0)
    {
      return(NA)
    }
    names(EC50) <- names(EC50.boot) <- distinct
    fitflags <- unlist(
      lapply(
        1:length(EC50),
        function(x) EC50[[x]]$fitFlag
      )
    )
    if (all(isFALSE(fitflags)))
    {
      drFit <- NA
    } else
    {
      drFit <- list(
        raw.data = FitData, drTable = EC50.table,
        drBootSplines = EC50.boot, drFittedSplines = EC50,
        control = control
      )
    }
  } else
  {
    if (length(EC50) ==
        0)
    {
      return(NA)
    }
    names(EC50) <- distinct
    fitflags <- unlist(
      lapply(
        1:length(EC50),
        function(x) EC50[[x]]$fitFlag
      )
    )
    if (all(isFALSE(fitflags)))
    {
      drFit <- NA
    } else
    {
      drFit <- list(
        raw.data = FitData, drTable = EC50.table,
        drFittedModels = EC50, control = control
      )
    }
  }

  class(drFit) <- "drFit"
  invisible(drFit)
}

#' Perform a smooth spline fit on response vs. concentration data of a single sample to determine the EC50.
#'
#' \code{growth.drFitSpline} performs a smooth spline fit determines the EC50 as the concentration
#' at the half-maximum value of the response parameter of the spline.
#'
#' @param conc Vector of concentration values.
#' @param test Vector of response parameter values of the same length as \code{conc}.
#' @param drID (Character) The name of the analyzed condition
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#'
#' @family dose-response analysis functions
#'
#' @return A \code{drFitSpline} object.
#' \item{raw.conc}{Raw data provided to the function as \code{conc}.}
#' \item{raw.test}{Raw data for the response parameter provided to the function as \code{test}.}
#' \item{drID}{(Character) Identifies the tested condition}
#' \item{fit.conc}{Fitted concentration values.}
#' \item{fit.test}{Fitted response values.}
#' \item{spline}{\code{smooth.spline} object generated by the \code{\link{smooth.spline}} function.}
#' \item{spline.low}{\code{x} and \code{y} values of \code{\link{lowess}} spline fit on raw data. Used to call \code{\link{smooth.spline}}.}
#' \item{parameters}{List of parameters estimated from dose response curve fit.}
#' \itemize{
#' \item \code{EC50}: Concentration at half-maximal response.
#' \item \code{yEC50}: Response value related to EC50.
#' \item \code{EC50.orig}: EC50 value in original scale, if a transformation was applied.
#' \item \code{yEC50.orig}: Response value for EC50 in original scale, if a transformation was applied.
#' }
#' \item{fitFlag}{(Logical) Indicates whether a spline could fitted successfully to data.}
#' \item{reliable}{(Logical) Indicates whether the performed fit is reliable (to be set manually).}
#' \item{control}{Object of class \code{grofit.control} containing list of options passed to the function as \code{control}.}
#' Use \code{\link{plot.drFitSpline}} to visualize the spline fit.
#'
#' @details During the spline fit with \code{\link{smooth.spline}}, higher weights are
#'   assigned to data points with a concentration value of 0, as well as to x-y pairs with
#'   a response parameter value of 0 and pairs at concentration values before
#'   zero-response parameter values.
#'
#' @references Matthias Kahm, Guido Hasenbrink, Hella Lichtenberg-Frate, Jost Ludwig, Maik Kschischo (2010). _grofit: Fitting Biological Growth Curves with R_. Journal of Statistical Software, 33(7), 1-21. DOI: 10.18637/jss.v033.i07
#' @references Christian Ritz, Florent Baty, Jens C. Streibig, Daniel Gerhard (2015). _Dose-Response Analysis Using R_. PLoS ONE 10(12): e0146021. DOI: 10.1371/journal.pone.0146021
#'
#' @export
#' @examples
#' conc <- c(0, rev(unlist(lapply(1:18, function(x) 10*(2/3)^x))),10)
#' response <- c(1/(1+exp(-0.7*(4-conc[-20])))+rnorm(19)/50, 0)
#'
#' TestRun <- growth.drFitSpline(conc, response, drID = 'test',
#'               control = growth.control(log.x.dr = TRUE, smooth.dr = 0.8))
#'
#' print(summary(TestRun))
#'
#' plot(TestRun)
growth.drFitSpline <- function(
    conc, test, drID = "undefined", control = growth.control()
)
{
  if (methods::is(control) !=
      "grofit.control" && methods::is(control) !=
      "fl.control")
    stop(
      "control must be of class grofit.control or fl.control!"
    )
  test.nm <- names(test)[1]
  test <- as.vector(as.numeric(as.matrix(test)))
  conc <- as.vector(as.numeric(as.matrix(conc)))
  if (is.vector(conc) ==
      FALSE || is.vector(test) ==
      FALSE)
    stop(
      "growth.drFitSpline: dose or response data must be a vector !"
    )
  if (control$neg.nan.act == FALSE)
  {
    missings <- is.na(conc) |
      is.na(test) |
      !is.numeric(conc) |
      !is.numeric(test)
    conc <- conc[!missings]
    test <- test[!missings]
    negs <- (conc < 0) | (test < 0)
    conc <- conc[!negs]
    test <- test[!negs]
  } else
  {
    if (sum(
      is.na(conc) |
      is.na(test)
    ))
      stop(
        "growth.drFitSpline: NA values encountered. Program terminated"
      )
    if ((sum((conc < 0)) >
         0) | (sum((test < 0)) >
               0))
      stop(
        "growth.drFitSpline: Negative values encountered. Program terminated"
      )
    if ((FALSE %in% is.numeric(conc)) ||
        (FALSE %in% is.numeric(test)))
      stop(
        "growth.drFitSpline: Non numeric values encountered. Program terminated"
      )
  }
  if (length(test) <
      6)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drFitSpline: There is not enough valid data. Must have at least 6 unique values!"
      )
    drFitSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, EC50.orig = NA,
        yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitSpline) <- "drFitSpline"
    return(drFitSpline)
  }
  if (length(test) <
      control$dr.have.atleast)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drFitSpline: number of valid data points is below the number specified in 'dr.have.atleast'. See growth.control()."
      )
    drFitSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, EC50.orig = NA,
        yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitSpline) <- "drFitSpline"
    return(drFitSpline)
  }
  if (control$log.x.dr == TRUE)
    conc.log <- log(conc + 1)
  if (control$log.y.dr == TRUE)
    test.log <- log(test + 1)
  if (control$log.x.dr == TRUE)
  {
    conc.fit <- log(conc + 1)
  } else
  {
    conc.fit <- conc
  }
  if (control$log.y.dr == TRUE)
  {
    test.fit <- log(test + 1)
  } else
  {
    test.fit <- test
  }
  spltest <- NULL
  fitFlag <- TRUE

  # perform first lowess spline fit followed by
  # smooth.spline onf lowess results
  try(
    spltest.low <- lowess(conc.fit, test.fit, f = 0.1),
    silent = TRUE
  )
  try(
    spltest.low$y <- spltest.low$y[!duplicated(spltest.low$x)],
    silent = TRUE
  )
  try(
    spltest.low$x <- spltest.low$x[!duplicated(spltest.low$x)],
    silent = TRUE
  )

  # assign high weights to concentration values
  # of 0, response values of 0, and response
  # values immediately before 0
  weights <- sapply(
    1:(length(spltest.low$x) -
         1), function(i) ifelse(
           spltest.low$x[[i]] == 0 | spltest.low$y[[i]] ==
             0 | spltest.low$y[[i + 1]] == 0, 1,
           0.05
         )
  )
  weights <- c(
    weights, ifelse(
      spltest.low$x[length(spltest.low$x)] ==
        0 | spltest.low$y[length(spltest.low$y)] ==
        0, 1, 0.05
    )
  )

  try(
    spltest.smooth <- smooth.spline(
      spltest.low$x, spltest.low$y, spar = control$smooth.dr,
      keep.data = FALSE, w = weights
    ),
    silent = TRUE
  )

  if (!exists("spltest.smooth") ||
      is.null(spltest.smooth) ==
      TRUE)
  {
    message(
      "Spline could not be fitted in dose-response analysis!\n"
    )
    fitFlag <- FALSE
    if (is.null(control$smooth.dr) ==
        TRUE)
    {
      message(
        "This might be caused by usage of smoothing parameter 'smooth.dr = NULL'. Re-running the function might solve the problem. If not, please specify 'smooth.dr'.\n"
      )
    }
    drFitSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, EC50.orig = NA,
        yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitSpline) <- "drFitSpline"
    return(drFitSpline)
  }
  conc.min <- min(conc.fit)
  conc.max <- max(conc.fit)
  c.pred <- seq(conc.min, conc.max, length.out = 1000)
  ytest <- stats::predict(spltest.smooth, c.pred)
  yEC.test <- (max(ytest$y) -
                 min(ytest$y))/2 +
    min(ytest$y)
  last.test <- max(ytest$y)
  kec.test <- 1
  for (k in 1:(length(c.pred) -
               1))
  {
    d1 <- (ytest$y[k] - yEC.test)
    d2 <- (ytest$y[k + 1] - yEC.test)
    if (((d1 <= 0) && (d2 >= 0)) | ((d1 >= 0) &&
                                    (d2 <= 0)))
    {
      kec.test <- k
      break
    }
  }
  EC.test <- c.pred[kec.test]
  if (control$suppress.messages == FALSE)
  {
    cat(
      "\n\n=== Dose response curve estimation ================\n"
    )
    cat(
      "--- EC 50 -----------------------------------------\n"
    )
    cat(paste("-->", as.character(drID)))
    cat("\n")
    cat(
      paste(
        c("xEC50", "yEC50"),
        c(EC.test, yEC.test)
      )
    )
  }
  if ((control$log.x.dr == TRUE) && (control$log.y.dr ==
                                     FALSE))
  {
    if (control$suppress.messages == FALSE)
    {
      cat("\n--> Original scale \n")
      cat(
        paste(
          c("xEC50", "yEC50"),
          c(
            exp(EC.test) -
              1, yEC.test
          )
        )
      )
    }
    EC.orig <- c(
      exp(EC.test) -
        1, yEC.test
    )
  } else
  {
    if ((control$log.x.dr == FALSE) && (control$log.y.dr ==
                                        TRUE))
    {
      if (control$suppress.messages == FALSE)
      {
        cat("\n--> Original scale \n")
        cat(
          paste(
            c("xEC50", "yEC50"),
            c(
              EC.test, exp(yEC.test) -
                1
            )
          )
        )
      }
      EC.orig <- c(
        EC.test, exp(yEC.test) -
          1
      )
    } else
    {
      if ((control$log.x.dr == TRUE) && (control$log.y.dr ==
                                         TRUE))
      {
        if (control$suppress.messages == FALSE)
        {
          cat("\n--> Original scale \n")
          cat(
            paste(
              c("xEC50", "yEC50"),
              c(
                exp(EC.test) -
                  1, exp(yEC.test) -
                  1
              )
            )
          )
        }
        EC.orig <- c(
          exp(EC.test) -
            1, exp(yEC.test) -
            1
        )
      } else
      {
        if ((control$log.x.dr == FALSE) &&
            (control$log.y.dr == FALSE))
        {
          EC.orig <- c(EC.test, yEC.test)
        }
      }
    }
  }
  if (control$suppress.messages == FALSE)
  {
    cat("\n\n\n")
  }
  drFitSpline <- list(
    raw.conc = conc, raw.test = test, drID = drID,
    fit.conc = ytest$x, fit.test = ytest$y, spline = spltest.smooth,
    spline.low = spltest.low, parameters = list(
      EC50 = EC.test[1], yEC50 = yEC.test, EC50.orig = EC.orig[1],
      yEC50.orig = EC.orig[2], test = test.nm
    ),
    fitFlag = fitFlag, reliable = NULL, control = control
  )
  class(drFitSpline) <- "drFitSpline"
  invisible(drFitSpline)
}

#' Perform a smooth spline fit on response vs. concentration data of a single sample
#'
#' \code{growth.drBootSpline} resamples the values in a dataset with replacement and performs a spline fit for each bootstrap sample to determine the EC50.
#'
#' @param conc Vector of concentration values.
#' @param test Vector of response parameter values of the same length as \code{conc}.
#' @param drID (Character) The name of the analyzed sample.
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#'
#' @family dose-response analysis functions
#'
#' @return An object of class \code{drBootSpline} containing a distribution of growth parameters and
#'   a \code{drFitSpline} object for each bootstrap sample. Use \code{\link{plot.drBootSpline}}
#'   to visualize all bootstrapping splines as well as the distribution of EC50.
#' \item{raw.conc}{Raw data provided to the function as \code{conc}.}
#' \item{raw.test}{Raw data for the response parameter provided to the function as \code{test}.}
#' \item{drID}{(Character) Identifies the tested condition.}
#' \item{boot.conc}{Table of concentration values per column, resulting from each spline fit of the bootstrap.}
#' \item{boot.test}{Table of response values per column, resulting from each spline fit of the bootstrap.}
#' \item{boot.drSpline}{List containing all \code{drFitSpline} objects generated by the call of \code{\link{growth.drFitSpline}}.}
#' \item{ec50.boot}{Vector of estimated EC50 values from each bootstrap entry.}
#' \item{ec50y.boot}{Vector of estimated response at EC50 values from each bootstrap entry.}
#' \item{BootFlag}{(Logical) Indicates the success of the bootstrapping operation.}
#' \item{control}{Object of class \code{grofit.control} containing list of options passed to the function as \code{control}.}
#'
#' @references Matthias Kahm, Guido Hasenbrink, Hella Lichtenberg-Frate, Jost Ludwig, Maik Kschischo (2010). _grofit: Fitting Biological Growth Curves with R_. Journal of Statistical Software, 33(7), 1-21. DOI: 10.18637/jss.v033.i07
#'
#' @export
#'
#' @examples
#' conc <- c(0, rev(unlist(lapply(1:18, function(x) 10*(2/3)^x))),10)
#' response <- c(1/(1+exp(-0.7*(4-conc[-20])))+rnorm(19)/50, 0)
#'
#' TestRun <- growth.drBootSpline(conc, response, drID = 'test',
#'                control = growth.control(log.x.dr = TRUE, smooth.dr = 0.8,
#'                                         nboot.dr = 50))
#'
#' print(summary(TestRun))
#' plot(TestRun, combine = TRUE)
#'
growth.drBootSpline <- function(
    conc, test, drID = "undefined", control = growth.control()
)
{
  test <- as.vector(as.numeric(as.matrix(test)))
  conc <- as.vector(as.numeric(as.matrix(conc)))
  if (is.vector(conc) ==
      FALSE || is.vector(test) ==
      FALSE)
    stop("Need concentration and treatment !")
  if (methods::is(control) !=
      "grofit.control" && methods::is(control) !=
      "fl.control")
    stop(
      "control must be of class grofit.control or fl.control!"
    )
  if (control$nboot.dr == 0)
    stop(
      "Number of bootstrap samples is zero! See ?growth.control or ?fl.control."
    )
  if (control$neg.nan.act == FALSE)
  {
    missings <- is.na(conc) |
      is.na(test) |
      !is.numeric(conc) |
      !is.numeric(test)
    conc <- conc[!missings]
    test <- test[!missings]
    negs <- (conc < 0) | (test < 0)
    conc <- conc[!negs]
    test <- test[!negs]
  }
  # Test if there are enough unique x-values
  # (conc) to perform spline fit
  if (length(unique(conc)) <
      4)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drBootSpline: There are not enough concentration values. Must have at least 4 unique values!"
      )
    drBootSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      boot.conc = NA, boot.test = NA, boot.drSpline = NA,
      ec50.boot = NA, ec50y.boot = NA, bootFlag = FALSE,
      control = control
    )
    class(drBootSpline) <- "drBootSpline"
    return(drBootSpline)
  } else
  {
    if (sum(
      is.na(conc) |
      is.na(test)
    ))
      stop("NA values encountered. Program terminated")
    if ((sum((conc < 0)) >
         0) | (sum((test < 0)) >
               0))
      stop(
        "growth.drFitSpline: Negative values encountered. Program terminated"
      )
    if ((FALSE %in% is.numeric(conc)) ||
        (FALSE %in% is.numeric(test)))
      stop(
        "growth.drFitSpline: Non numeric values encountered. Program terminated"
      )
  }
  if (length(test) <
      6)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drBootSpline: There is not enough valid data. Must have at least 6 unique values!"
      )
    drBootSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      boot.conc = NA, boot.test = NA, boot.drSpline = NA,
      ec50.boot = NA, ec50y.boot = NA, bootFlag = FALSE,
      control = control
    )
    class(drBootSpline) <- "drBootSpline"
    return(drBootSpline)
  }
  if (length(test) <
      control$dr.have.atleast)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drBootSpline: number of valid data points is below the number specified in 'dr.have.atleast'. See growth.control()."
      )
    drBootSpline <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      boot.conc = NA, boot.test = NA, boot.drSpline = NA,
      ec50.boot = NA, ec50y.boot = NA, bootFlag = FALSE,
      control = control
    )
    class(drBootSpline) <- "drBootSpline"
    return(drBootSpline)
  }

  # /// transformation of data...
  if (control$log.x.dr == TRUE)
    conc.log <- log(conc + 1)
  if (control$log.y.dr == TRUE)
    test.log <- log(test + 1)
  if (control$log.x.dr == TRUE)
    conc.boot <- log(conc + 1) else conc.boot <- conc
  if (control$log.y.dr == TRUE)
    test.boot <- log(test + 1) else test.boot <- test

  # /// Initialize some variables
  boot.x <- array(NA, c(control$nboot.dr, 1000))
  boot.y <- array(NA, c(control$nboot.dr, 1000))
  ECtest.boot <- seq(0, 0, length.out = control$nboot.dr)
  y.EC50.boot <- seq(0, 0, length.out = control$nboot.dr)
  splinefit <- list()
  sa <- seq(1, length(conc.boot))

  # /// begin bootstrapping
  for (b in 1:control$nboot.dr)
  {
    s <- sample(
      sa, length(conc.boot),
      replace = TRUE
    )
    s.conc <- conc.boot[s]
    while (length(unique(s.conc)) <
           4)
    {
      s <- sample(
        sa, length(conc.boot),
        replace = TRUE
      )
      s.conc <- conc.boot[s]
    }
    s.test <- test.boot[s]
    spltest <- NULL
    control.changed <- control
    control.changed$suppress.messages <- TRUE
    splinefit[[b]] <- growth.drFitSpline(s.conc, s.test, drID, control.changed)
    spltest <- splinefit[[b]]$spline
    boot.x[b, 1:length(splinefit[[b]]$fit.conc)] <- splinefit[[b]]$fit.conc
    boot.y[b, 1:length(splinefit[[b]]$fit.test)] <- splinefit[[b]]$fit.test
    if (is.null(spltest) ==
        TRUE)
    {
      message(
        "Spline could not be fitted in dose-response analysis!!\n"
      )
      if (is.null(control$smooth.dr) ==
          TRUE)
      {
        message(
          "This might be caused by usage of smoothing parameter 'smooth.dr = NULL'.\n"
        )
      }
      stop("Error in drBootSpline")
    }
    ECtest.boot[b] <- splinefit[[b]]$parameters$EC50
    y.EC50.boot[b] <- splinefit[[b]]$parameters$yEC50
  }
  ECtest.boot[which(!is.finite(ECtest.boot))] <- NA
  if (control$clean.bootstrap == TRUE)
    ECtest.boot[which(ECtest.boot < 0)] <- NA
  m.test <- mean(ECtest.boot, na.rm = TRUE)
  s.test <- sd(ECtest.boot, na.rm = TRUE)
  if (control$suppress.messages == FALSE)
  {
    cat(
      "=== Bootstrapping of dose response curve ==========\n"
    )
    cat(
      "--- EC 50 -----------------------------------------\n"
    )
    cat("\n")
    cat(
      paste(
        "Mean  : ", as.character(m.test),
        "StDev : ", as.character(s.test),
        "\n"
      )
    )
    cat(
      paste(
        "90% CI: ", as.character(
          c(
            m.test - 1.645 * s.test/control$nboot.dr,
            m.test + 1.645 * s.test/control$nboot.dr
          )
        )
      )
    )
    cat("\n")
    cat(
      paste(
        "95% CI: ", as.character(
          c(
            m.test - 1.96 * s.test/control$nboot.dr,
            m.test + 1.96 * s.test/control$nboot.dr
          )
        )
      )
    )
    cat("\n\n")
  }
  EC50 <- data.frame(
    EC50.boot = m.test, EC50.sd = s.test, CI.90.lo = m.test -
      1.645 * s.test, CI.90.up = m.test + 1.645 *
      s.test, CI.95.lo = m.test - 1.96 * s.test,
    CI.95.up = m.test + 1.96 * s.test
  )
  if (control$log.x.dr == TRUE && control$suppress.messages ==
      FALSE)
  {
    cat("\n")
    cat(
      "--- EC 50 in original scale -----------------------\n"
    )
    cat("\n")
    cat(
      paste(
        "Mean  : ", as.character(
          exp(m.test) -
            1
        ),
        "\n"
      )
    )
    cat(
      paste(
        "90% CI: ", as.character(
          c(
            exp(m.test - 1.645 * s.test/control$nboot.dr) -
              1, exp(m.test + 1.645 * s.test/control$nboot.dr) -
              1
          )
        )
      )
    )
    cat("\n")
    cat(
      paste(
        "95% CI: ", as.character(
          c(
            exp(m.test - 1.96 * s.test/control$nboot.dr) -
              1, exp(m.test + 1.96 * s.test/control$nboot.dr) -
              1
          )
        )
      )
    )
    cat("\n\n")
  } # if (control$log.x.dr == TRUE && control$suppress.messages == FALSE)
  drBootSpline <- list(
    raw.conc = conc, raw.test = test, drID = drID,
    boot.conc = boot.x, boot.test = boot.y, boot.drSpline = splinefit,
    ec50.boot = ECtest.boot, ec50y.boot = y.EC50.boot,
    bootFlag = TRUE, control = control
  )
  class(drBootSpline) <- "drBootSpline"
  invisible(drBootSpline)
}


#' Fit various models to response vs. concentration data of a single sample to determine the EC50.
#'
#' @param conc Vector of concentration values.
#' @param test Vector of response parameter values of the same length as \code{conc}.
#' @param drID (Character) The name of the analyzed condition
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#'
#' @return A \code{drFitModel} object.
#' @references Christian Ritz, Florent Baty, Jens C. Streibig, Daniel Gerhard (2015). _Dose-Response Analysis Using R_. PLoS ONE 10(12): e0146021. DOI: 10.1371/journal.pone.0146021
#' @export
#'
#' @import drc
#'
#' @examples
#' conc <- c(0, rev(unlist(lapply(1:18, function(x) 10*(2/3)^x))),10)
#' response <- c(1/(1+exp(-0.7*(4-conc[-20])))+rnorm(19)/50, 0)
#'
#' TestRun <- growth.drFitModel(conc, response, drID = 'test')
#'
#' print(summary(TestRun))
#' plot(TestRun)
growth.drFitModel <- function(
    conc, test, drID = "undefined", control = growth.control()
)
{
  if (methods::is(control) !=
      "grofit.control" && methods::is(control) !=
      "fl.control")
    stop(
      "growth.drFitModel: control must be of class grofit.control or fl.control!"
    )
  test.nm <- control$dr.parameter
  if (test.nm != control$dr.parameter && control$dr.parameter !=
      "mu.linfit")
    test.nm <- control$dr.parameter
  test <- as.vector(as.numeric(as.matrix(test)))
  conc <- as.vector(as.numeric(as.matrix(conc)))
  models <- control$dr.model
  if (is.vector(conc) ==
      FALSE || is.vector(test) ==
      FALSE)
    stop(
      "growth.drFitModel: dose or response data must be a vector !"
    )
  if (control$neg.nan.act == FALSE)
  {
    missings <- is.na(conc) |
      is.na(test) |
      !is.numeric(conc) |
      !is.numeric(test)
    conc <- conc[!missings]
    test <- test[!missings]
    negs <- (conc < 0) | (test < 0)
    conc <- conc[!negs]
    test <- test[!negs]
  } else
  {
    if (sum(
      is.na(conc) |
      is.na(test)
    ))
      stop(
        "growth.drFitModel: NA values encountered. Program terminated"
      )
    if ((sum((conc < 0)) >
         0) | (sum((test < 0)) >
               0))
      stop(
        "growth.drFitModel: Negative values encountered. Program terminated"
      )
    if ((FALSE %in% is.numeric(conc)) ||
        (FALSE %in% is.numeric(test)))
      stop(
        "growth.drFitModel: Non numeric values encountered. Program terminated"
      )
  }
  if (length(test) <
      6)
  {
    if (control$suppress.messages == FALSE)
      message(
        "growth.drFitModel: There is not enough valid data. Must have at least 6 unique values!"
      )
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, Vmax = NA, Km = NA,
        EC50.orig = NA, yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitModel"
    return(drFitModel)
  }
  if (length(test) <
      control$dr.have.atleast)
  {
    if (control$suppress.messages == FALSE)
      message(
        "growth.drFitModel: number of valid data points is below the number specified in 'dr.have.atleast'. See growth.control()."
      )
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, Vmax = NA, Km = NA,
        EC50.orig = NA, yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitModel"
    return(drFitModel)
  }
  if (control$dr.method == "model.MM" && !all(models %in% c("MM.2", "MM.3")))
  {
    models <- c("MM.2", "MM.3")
  }
  # Perform model fits
  model.fits <- list()
  for (i in 1:length(models))
  {
    model.fits[[i]] <- try(
      suppressWarnings(
        invisible(
          drc::drm(
            test ~ as.numeric(as.character(conc)),
            fct = get(models[i])(),
            control = drc::drmc(
              errorm = FALSE, noMessage = TRUE,
              otrace = TRUE
            )
          )
        )
      ),
      silent = T
    )
  }
  models <- models[!unlist(
    lapply(
      1:length(model.fits),
      function(x) class(model.fits[[x]])
    )
  ) %in%
    c("try-error", "list")]
  model.fits <- model.fits[!unlist(
    lapply(
      1:length(model.fits),
      function(x) class(model.fits[[x]])
    )
  ) %in%
    c("try-error", "list")]
  names(model.fits) <- models

  if (length(model.fits) <
      1)
  {
    if (control$suppress.messages == FALSE)
      message(
        "growth.drFitModel: No DR model could be fit to the test and concentration data."
      )
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, Vmax = NA, Km = NA,
        EC50.orig = NA, yEC50.orig = NA, test = test.nm
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitModel"
    return(drFitModel)
  }
  # select best fitting model
  model.AIC <- lapply(
    1:length(model.fits),
    function(x) AIC(model.fits[[x]])
  )
  names(model.AIC) <- models

  best.model.ndx <- which.min(unlist(model.AIC))
  best.model.nm <- models[best.model.ndx]
  best.model <- model.fits[[best.model.ndx]]

  # get EC50 value
  if (any(grep("BC", best.model.nm)))
  {
    ec50 <- drc::ED(
      best.model, 50, lower = 0.1, upper = 1000,
      interval = "fls", display = FALSE
    )
  } else if (any(grep("NEC", best.model.nm)))
  {
    ec50 <- best.model$coefficients[4]
  } else
  {
    ec50 <- drc::ED(
      best.model, 50, interval = "delta", display = FALSE
    )
  }
  # get response at ec50
  dataList <- best.model[["dataList"]]
  dose <- dataList[["dose"]]

  concgrid <- seq(
    min(dose),
    max(dose),
    length = 200
  )
  respgrid <- best.model$curve[[1]](concgrid)

  y.ec50 <- drc::PR(object = best.model, xVec = ec50[1])

  # drc.models <- drc::getMeanFunctions(display
  # = FALSE) drc.models.nm <-
  # c(unlist(lapply(1:length(drc.models),
  # function(x) drc.models[[x]][1])), 'NEC.4')
  # drc.models.descr <-
  # c(unlist(lapply(1:length(drc.models),
  # function(x) drc.models[[x]][2])), 'model
  # for estimation of\nno effect concentration
  # (NEC)')
  if (control$dr.method == "model.MM")
  {
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = concgrid, fit.test = respgrid,
      model = best.model, parameters = list(
        EC50 = ec50, yEC50 = y.ec50, Vmax = coef(best.model)[grep("d:", names(coef(best.model)))],
        Km = coef(best.model)[grep("e:", names(coef(best.model)))],
        test = test.nm, model = best.model.nm
      ),
      fitFlag = TRUE, reliable = NULL, control = control
    )

  } else
  {
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = concgrid, fit.test = respgrid,
      model = best.model, parameters = list(
        EC50 = ec50, yEC50 = y.ec50, Vmax = NA,
        Km = NA, test = test.nm, model = best.model.nm
      ),
      fitFlag = TRUE, reliable = NULL, control = control
    )
  }

  class(drFitModel) <- "drFitModel"
  invisible(drFitModel)
}

#' Fit a biosensor model (Meyer et al., 2019) to response vs. concentration data
#'
#' @param flTable A dataframe containing the data for the dose-response model estimation. Such table of class \code{flTable} can be obtained by running \code{\link{flFit}} with \code{dr.method = 'model'} as argument in the \code{fl.control} object.
#' @param control A \code{fl.control} object created with \code{\link{fl.control}}, defining relevant fitting options.
#' @param dr.method (Character) Perform either a smooth spline fit on response parameter vs. concentration data (\code{'spline'}) or fit a biosensor response model with \code{'model'} (proposed by Meyer et al., 2019).
#' @param dr.parameter (Character or numeric) The response parameter in the output table to be used for creating a dose response curve. See \code{\link{fl.drFit}} for further details. Default: \code{'max_slope.spline'}, which represents the maximum slope of the spline fit Typical options include: \code{'max_slope.linfit'}, \code{'dY.linfit'}, \code{'max_slope.spline'}, and \code{'dY.spline'}.
#'
#' @return An object of class \code{drFit}.
#' \item{raw.data}{Data that passed to the function as \code{flTable}.}
#' \item{drTable}{Dataframe containing condition identifiers, fit options, and results of the dose-response analysis.}
#' \item{drFittedModels}{List of all \code{drFitModel} objects generated by the call of \code{\link{fl.drFitModel}} for each distinct experiment.}
#' \item{control}{Object of class \code{fl.control} created with the call of \code{\link{fl.control}}.}
#'
#' @export
#'
#' @details
#' Common response parameters used in dose-response analysis:<br><br><b>Linear fit:</b><br>- max_slope.linfit: Fluorescence increase rate<br>- lambda.linfit: Lag time<br>- dY.linfit: Maximum Fluorescence - Minimum Fluorescence<br>- A.linfit: Maximum fluorescence<br><br><b>Spline fit:</b><br>- max_slope.spline: Fluorescence increase rate<br>- lambda.spline: Lag time<br>- dY.spline: Maximum Fluorescence - Minimum Fluorescence<br>- A.spline: Maximum fluorescence<br>- integral.spline: Integral<br><br><b>Parametric fit:</b><br>- max_slope.model: Fluorescence increase rate<br>- lambda.model: Lag time<br>- dY.model: Maximum Fluorescence - Minimum Fluorescence<br>- A.model: Maximum fluorescence<br>- integral.model: Integral'
#'
#' @references Meyer, A.J., Segall-Shapiro, T.H., Glassey, E. et al. _Escherichia coli “Marionette” strains with 12 highly optimized small-molecule sensors._ Nat Chem Biol 15, 196–204 (2019). DOI: 10.1038/s41589-018-0168-3
#'
#' @examples
#' \donttest{
#' # Load example dataset
#' input <- read_data(data.fl = system.file('lac_promoters_fluorescence.txt', package = 'QurvE'),
#'                    csvsep.fl = "\t")
#'
#' # Run fluorescence curve analysis workflow
#' fitres <- flFit(fl_data = input$fluorescence,
#'                 time = input$time,
#'                 parallelize = FALSE,
#'                 control = fl.control(x_type = 'time', norm_fl = FALSE,
#'                                      suppress.messages = TRUE))
#'
#' # Perform dose-response analysis
#' drFit <- fl.drFit(flTable = fitres$flTable,
#'                   control = fl.control(dr.method = 'model',
#'                                        dr.parameter = 'max_slope.linfit'))
#'
#' # Inspect results
#' summary(drFit)
#' plot(drFit)
#' }
fl.drFit <- function(
    flTable, control = fl.control(
      dr.method = "model", dr.parameter = "max_slope.spline"
    )
)
{
  if (is(control) !=
      "fl.control")
    stop("control must be of class fl.control!")
  EC50.table <- NULL
  all.EC50 <- NA
  if (is.character(control$dr.parameter))
  {
    dr.parameter <- match(control$dr.parameter, colnames(flTable))
  }
  flTable <- flTable[!is.na(flTable[, 1]),
  ]
  table.tests <- table(
    (flTable[, 1])[which(
      (flTable[, 4] == TRUE) & (is.na(flTable[, dr.parameter]) ==
                                  FALSE)
    )]
  )
  distinct <- names(table.tests)
  EC50 <- list()
  validdata <- cbind(
    as.character(distinct),
    table.tests
  )
  colnames(validdata) <- c("TestID", "Number")
  rownames(validdata) <- rep("     ", times = dim(validdata)[1])
  if (control$suppress.messages == FALSE)
  {
    cat("\n")
    cat(
      "=== Dose-Response Estimation via Model Fit ==============================\n"
    )
    cat(
      "---------------------------------------------------\n"
    )
    cat("--> Checking data ...\n")
    cat(
      paste(
        "--> Number of distinct tests found:",
        as.character(length(distinct))
      ),
      "\n"
    )
    cat("--> Valid datasets per test: \n")
    print(validdata, quote = FALSE)
  }
  if (TRUE %in% (table.tests < control$dr.have.atleast))
  {
    warning(
      paste(
        "Warning: following tests have not enough ( <",
        as.character(control$dr.have.atleast - 1),
        ") datasets:\n"
      )
    )
    warning(
      paste(
        distinct[(table.tests < control$dr.have.atleast)],
        sep = "\n"
      )
    )
    warning("These tests will not be regarded\n")
    distinct <- distinct[table.tests >= control$dr.have.atleast]
  }
  if ((length(distinct)) ==
      0)
  {
    warning(
      paste(
        "There are no tests having enough ( >",
        as.character(control$dr.have.atleast - 1),
        ") datasets!\n"
      )
    )
    drFitfl <- list(
      raw.data = flTable, drTable = NA, drFittedModels = NA,
      control = control
    )
    class(drFitfl) <- "drFitfl"
    return(drFitfl)
  } else
  {
    skip <- c()
    for (i in 1:length(distinct))
    {
      conc <- factor(
        (flTable[, 3])[which(flTable[, 1] == distinct[i])]
      )
      if (length(levels(conc)) <
          4)
      {
        message(
          paste0(
            distinct[i], " does not have enough unique concentrations. A condition must have at least 4 different concentrations to be considered for dose-response analysis."
          )
        )
        skip <- c(skip, i)
        next
      }
      test <- (as.numeric(flTable[, dr.parameter]))[flTable[,
                                                            1] == distinct[i]]
      names(test) <- rep(
        names(flTable)[dr.parameter],
        length(test)
      )
      drID <- distinct[i]
      EC50[[i]] <- try(
        fl.drFitModel(conc, test, drID, control),
        silent = TRUE
      )
      description <- data.frame(
        Test = distinct[i], log.x = control$log.x.dr,
        log.y = control$log.y.dr
      )
      if (is(EC50[[i]]) !=
          "try-error")
      {
        out.row <- cbind(description, summary.drFitFLModel(EC50[[i]]))
      } else
      {
        out.row <- cbind(
          description, data.frame(
            yEC50 = NA, y.min = NA, y.max = NA,
            fc = NA, K = NA, n = NA, yEC50.orig = NA,
            K.orig = NA, test = NA
          )
        )
      }
      EC50.table <- rbind(EC50.table, out.row)

    }
  }
  class(EC50.table) <- c("drTable", "list")
  if (exists("skip") &&
      !is.null(skip))
  {
    distinct <- distinct[-skip]
    EC50 <- EC50[-skip]
  }
  names(EC50) <- distinct
  drFitfl <- list(
    raw.data = flTable, drTable = EC50.table, drFittedModels = EC50,
    control = control
  )
  class(drFitfl) <- "drFitfl"
  invisible(drFitfl)
}

#' Perform a biosensor model fit on response vs. concentration data of a single sample.
#'
#' \code{fl.drFitModel} fits the biosensor model proposed by Meyer et al. (2019) to the provided response (e.g., \code{max_slope.spline} vs. concentration data to determine the leakiness, sensitivity, induction fold-change, and cooperativity.
#'
#' @param conc Vector of concentration values.
#' @param test Vector of response parameter values of the same length as \code{conc}.
#' @param drID (Character) The name of the analyzed condition
#' @param control A \code{fl.control} object created with \code{\link{fl.control}}, defining relevant fitting options.
#'
#' @return A \code{drFitFLModel} object.
#' \item{raw.conc}{Raw data provided to the function as \code{conc}.}
#' \item{raw.test}{Raw data for the response parameter provided to the function as \code{test}.}
#' \item{drID}{(Character) Identifies the tested condition}
#' \item{fit.conc}{Fitted concentration values.}
#' \item{fit.test}{Fitted response values.}
#' \item{model}{\code{nls} object generated by the \code{\link{nlsLM}} function.}
#' \item{parameters}{List of parameters estimated from dose response curve fit.}
#' \itemize{
#' \item \code{yEC50}: Response value related to EC50.
#' \item \code{y.min}: Minimum fluorescence ('leakiness', if lowest concentration is 0).
#' \item \code{y.max}: Maximum fluorescence.
#' \item \code{fc}: Fold change (\code{y.max} divided by \code{y.min}).
#' \item \code{K}: Concentration at half-maximal response ('sensitivity').
#' \item \code{n}: Cooperativity.
#' \item \code{yEC50.orig}: Response value for EC50 in original scale, if a transformation was applied.
#' \item \code{K.orig}: K in original scale, if a transformation was applied.
#' \item \code{test.nm}: Test identifier extracted from \code{test}.
#' }
#' \item{fitFlag}{(Logical) Indicates whether a spline could fitted successfully to data.}
#' \item{reliable}{(Logical) Indicates whether the performed fit is reliable (to be set manually).}
#' \item{control}{Object of class \code{fl.control} created with the call of \code{\link{fl.control}}.}
#' Use \code{\link{plot.drFitModel}} to visualize the model fit.
#'
#' @export
#' @importFrom minpack.lm nlsLM
#'
#' @references Meyer, A.J., Segall-Shapiro, T.H., Glassey, E. et al. _Escherichia coli “Marionette” strains with 12 highly optimized small-molecule sensors._ Nat Chem Biol 15, 196–204 (2019). DOI: 10.1038/s41589-018-0168-3
#'
#' @examples
#' # Create concentration values via a serial dilution
#' conc <- c(0, rev(unlist(lapply(1:18, function(x) 10*(2/3)^x))),10)
#'
#' # Simulate response values via biosensor equation
#' response <- biosensor.eq(conc, y.min = 110, y.max = 6000, K = 0.5, n = 2) +
#'             0.01*6000*rnorm(10)
#'
#' # Perform fit
#' TestRun <- fl.drFitModel(conc, response, drID = 'test', control = fl.control())
#'
#' print(summary(TestRun))
#' plot(TestRun)
fl.drFitModel <- function(
    conc, test, drID = "undefined", control = fl.control()
)
{
  if (is(control) !=
      "fl.control")
    stop("control must be of class fl.control!")
  test.nm <- names(test)[1]
  test <- as.vector(as.numeric(as.matrix(test)))
  conc <- as.vector(as.numeric(as.matrix(conc)))
  if (is.vector(conc) ==
      FALSE || is.vector(test) ==
      FALSE)
    stop(
      "fl.drFitModel: dose or response data must be a vector !"
    )
  if (control$neg.nan.act == FALSE)
  {
    missings <- is.na(conc) |
      is.na(test) |
      !is.numeric(conc) |
      !is.numeric(test)
    conc <- conc[!missings]
    test <- test[!missings]
    negs <- (conc < 0) | (test < 0)
    conc <- conc[!negs]
    test <- test[!negs]
  } else
  {
    if (sum(
      is.na(conc) |
      is.na(test)
    ))
      stop(
        "fl.drFitModel: NA values encountered. Program terminated"
      )
    if ((sum((conc < 0)) >
         0) | (sum((test < 0)) >
               0))
      stop(
        "fl.drFitModel: Negative values encountered. Program terminated"
      )
    if ((FALSE %in% is.numeric(conc)) ||
        (FALSE %in% is.numeric(test)))
      stop(
        "fl.drFitModel: Non numeric values encountered. Program terminated"
      )
  }
  if (length(test) <
      4)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drFitModel: There is not enough valid data. Must have at least 4 unique values!"
      )
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, EC50.orig = NA,
        yEC50.orig = NA
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitFLModel"
    return(drFitModel)
  }
  if (length(test) <
      control$dr.have.atleast)
  {
    if (control$suppress.messages == FALSE)
      message(
        "drFitModel: number of valid data points is below the number specified in 'dr.have.atleast'. See fl.control()."
      )
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, spline = NA,
      parameters = list(
        EC50 = NA, yEC50 = NA, EC50.orig = NA,
        yEC50.orig = NA
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitFLModel"
    return(drFitModel)
  }
  if (control$log.x.dr == TRUE)
    conc.log <- log(conc + 1)
  if (control$log.y.dr == TRUE)
    test.log <- log(test + 1)
  if (control$log.x.dr == TRUE)
  {
    conc.fit <- log(conc + 1)
  } else
  {
    conc.fit <- conc
  }
  if (control$log.y.dr == TRUE)
  {
    test.fit <- log(test + 1)
  } else
  {
    test.fit <- test
  }
  test.fit <- test.fit[order(conc.fit)]
  conc.fit <- conc.fit[order(conc.fit)]
  y.min <- mean(test.fit[which(conc.fit == conc.fit[1])])
  fitFlag <- TRUE
  n_candidates <- seq(0, 1, 0.01)
  i <- 1
  # plot(conc.fit, test.fit) title(drID)
  df <- data.frame(x = conc.fit, test.fit = test.fit)
  y.model <- list()
  y.model[["convInfo"]] <- list()
  y.model$convInfo[["isConv"]] <- FALSE
  while (y.model$convInfo$isConv == FALSE && i <
         100)
  {
    i <- i + 1
    try(
      suppressWarnings(
        y.model <- minpack.lm::nlsLM(
          test.fit ~ biosensor.eq(
            x = conc.fit, y.min, y.max, K,
            n
          ),
          start = initbiosensor(x = conc.fit, y = test.fit, n = n_candidates[i])
        )
      ),
      silent = T
    )
    # try( suppressWarnings(y.model <-
    # nls(test.fit ~ biosensor.eq(x=conc.fit,
    # y.min, y.max, K, n), start =
    # initbiosensor(x=conc.fit, y=test.fit, n
    # = n_candidates[i]))), silent = F )
  }
  if (y.model$convInfo$isConv == FALSE)
  {
    if (control$suppress.messages == FALSE)
    {
      warning(
        "Model could not be fitted in dose-response analysis!\n"
      )
    }
    drFitModel <- list(
      raw.conc = conc, raw.test = test, drID = drID,
      fit.conc = NA, fit.test = NA, model = y.model,
      parameters = list(
        yEC50 = NA, y.min = NA, y.max = NA,
        fc = NA, K = NA, n = NA, yEC50.orig = NA,
        K.orig = NA, test = NA
      ),
      fitFlag = FALSE, reliable = NULL, control = control
    )
    class(drFitModel) <- "drFitFLModel"
    return(drFitModel)
  }
  # lines(conc.fit, biosensor.eq(x=conc.fit,
  # y.min=coef(y.model)[1],
  # y.max=coef(y.model)[2], K=coef(y.model)[3],
  # n=coef(y.model)[4]), col = 'red')
  m <- summary(y.model)
  par <- m$parameters
  # y.min <- par[1,1]
  y.max <- par[1, 1]
  fc <- par[1, 1]/y.min  # fold-change
  K <- par[2, 1]  # sensitivity
  n <- par[3, 1]  # cooperativity

  x_fit <- seq(
    0, max(conc.fit),
    length.out = 1000
  )
  y_fit <- biosensor.eq(
    x_fit, y.min = y.min, y.max = par[1, 1], K = par[2,
                                                     1], n = par[3, 1]
  )
  # plot(conc.fit, test.fit) lines(xin,
  # biosensor.eq(xin, y.min=par[1,1],
  # y.max=par[2,1], K=par[3,1], n=par[4,1]))

  # /// estimating EC 50 values
  yEC.test <- y.min + (y.max - y.min) * (K^n/(K^n +
                                                K^n))
  EC.test <- K
  if (control$suppress.messages == FALSE)
  {
    cat(
      "\n\n=== Dose response curve estimation ================\n"
    )
    cat(
      "--- EC 50 -----------------------------------------\n"
    )
    cat(paste("-->", as.character(drID)))
    cat("\n")
    cat(
      paste(
        c(
          "sensitivity:", "yEC50:", "fold change:",
          "leakiness:"
        ),
        c(
          signif(EC.test, digits = 3),
          round(yEC.test),
          round(fc, 2),
          round(y.min, 1)
        ),
        collapse = " | "
      )
    )
  }
  if ((control$log.x.dr == TRUE) && (control$log.y.dr ==
                                     FALSE))
  {
    if (control$suppress.messages == FALSE)
    {
      cat("\n--> Original scale \n")
      cat(
        paste(
          c("xEC50", "yEC50"),
          c(
            exp(EC.test) -
              1, yEC.test
          )
        )
      )
    }
    EC.orig <- c(
      exp(EC.test) -
        1, yEC.test
    )
  } else
  {
    if ((control$log.x.dr == FALSE) && (control$log.y.dr ==
                                        TRUE))
    {
      if (control$suppress.messages == FALSE)
      {
        cat("\n--> Original scale \n")
        cat(
          paste(
            c("xEC50", "yEC50"),
            c(
              EC.test, exp(yEC.test) -
                1
            )
          )
        )
      }
      EC.orig <- c(
        EC.test, exp(yEC.test) -
          1
      )
    } else
    {
      if ((control$log.x.dr == TRUE) && (control$log.y.dr ==
                                         TRUE))
      {
        if (control$suppress.messages == FALSE)
        {
          cat("\n--> Original scale \n")
          cat(
            paste(
              c("xEC50", "yEC50"),
              c(
                exp(EC.test) -
                  1, exp(yEC.test) -
                  1
              )
            )
          )
        }
        EC.orig <- c(
          exp(EC.test) -
            1, exp(yEC.test) -
            1
        )
      } else
      {
        if ((control$log.x.dr == FALSE) &&
            (control$log.y.dr == FALSE))
        {
          EC.orig <- c(EC.test, yEC.test)
        }
      }
    }
  }
  if (control$suppress.messages == FALSE)
  {
    cat("\n\n\n")
  }
  drFitModel <- list(
    raw.conc = conc, raw.test = test, drID = drID,
    fit.conc = x_fit, fit.test = y_fit, model = y.model,
    parameters = list(
      yEC50 = yEC.test, y.min = y.min, y.max = y.max,
      fc = fc, K = K, n = n, yEC50.orig = EC.orig[2],
      K.orig = EC.orig[1], test = test.nm
    ),
    fitFlag = fitFlag, reliable = NULL, control = control
  )
  class(drFitModel) <- "drFitFLModel"
  invisible(drFitModel)
}