R/multipleFitFunction_h0.R

In sicegar: Analysis of Single-Cell Viral Growth Curves

#' @title multiple fit function h0.
#'
#' @param dataInput A dataframe or a list containing the dataframe. The dataframe should be composed of at least two columns. One represents time, and the other represents intensity. The data should be normalized with the normalize data function sicegar::normalizeData() before being imported into this function.
#' @param model Type of fit model that will be used. Can be "sigmoidal", or "double_sigmoidal".
#' @param n_runs_max This number indicates the upper limit of the fitting attempts. Default is 500.
#' @param n_runs_min This number indicates the lower limit of the successful fitting attempts. It should be smaller than the upper limit of the fitting attempts (n_runs_max). Default is 20.
#' @param showDetails Logical if TRUE prints details of intermediate steps of individual fits (Default is FALSE).
#' @param dataInputName Name of data set (Default is 'NA').
#' @param ... All other arguments that model functions ("sigmoidalFitFunction" and, "doublesigmoidalFitFunction") may need.
#'
#' @description Calls the fitting algorithms to fit the data multiple times with starting from different randomly generated initial parameters in each run. Multiple attempts at fitting the data are necessary to avoid local minima.
#' @return Returns the parameters related with the model fitted for the input data.
#' @export
#'
#' @examples
#' # runif() is used here for consistency with previous versions of the sicegar package. However,
#' # rnorm() will generate symmetric errors, producing less biased numerical parameter estimates.
#' # We recommend errors generated with rnorm() for any simulation studies on sicegar.
#'# Example 1 (sigmoidal function with normalization)

#'time <- seq(3, 24, 0.5)
#'
#'#simulate intensity data and add noise
#'noise_parameter <- 2.5
#'intensity_noise <- stats::runif(n = length(time), min = 0, max = 1) * noise_parameter
#'intensity <- sigmoidalFitFormula_h0(time, maximum = 4, slopeParam = 1, midPoint = 8, h0 = 1)
#'intensity <- intensity + intensity_noise
#'
#'dataInput <- data.frame(intensity = intensity, time = time)
#'normalizedInput <- normalizeData(dataInput, dataInputName = "sample001")
#'parameterVector <- multipleFitFunction_h0(dataInput = normalizedInput,
#'                                       model = "sigmoidal",
#'                                       n_runs_min = 20,
#'                                       n_runs_max = 500)
#'
#'#Check the results
#'if(parameterVector$isThisaFit){
#'   intensityTheoretical <- sigmoidalFitFormula_h0(time,
#'                             maximum = parameterVector$maximum_Estimate,
#'                             slopeParam = parameterVector$slopeParam_Estimate,
#'                             midPoint = parameterVector$midPoint_Estimate,
#'                             h0 = parameterVector$h0_Estimate)
#'
#'  comparisonData <- cbind(dataInput, intensityTheoretical)
#'
#'  print(parameterVector$residual_Sum_of_Squares)
#'
#'  require(ggplot2)
#'  ggplot(comparisonData)+
#'    geom_point(aes(x = time, y = intensity)) +
#'    geom_line(aes(x = time, y = intensityTheoretical), color = "orange") +
#'    expand_limits(x = 0, y = 0)
#' }
#'
#'
#'
#'if(!parameterVector$isThisaFit){
#'   print(parameterVector)
#'}
#'

#'# Example 2 (doublesigmoidal function with normalization)
#'time <- seq(3, 24, 0.1)
#'
#'#simulate intensity data with noise
#'noise_parameter <- 0.2
#'intensity_noise <- stats::runif(n = length(time), min = 0, max = 1) * noise_parameter
#'intensity <- doublesigmoidalFitFormula_h0(time,
#'                                       finalAsymptoteIntensityRatio = .3,
#'                                       maximum = 4,
#'                                       slope1Param = 1,
#'                                       midPoint1Param = 7,
#'                                       slope2Param = 1,
#'                                       midPointDistanceParam = 8,
#'                                       h0 = 1)
#'intensity <- intensity + intensity_noise
#'
#'dataInput <- data.frame(intensity = intensity, time = time)
#'normalizedInput <- normalizeData(dataInput)
#'parameterVector <- multipleFitFunction_h0(dataInput = normalizedInput,
#'                            dataInputName="sample001",
#'                            model = "doublesigmoidal",
#'                            n_runs_min = 20,
#'                            n_runs_max = 500,
#'                            showDetails = FALSE)
#'
#'
#'#Check the results
#'if(parameterVector$isThisaFit){
#'  intensityTheoretical <-
#'        doublesigmoidalFitFormula_h0(
#'                time,
#'                finalAsymptoteIntensityRatio = parameterVector$finalAsymptoteIntensityRatio_Estimate,
#'                maximum = parameterVector$maximum_Estimate,
#'                slope1Param = parameterVector$slope1Param_Estimate,
#'                midPoint1Param = parameterVector$midPoint1Param_Estimate,
#'                slope2Param = parameterVector$slope2Param_Estimate,
#'                midPointDistanceParam = parameterVector$midPointDistanceParam_Estimate,
#'                h0 = parameterVector$h0_Estimate)
#'
#'  comparisonData <- cbind(dataInput, intensityTheoretical)
#'
#'  require(ggplot2)
#'  ggplot(comparisonData) +
#'    geom_point(aes(x = time, y = intensity)) +
#'    geom_line(aes(x = time, y = intensityTheoretical), color = "orange") +
#'    expand_limits(x = 0, y = 0)
#'  }
#'
#'if(!parameterVector$isThisaFit){
#'   print(parameterVector)
#'   }
#'
#'
multipleFitFunction_h0 <-
  function(dataInput,
           dataInputName = NA,
           model,
           n_runs_min = 20,
           n_runs_max = 500,
           showDetails = FALSE, ...)
  {
    dataInputCheck <- dataCheck(dataInput)

    if(!(model %in% c("sigmoidal", "doublesigmoidal")) ){
      stop("model should be one of sigmoidal, doublesigmoidal")
    }

    counterBetterFit <- 0
    counterCorrectFit <- 0
    counterTotalFit <- 0
    residual_Sum_of_Squares_min <- Inf
    storedModelOutput <- c()
    storedModelOutput$residual_Sum_of_Squares <- Inf

    while(counterCorrectFit < n_runs_min & counterTotalFit < n_runs_max){

      counterTotalFit <- counterTotalFit+1
      if(model == "sigmoidal"){
        modelOutput <- sigmoidalFitFunction_h0(dataInput = dataInput,
                                            tryCounter = counterTotalFit, ...)
      }

      if(model == "doublesigmoidal"){
        modelOutput <- doublesigmoidalFitFunction_h0(dataInput = dataInput,
                                                  tryCounter = counterTotalFit, ...)
      }

      if(is.na(dataInputName)){

        isalist <- (is.list(dataInput) & !is.data.frame(dataInput))
        if(isalist){
          modelOutput$dataInputName <- as.character(dataInput$dataInputName)
        }

        if(!isalist){
          modelOutput$dataInputName <- NA
        }
      }

      if(!is.na(dataInputName)){

        isalist <- (is.list(dataInput) & !is.data.frame(dataInput))
        if(isalist){

          if(is.na(dataInput$dataInputName)){

            modelOutput$dataInputName <- as.character(dataInputName)
          }

          if(!is.na(dataInput$dataInputName))
          {
            if(dataInput$dataInputName != dataInputName){
              stop("the input data has already have a name")
            }
            if(dataInput$dataInputName == dataInputName){
              modelOutput$dataInputName <- as.character(dataInputName)
            }
          }
        }
        if(!isalist){
          modelOutput$dataInputName <- as.character(dataInputName)
        }
      }


      if(modelOutput[["isThisaFit"]]){
        counterCorrectFit <- counterCorrectFit + 1

        if(residual_Sum_of_Squares_min > modelOutput$residual_Sum_of_Squares){
          counterBetterFit <- counterBetterFit + 1
          residual_Sum_of_Squares_min <- modelOutput$residual_Sum_of_Squares
          storedModelOutput <- modelOutput
        }
      }

      if(showDetails){
        print(c(betterFit = counterBetterFit,
                correctFit = counterCorrectFit,
                totalFit = counterTotalFit,
                currentOutput = modelOutput$residual_Sum_of_Squares,
                bestOutput = storedModelOutput$residual_Sum_of_Squares))}

    }

    # add number of independent runs to outputs
    # might be important for quality checks

    storedModelOutput$betterFit <- counterBetterFit
    storedModelOutput$correctFit <- counterCorrectFit
    storedModelOutput$totalFit <- counterTotalFit

    return(storedModelOutput)
  }