R/doublesigmoidalFitFunctions_h0.R

Defines functions doublesigmoidalRenormalizeParameters_h0 f_slope_doublesigmoidal_h0 f_argmax_doublesigmoidal_h0 doublesigmoidalFitFormula_h0 doublesigmoidalFitFunction_h0

Documented in doublesigmoidalFitFormula_h0 doublesigmoidalFitFunction_h0

#' @title Double sigmoidal fit function with 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 tryCounter A counter that shows the number of times the data was fit via maximum likelihood function.
#' @param startList The vector containing the initial set of parameters that algorithm tries for the first fit attempt for the relevant parameters. The vector composes of seven elements; 'finalAsymptoteIntensityRatio', 'maximum', 'slope1Param', 'midPoint1Param' , 'slope2Param', 'midPointDistanceParam', and 'h0'. Detailed explanations of those parameters can be found in vignettes. Defaults are  finalAsymptoteIntensityRatio = 0, maximum = 1, slope1Param = 1, midPoint1Param = 0.33, slope2Param = 1, midPointDistanceParam=0.29, and h0 = 0. The numbers are in normalized time intensity scale.
#' @param lowerBounds The lower bounds for the randomly generated start parameters. The vector is composed of seven elements; 'finalAsymptoteIntensityRatio', 'maximum', 'slope1Param', 'midPoint1Param' , 'slope2Param', 'midPointDistanceParam', and 'h0'. Detailed explanations of those parameters can be found in vignettes. Defaults are finalAsymptoteIntensityRatio = 0, maximum = 0.3, slope1Param = .01, midPoint1Param = -0.52, slope2Param = .01, midPointDistanceParam = 0.04, and h0 = -0.1. The numbers are in normalized time intensity scale.
#' @param upperBounds The upper bounds for the randomly generated start parameters. The vector is composed of seven elements; 'finalAsymptoteIntensityRatio', 'maximum', 'slope1Param', 'midPoint1Param' , 'slope2Param', 'midPointDistanceParam', and 'h0'. Detailed explanations of those parameters can be found in vignettes. Defaults are finalAsymptoteIntensityRatio = 1, maximum = 1.5, slope1Param = 180, midPoint1Param = 1.15, slope2Param = 180, midPointDistanceParam = 0.63, and h0 = 0.3. The numbers are in normalized time intensity scale.
#' @param min_Factor Defines the minimum step size used by the fitting algorithm. Default is 1/2^20.
#' @param n_iterations Defines maximum number of iterations used by the fitting algorithm. Default is 1000.
#'
#' @description The function fits a double sigmoidal curve to given data by using likelihood maximization (LM) algorithm and provides the parameters (maximum, final asymptote intensity, slope1Param, midpoint1Param, slope2Param, and midpoint distance) describing the double sigmoidal fit as output. It also contains information about the goodness of fits such as AIC, BIC, residual sum of squares, and log likelihood.
#' @return Returns the fitted parameters and goodness of fit metrics.
#' @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.
#'time=seq(3, 24, 0.1)
#'
#'#simulate intensity data and add 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 <- doublesigmoidalFitFunction_h0(normalizedInput, tryCounter = 1)
#'
#'
#'#Check the results
#'# doublesigmoidalFitFunction_h0() is run on the startList param values (because 'tryCounter = 1')
#'# use multipleFitFunction() for multiple random starts in order to optimize
#'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)) +
#'    expand_limits(x = 0, y = 0)}
#'
#'if(!parameterVector$isThisaFit) {print(parameterVector)}
#'
doublesigmoidalFitFunction_h0 <- function(dataInput,
                                       tryCounter,
                                       startList=list(finalAsymptoteIntensityRatio = 0,
                                                      maximum = 1,
                                                      slope1Param = 1,
                                                      midPoint1Param = 0.33,
                                                      slope2Param = 1,
                                                      midPointDistanceParam = 0.29,
                                                      h0 = 0.1),
                                       lowerBounds=c(finalAsymptoteIntensityRatio = 0,
                                                     maximum = 0.3,
                                                     slope1Param = .01,
                                                     midPoint1Param = -0.52,
                                                     slope2Param = .01,
                                                     midPointDistanceParam = 0.04,
                                                     h0 =  0),
                                       upperBounds=c(finalAsymptoteIntensityRatio = 1,
                                                     maximum = 1.5,
                                                     slope1Param = 180,
                                                     midPoint1Param = 1.15,
                                                     slope2Param = 180,
                                                     midPointDistanceParam = 0.63,
                                                     h0 = 0.3),
                                       min_Factor = 1/2^20,
                                       n_iterations = 1000)
{

  isalist = (is.list(dataInput) & !is.data.frame(dataInput))
  if(isalist){
    dataFrameInput <- dataInput$timeIntensityData
  }
  isadataframe <- (is.data.frame(dataInput))
  if(isadataframe){
    dataFrameInput <- dataInput
  }

  if(tryCounter==1){
    counterDependentStartList <- startList
  }
  if(tryCounter!=1){
    randomVector <- stats::runif(length(startList), 0, 1)
    names(randomVector) <- c("finalAsymptoteIntensityRatio",
                             "maximum",
                             "slope1Param",
                             "midPoint1Param",
                             "slope2Param",
                             "midPointDistanceParam",
                             "h0")
    counterDependentStartVector <- randomVector * (upperBounds - lowerBounds) + lowerBounds
    counterDependentStartList <- as.list(counterDependentStartVector)
  }


  theFitResult <- try(minpack.lm::nlsLM(intensity ~ doublesigmoidalFitFormula_h0(time,
                                                                                       finalAsymptoteIntensityRatio,
                                                                                       maximum,
                                                                                       slope1Param,
                                                                                       midPoint1Param,
                                                                                       slope2Param,
                                                                                       midPointDistanceParam,
                                                                                       h0),
                                        dataFrameInput,
                                        start = counterDependentStartList,
                                        control = list(maxiter = n_iterations, minFactor = min_Factor),
                                        lower = lowerBounds,
                                        upper = upperBounds,
                                        trace=F), silent = TRUE)

  if(!inherits(theFitResult, "try-error"))
  {
    parameterMatrix <- summary(theFitResult)$parameters
    colnames(parameterMatrix) <- c("Estimate", "Std_Error", "t_value", "Pr_t")

    parameterVector <- c(t(parameterMatrix))
    names(parameterVector) <- c("finalAsymptoteIntensityRatio_N_Estimate", "finalAsymptoteIntensityRatio_Std_Error", "finalAsymptoteIntensityRatio_t_value", "finalAsymptoteIntensityRatio_Pr_t",
                                "maximum_N_Estimate", "maximum_Std_Error", "maximum_t_value", "maximum_Pr_t",
                                "slope1Param_N_Estimate", "slope1Param_Std_Error", "slope1Param_t_value", "slope1Param_Pr_t",
                                "midPoint1Param_N_Estimate", "midPoint1Param_Std_Error", "midPoint1Param_t_value", "midPoint1Param_Pr_t",
                                "slope2Param_N_Estimate", "slope2Param_Std_Error", "slope2Param_t_value", "slope2Param_Pr_t",
                                "midPointDistanceParam_N_Estimate", "midPointDistanceParam_Std_Error", "midPointDistanceParam_t_value", "midPointDistanceParam_Pr_t",
                                "h0_N_Estimate","h0_Std_Error","h0_t_value","h0_Pr_t"   # ← appended
                               )

    parameterVector <- c(parameterVector,
                         residual_Sum_of_Squares = sum((as.vector(stats::resid(theFitResult)))^2)[1],
                         log_likelihood = as.vector(stats::logLik(theFitResult))[1],
                         AIC_value = as.vector(stats::AIC(theFitResult))[1],
                         BIC_value = as.vector(stats::BIC(theFitResult))[1])

    parameterList <- as.list(parameterVector)
    parameterList$isThisaFit <- TRUE
    parameterList$startVector <- counterDependentStartList
    if(isalist){
      parameterList$dataScalingParameters <- as.list(dataInput$dataScalingParameters)
    }
    parameterList$model <- as.character("doublesigmoidal")
    parameterList$additionalParameters <- FALSE

    parameterDf <- as.data.frame(parameterList)
    #Renormalize Parameters
    parameterDf <- doublesigmoidalRenormalizeParameters_h0(parameterDf, isalist)

  }

  if(inherits(theFitResult, "try-error"))
  {
    parameterVector <- rep(NA, 28)
    names(parameterVector) <- c("finalAsymptoteIntensityRatio_N_Estimate", "finalAsymptoteIntensityRatio_Std_Error", "finalAsymptoteIntensityRatio_t_value", "finalAsymptoteIntensityRatio_Pr_t",
                                "maximum_N_Estimate", "maximum_Std_Error", "maximum_t_value", "maximum_Pr_t",
                                "slope1Param_N_Estimate", "slope1Param_Std_Error", "slope1Param_t_value", "slope1Param_Pr_t",
                                "midPoint1Param_N_Estimate", "midPoint1Param_Std_Error", "midPoint1Param_t_value", "midPoint1Param_Pr_t",
                                "slope2Param_N_Estimate", "slope2Param_Std_Error", "slope2Param_t_value", "slope2Param_Pr_t",
                                "midPointDistanceParam_N_Estimate", "midPointDistanceParam_Std_Error", "midPointDistanceParam_t_value", "midPointDistanceParam_Pr_t",
                                "h0_N_Estimate","h0_Std_Error","h0_t_value","h0_Pr_t"   # ← appended
                                )

    parameterVector <- c(parameterVector,
                         residual_Sum_of_Squares = Inf,
                         log_likelihood = NA,
                         AIC_value = NA,
                         BIC_value = NA)

    parameterList <- as.list(parameterVector)
    parameterList$isThisaFit <- FALSE
    parameterList$startVector <- counterDependentStartList
    if(isalist) {
      parameterList$dataScalingParameters <- as.list(dataInput$dataScalingParameters)
    }
    parameterList$model <- "doublesigmoidal"

    parameterDf <- as.data.frame(parameterList)
    #Renormalize Parameters
    parameterDf <- doublesigmoidalRenormalizeParameters_h0(parameterDf, isalist)

  }

  return(parameterDf)

}
#**************************************


#**************************************
#' @title Double Sigmoidal Formula with h0
#'
#' @param x the "time" column of the dataframe
#' @param finalAsymptoteIntensityRatio this is the ratio between the final asymptote intensity and maximum intensity of the fitted curve.
#' @param maximum the maximum intensity that the double sigmoidal function reaches.
#' @param slope1Param the slope parameter of the sigmoidal function at the steepest point in the exponential phase of the viral production.
#' @param midPoint1Param the x axis value of the steepest point in the function.
#' @param slope2Param the slope parameter of the sigmoidal function at the steepest point in the lysis phase. i.e when the intensity is decreasing.
#' @param midPointDistanceParam the distance between the time of steepest increase and steepest decrease in the intensity data. In other words the distance between the x axis values of arguments of slope1Param and slope2Param.
#' @param h0 the lower asymptote (baseline) intensity
#'
#' @description Calculates intensities using the double sigmoidal model fit and the parameters (maximum, final asymptote intensity, slope1Param, midpoint1Param, slope2Param, and midpoint distance).
#' @return Returns the predicted intensities for the given time points with the double sigmoidal fitted parameters for the double sigmoidal fit.
#' @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.
#'time <- seq(3, 24, 0.1)
#'
#'#simulate intensity data and add 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 <- doublesigmoidalFitFunction_h0(normalizedInput, tryCounter = 1)
#'
#'
#'#Check the results
#'# doublesigmoidalFitFunction_h0() is run on the startList param values (because 'tryCounter = 1')
#'# use multipleFitFunction() for multiple random starts in order to optimize
#'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)) +
#'    expand_limits(x = 0, y = 0)
#'    }
#'
#'if(!parameterVector$isThisaFit){
#'   print(parameterVector)
#'   }
#'
#'
doublesigmoidalFitFormula_h0 <- function(x,
                                      finalAsymptoteIntensityRatio,
                                      maximum,
                                      slope1Param,
                                      midPoint1Param,
                                      slope2Param,
                                      midPointDistanceParam,
                                      h0)
{
  if(slope1Param < 0){
    stop("slope1Param should be a positive number")
  }
  if(slope2Param < 0){
    stop("slope2Param should be a positive number. It is the absolute value of the second slopeParam")
  }
  if(midPointDistanceParam < 0){
    stop("midPointDistanceParam should be a positive number. It is the distance between two steepest points of exponential phase and lysis")
  }
  if(finalAsymptoteIntensityRatio < 0 | finalAsymptoteIntensityRatio > 1){
    stop("finalAsymptoteIntensityRatio should be a number between 0 and 1")
  }
  if(maximum < 0){
    stop("maximum should be a positive number")
  }

  optimizeIntervalMin <- midPoint1Param - 2 * midPointDistanceParam
  optimizeIntervalMax <- midPoint1Param + 3 * midPointDistanceParam
  xmax <- stats::optimize(f1,
                          c(optimizeIntervalMin, optimizeIntervalMax),
                          tol = 0.0001,
                          B1 = slope1Param, M1 = midPoint1Param, B2 = slope2Param, L = midPointDistanceParam, maximum = TRUE);
  argumentt <- xmax$maximum;
  constt <- f0(argumentt, slope1Param, midPoint1Param, slope2Param, midPointDistanceParam);
  y <- f2_h0(x,
          finalAsymptoteIntensityRatio, maximum,
          slope1Param, midPoint1Param,
          slope2Param, midPointDistanceParam,
          constt, argumentt, h0);
  return(y)
}
#**************************************


#**************************************
# Those four functions are necessary internal steps for double sigmoidal model
# All four are internal functions
# Name convention is different than the rest of the package
#     A2 correspond to finalAsymptoteIntensityRatio
#     Ka correspond to maximum
#     B1 correspond to slope1Param
#     B2 correspond to slope2Param
#     M1 correspond to midpoint1Param
#     L correspond to midPointDistanceParam
#     argument correspond to maximum_x (i.e the time when intensity reaches maximum with
#                                       respect to the double sigmoidal model)
#     const correspond to the maximum value the function reaches when the input data intensity and
#           is normalized between 0,1

# f0, f1, and f0mid are unchanged and are available in doublesigmoidlFitFunctions.R

#     f0 <- function (x, B1, M1, B2, L) {
#       (1/( ( 1+exp(-B1 * (x - M1)) ) * ( 1 + exp(B2 * (x - (M1 + L))) ) ))
#     }
#     f1 <- function (x, B1, M1, B2, L) {
#       log((1 / ( ( 1 + exp(-B1 * (x - M1)) ) * ( 1 + exp(B2 * (x - (M1 + L))) ) )))
#     }
f2_h0 <- function (x, A2, Ka, B1, M1, B2, L, const, argument, h0){
  fBasics::Heaviside(x - argument) * (f0(x, B1, M1, B2, L) * ((Ka - A2 * Ka)/(const)) + A2 * Ka) +
    (1 - fBasics::Heaviside(x - argument)) * (f0(x, B1, M1, B2, L) * ((Ka - h0)/(const)) + h0)
}
#     f0mid <- function (x, B1, M1, B2, L,const) {
#       -const / 2 + 1 / ( ( 1 + exp(-B1 * (x-M1)) ) * ( 1 + exp(B2*(x-(M1+L))) ) )
#     }
#**************************************


#**************************************
# @title f_argmax_doublesigmoidal_h0
#
# @param parameterVector is a dataFrame that includes columns with names:
#
# *"dataScalingParameters.timeRange" that represents maximum time of intensity measurements,
#
# parameters of double sigmoidal model:
#  *"slope1Param_Estimate",
#  *"midPoint1Param_Estimate",
#  *"slope2Param_Estimate",
#  *"midPointDistanceParam_Estimate".
#  (Note: additional parameters in Df does not cause problems)
#
# finally a column with the used model information model="doublesigmoidal".
# @description This function is designed to compensate for the small discrepancy in numerical parameters in the double sigmoidal model. This function is called by numericalReCalculation.
# @return return the x value of the maximum in between time 0 and maximum time measured in sample intensities. In other words the time that the double sigmoidal function reaches its maximum.
#
# @examples
# # Related example
#
## Initial Command to Reset the System
#rm(list = ls())
#if (is.integer(dev.list())){dev.off()}
#cat("\014")
#
## Generate a parameters DF
#parameterDf = data.frame(dataScalingParameters.timeRange = 24,
#                         slope1Param_Estimate = 1,
#                         midPoint1Param_Estimate = 7,
#                         slope2Param_Estimate = 1,
#                         midPointDistanceParam_Estimate = 8,
#                         model = "doublesigmoidal")
#
# # Find the x value of the maximum for this model
# # with in the range of [0, dataScalingParameters.timeRange]
# maximum_x = f_argmax_doublesigmoidal_h0(parameterDf)
#
f_argmax_doublesigmoidal_h0 <- function(parameterVector){

  slope1Param <- parameterVector$slope1Param_N_Estimate
  slope2Param <- parameterVector$slope2Param_N_Estimate
  midPointDistanceParam <- parameterVector$midPointDistanceParam_N_Estimate
  finalAsymptoteIntensityRatio <- parameterVector$finalAsymptoteIntensityRatio_N_Estimate
  maximum <- parameterVector$maximum_N_Estimate
  midPoint1Param <- parameterVector$midPoint1Param_N_Estimate
  timeRange <- parameterVector$dataScalingParameters.timeRange

  if(slope1Param <0 ){
    stop("slope1Param should be a positive number")
  }
  if(slope2Param < 0){
    stop("slope2Param should be a positive number. It is the absolute value of the slope2Param")
  }
  if(midPointDistanceParam < 0){
    stop("midPointDistanceParam should be a positive number. It is the distance between two steepest points of exponential phase and lysis")
  }
  if(finalAsymptoteIntensityRatio < 0 | finalAsymptoteIntensityRatio > 1){
    stop("finalAsymptoteIntensityRatio should be a number between 0 and 1")
  }
  if(maximum < 0){
    stop("maximum should be a positive number")
  }

  optimizeIntervalMin <- midPoint1Param - 2 * midPointDistanceParam
  optimizeIntervalMax <- midPoint1Param + 3 * midPointDistanceParam
  xmax <- stats::optimize(f1,
                          c(optimizeIntervalMin, optimizeIntervalMax),
                          tol = 0.0001,
                          B1 = slope1Param, M1 = midPoint1Param, B2 = slope2Param, L = midPointDistanceParam, maximum = TRUE);
  argumentt <- xmax$maximum * timeRange;
  return(argumentt)
}
#**************************************


#**************************************
# f_mid1_doublesigmoidal is unchanged and is available in doublesigmoidalFitFunctions.R

#     # @title f_mid1_doublesigmoidal
#     #
#     # @param parameterDf is a dataFrame that includes columns with names:
#     #
#     # *"dataScalingParameters.timeRange" that represents maximum time of intensity measurements,
#     #
#     # parameters of double sigmoidal model:
#     #  *"slope1Param_Estimate",
#     #  *"midPoint1Param_Estimate",
#     #  *"slope2Param_Estimate",
#     #  *"midPointDistanceParam_Estimate".
#     #  (Note: additional parameters in Df does not cause problems)
#     #
#     # finally a column with the used model information model="doublesigmoidal".
#     #
#     # @description  This function is designed to compensate for the small discrepancy the in numerical parameters in double sigmoidal model. This function is called by numericalReCalculation.
#     # @return returns the x values of midpoint1Param, by assuming the y of midpoint1Param is the half of the maximum value it reaches.
#     # @export
#     #
#     # @examples
#     # # Related example
#     #
#     ## Initial Command to Reset the System
#     #rm(list = ls())
#     #if (is.integer(dev.list())){dev.off()}
#     #cat("\014")
#     #
#     ## Generate a parameters DF
#     #parameterDf <- data.frame(dataScalingParameters.timeRange = 24,
#     #                          slope1Param_Estimate = 1,
#     #                          midPoint1Param_Estimate = 7,
#     #                          slope2Param_Estimate = 1,
#     #                          midPointDistanceParam_Estimate = 8,
#     #                          model = "doublesigmoidal")
#     #
#     # # Find the x value of the maximum for this model
#     # # with in the range of [0, dataScalingParameters.timeRange]
#     # midPoint1Param_x <- f_mid1_doublesigmoidal(parameterDf)
#     #

#     f_mid1_doublesigmoidal <- function(parameterDf){
#
#       max_x <- parameterDf$dataScalingParameters.timeRange
#       xmax <- stats::optimize(f1,
#                               interval = c(-1.125 * max_x, max_x * 3),
#                               tol = 0.0001,
#                               B1 = parameterDf$slope1Param_Estimate,
#                               M1 = parameterDf$midPoint1Param_Estimate,
#                               B2 = parameterDf$slope2Param_Estimate,
#                               L = parameterDf$midPointDistanceParam_Estimate,
#                               maximum = TRUE);
#       argumentt <- xmax$maximum;
#       constt <- f0(argumentt,
#                    B1 = parameterDf$slope1Param_Estimate,
#                    M1 = parameterDf$midPoint1Param_Estimate,
#                    B2 = parameterDf$slope2Param_Estimate,
#                    L = parameterDf$midPointDistanceParam_Estimate);
#
#       mid1x <- try(expr = stats::uniroot(f0mid,
#                                          interval = c(-1.125 * max_x, argumentt),
#                                          tol = 0.0001,
#                                          B1 = parameterDf$slope1Param_Estimate,
#                                          M1 = parameterDf$midPoint1Param_Estimate,
#                                          B2 = parameterDf$slope2Param_Estimate,
#                                          L = parameterDf$midPointDistanceParam_Estimate, const = constt), silent = TRUE);
#       if(inherits(mid1x, "try-error"))
#       {
#         rangeList <- seq(-1,30)
#         rangeListCounter <- -1
#         while(inherits(mid1x, "try-error") && rangeListCounter <= 30)
#         {
#           #print(rangeListCounter)
#           rangeListSelection <- rangeList[rangeListCounter + 2];
#           mid1x <- try(expr = stats::uniroot(f0mid,
#                                              interval = c((-1) * (2^rangeListSelection) * max_x, argumentt),
#                                              tol = 0.0001,
#                                              B1 = parameterDf$slope1Param_Estimate,
#                                              M1 = parameterDf$midPoint1Param_Estimate,
#                                              B2 = parameterDf$slope2Param_Estimate,
#                                              L = parameterDf$midPointDistanceParam_Estimate, const = constt), silent = TRUE);
#           rangeListCounter <- rangeListCounter + 1
#         }
#       }
#
#       mid1x <- mid1x$root
#
#       return(mid1x)
#     }
#**************************************


#**************************************
# f_mid2_doublesigmoidal is unchanged and is available in doublesigmoidalFitFunctions.R

#     # @title f_mid2_doublesigmoidal
#     #
#     # @param parameterDf is a dataFrame that includes columns with names:
#     #
#     # *"dataScalingParameters.timeRange" that represents maximum time of intensity measurements,
#     #
#     # parameters of double sigmoidal model:
#     #  *"slope1Param_Estimate",
#     #  *"midPoint1Param_Estimate",
#     #  *"slope2Param_Estimate",
#     #  *"midPointDistanceParam_Estimate".
#     #  (Note: additional parameters in Df does not cause problems)
#     #
#     # finally a column with the used model information model="doublesigmoidal".
#     #
#     # @description This function is designed to compensate for the small discrepancy in numerical parameters in double sigmoidal model. This function is called by numericalReCalculation.
#     # @return returns the x values of midpoint2, assuming the y value of midpoint2 is half of the distance between the maximum value it reaches and the final asymptote.
#     # @export
#     #
#     # @examples
#     # # Related example
#     #
#     ## Initial Command to Reset the System
#     #rm(list = ls())
#     #if (is.integer(dev.list())){dev.off()}
#     #cat("\014")
#     #
#     ## Generate a parameters DF
#     #parameterDf <- data.frame(dataScalingParameters.timeRange = 24,
#     #                          slope1Param_Estimate = 1,
#     #                          midPoint1Param_Estimate = 7,
#     #                          slope2Param_Estimate = 1,
#     #                          midPointDistanceParam_Estimate = 8,
#     #                          model = "doublesigmoidal")
#     #
#     # # Find the x value of the maximum for this model
#     # # with in the range of [0, dataScalingParameters.timeRange]
#     # midPoint2_x <- f_mid2_doublesigmoidal(parameterDf)
#     #
#     f_mid2_doublesigmoidal <- function(parameterDf){ # also untouched for now, but also a place we added try-error stop...
#
#       max_x <- parameterDf$dataScalingParameters.timeRange
#       xmax <- stats::optimize(f1,
#                               interval = c(-1.125 * max_x, max_x * 3),
#                               tol = 0.0001,
#                               B1 = parameterDf$slope1Param_Estimate,
#                               M1 = parameterDf$midPoint1Param_Estimate,
#                               B2 = parameterDf$slope2Param_Estimate,
#                               L = parameterDf$midPointDistanceParam_Estimate,
#                               maximum = TRUE);
#
#       argumentt <- xmax$maximum;
#       constt <- f0(argumentt,
#                    B1 = parameterDf$slope1Param_Estimate,
#                    M1 = parameterDf$midPoint1Param_Estimate,
#                    B2 = parameterDf$slope2Param_Estimate ,
#                    L = parameterDf$midPointDistanceParam_Estimate);
#
#       mid2x <- try(stats::uniroot(f0mid,
#                                   interval = c(argumentt, max_x * (1.25)),
#                                   tol = 0.0001,
#                                   B1 = parameterDf$slope1Param_Estimate,
#                                   M1 = parameterDf$midPoint1Param_Estimate,
#                                   B2 = parameterDf$slope2Param_Estimate,
#                                   L = parameterDf$midPointDistanceParam_Estimate, const=constt), silent = TRUE);
#
#       if(inherits(mid2x, "try-error")){
#         rangeList <- seq(-1,30)
#         rangeListCounter <- -1
#         while(inherits(mid2x, "try-error") && rangeListCounter <= 30){
#           rangeListSelection <- rangeList[rangeListCounter + 2];
#           mid2x <- try(stats::uniroot(f0mid,
#                                       interval = c(argumentt, max_x * (2^rangeListSelection)),
#                                       tol = 0.0001, B1 = parameterDf$slope1Param_Estimate,
#                                       M1 = parameterDf$midPoint1Param_Estimate,
#                                       B2 = parameterDf$slope2Param_Estimate,
#                                       L = parameterDf$midPointDistanceParam_Estimate, const = constt), silent = TRUE);
#           rangeListCounter <- rangeListCounter + 1
#         }
#       }
#
#       mid2x <- mid2x$root
#
#       return(mid2x)
#     }

#**************************************


#**************************************
# @title f_slope_doublesigmoidal_h0
# @param x as a single time point or a sequence of time points that one wants to calculate the derivative of.
# @param parameterDf is a dataFrame that includes columns with names:
#
# *"dataScalingParameters.timeRange" which represents the maximum time of intensity measurements,
#
# parameters of double sigmoidal model:
#  *"finalAsymptoteIntensityRatio_Estimate",
#  *"maximum_Estimate",
#  *"slope1Param_Estimate",
#  *"midPoint1Param_Estimate",
#  *"slope2Param_Estimate",
#  *"midPointDistanceParam_Estimate".
#
# finally a column with the used model information model="doublesigmoidal".
# @param timeStep is the time step for the derivative calculation increment "h". Default value is 0.00001
#
#
# @description The function calculates the numerical slope of double sigmoidal function with given parameters by using 5 points derivative. It is designed to compensate for the small discrepancy in numerical parameters in the double sigmoidal model. The function is called by numericalReCalculation.
# @return return the numerical slope of double sigmoidal for the given x value.
# @export
#
# @examples
# # Related examples
#
## Example 1
## Initial Command to Reset the System
#rm(list = ls())
#if (is.integer(dev.list())){dev.off()}
#cat("\014")
#
## Generate a parameters DF
#parameterDataframe <- data.frame(dataScalingParameters.timeRange = 24,
#                                 finalAsymptoteIntensityRatio_Estimate = .3,
#                                 maximum_Estimate = 4,
#                                 slope1Param_Estimate = 1,
#                                 midPoint1Param_Estimate = 7,
#                                 slope2Param_Estimate = 1,
#                                 midPointDistanceParam_Estimate = 8,
#                                 model = "doublesigmoidal")
#
# # Find the x value of the maximum for this model
# # with in the range of [0, dataScalingParameters.timeRange]
# slopeParam_Estimate <- f_slope_doublesigmoidal_h0(x = 6.5,
#                                                parameterDf = parameterDataframe,
#                                                timeStep = 0.00001)
#
## Example 2
## x might be vector as well
## Initial Command to Reset the System
#rm(list = ls())
#if (is.integer(dev.list())){dev.off()}
#cat("\014")
#
## Generate a parameters DF
#parameterDataframe <- data.frame(dataScalingParameters.timeRange = 24,
#                                 finalAsymptoteIntensityRatio_Estimate = .3,
#                                 maximum_Estimate = 4,
#                                 slope1Param_Estimate = 1,
#                                 midPoint1Param_Estimate = 7,
#                                 slope2Param_Estimate = 1,
#                                 midPointDistanceParam_Estimate = 8,
#                                 model = "doublesigmoidal")
#
# # Find the x value of the maximum for this model
# # with in the range of [0, dataScalingParameters.timeRange]
# slopeParam_Estimate <- f_slope_doublesigmoidal_h0(x = seq(-6.5, 50, 0.5),
#                                                parameterDf = parameterDataframe,
#                                                timeStep = 0.00001)
## Note that some points in the x sequence are out of the range of
## time <- [0, dataScalingParameters.timeRange]
#
f_slope_doublesigmoidal_h0 <- function(x, parameterDf, timeStep = 0.00001){

  fxp2h <- doublesigmoidalFitFormula_h0(x = x + 2 * timeStep,
                                     finalAsymptoteIntensityRatio = parameterDf$finalAsymptoteIntensityRatio_Estimate,
                                     maximum = parameterDf$maximum_Estimate,
                                     slope1Param = parameterDf$slope1Param_Estimate,
                                     midPoint1Param = parameterDf$midPoint1Param_Estimate,
                                     slope2Param = parameterDf$slope2Param_Estimate,
                                     midPointDistanceParam = parameterDf$midPointDistanceParam_Estimate,
                                     h0 = parameterDf$h0_Estimate)
  fxph <- doublesigmoidalFitFormula_h0(x = x + timeStep,
                                    finalAsymptoteIntensityRatio = parameterDf$finalAsymptoteIntensityRatio_Estimate,
                                    maximum = parameterDf$maximum_Estimate,
                                    slope1Param = parameterDf$slope1Param_Estimate,
                                    midPoint1Param = parameterDf$midPoint1Param_Estimate,
                                    slope2Param = parameterDf$slope2Param_Estimate,
                                    midPointDistanceParam = parameterDf$midPointDistanceParam_Estimate,
                                    h0 = parameterDf$h0_Estimate)

  fxmh <- doublesigmoidalFitFormula_h0(x = x - timeStep,
                                    finalAsymptoteIntensityRatio = parameterDf$finalAsymptoteIntensityRatio_Estimate,
                                    maximum = parameterDf$maximum_Estimate,
                                    slope1Param = parameterDf$slope1Param_Estimate,
                                    midPoint1Param = parameterDf$midPoint1Param_Estimate,
                                    slope2Param = parameterDf$slope2Param_Estimate,
                                    midPointDistanceParam = parameterDf$midPointDistanceParam_Estimate,
                                    h0 = parameterDf$h0_Estimate)

  fxm2h <- doublesigmoidalFitFormula_h0(x = x - 2 * timeStep,
                                     finalAsymptoteIntensityRatio = parameterDf$finalAsymptoteIntensityRatio_Estimate,
                                     maximum = parameterDf$maximum_Estimate,
                                     slope1Param = parameterDf$slope1Param_Estimate,
                                     midPoint1Param = parameterDf$midPoint1Param_Estimate,
                                     slope2Param = parameterDf$slope2Param_Estimate,
                                     midPointDistanceParam = parameterDf$midPointDistanceParam_Estimate,
                                     h0 = parameterDf$h0_Estimate)

  der <- (-1 * fxp2h + 8 * fxph - 8 * fxmh + 1 * fxm2h) / (12 * timeStep)
  return(der)
}
#**************************************


#**************************************
# @title sigmoidalRenormalizeParameters_h0 (This is an internal function)
# @param parametersDf is the parameter data frame generated by sigmoidal fit function and
#        includes the parameters named
#        *maximum_N_Estimate (normalized Maximum Estimate)
#        *finalAsymptoteIntensityRatio_N_Estimate (normalized final asymptote intensity estimate)
#        *slope1Param_N_Estimate (normalized Slope1Param Estimate)
#        *midPoint1Param_N_Estimate (normalized Midpoint1Param Estimate)
#        *slope2Param_N_Estimate (normalized Slope2Param Estimate)
#        *midPointDistanceParam_N_Estimate (normalized distance between midpoints estimate)
#.       *h0_N_Estimate (normalized estimate of lower asymptote)
#        *dataScalingParameters.intensityRange the range of initial unnormalized intensity. Provided if the data is normalized
#        *parameterDF$dataScalingParameters.intensityMin the minimum of unnormalized intensity. Provided if the data is normalized
#        *parameterDF$dataScalingParameters.timeRange the maximum time that the experiment reach. Provided if the data is normalized
# @param isalist checks if the input provided is a list (i.e normalized) or a data frame (i.e not normalized)
# @details If the fit was done on normalized data, parameter estimates will be on normalized scale.
#          This function unnormalizes those parameter estimations.
doublesigmoidalRenormalizeParameters_h0 <- function(parameterDF, isalist)
{
  if(isalist){
    parameterDF$finalAsymptoteIntensityRatio_Estimate <- (parameterDF$finalAsymptoteIntensityRatio_N_Estimate * parameterDF$maximum_N_Estimate * parameterDF$dataScalingParameters.intensityRange + parameterDF$dataScalingParameters.intensityMin) /
      (parameterDF$maximum_N_Estimate  * parameterDF$dataScalingParameters.intensityRange + parameterDF$dataScalingParameters.intensityMin)
    parameterDF$maximum_Estimate <- parameterDF$maximum_N_Estimate * parameterDF$dataScalingParameters.intensityRange + parameterDF$dataScalingParameters.intensityMin
    parameterDF$slope1Param_Estimate <- parameterDF$slope1Param_N_Estimate / parameterDF$dataScalingParameters.timeRange
    parameterDF$midPoint1Param_Estimate <- parameterDF$midPoint1Param_N_Estimate * parameterDF$dataScalingParameters.timeRange
    parameterDF$slope2Param_Estimate <- parameterDF$slope2Param_N_Estimate / parameterDF$dataScalingParameters.timeRange
    parameterDF$midPointDistanceParam_Estimate <- parameterDF$midPointDistanceParam_N_Estimate * parameterDF$dataScalingParameters.timeRange
    parameterDF$h0_Estimate <- parameterDF$h0_N_Estimate * parameterDF$dataScalingParameters.intensityRange + parameterDF$dataScalingParameters.intensityMin
  }
  if(!isalist){
    parameterDF$finalAsymptoteIntensityRatio_Estimate <- parameterDF$finalAsymptoteIntensityRatio_N_Estimate
    parameterDF$maximum_Estimate <- parameterDF$maximum_N_Estimate
    parameterDF$slope1Param_Estimate <- parameterDF$slope1Param_N_Estimate
    parameterDF$midPoint1Param_Estimate <- parameterDF$midPoint1Param_N_Estimate
    parameterDF$slope2Param_Estimate <- parameterDF$slope2Param_N_Estimate
    parameterDF$midPointDistanceParam_Estimate <- parameterDF$midPointDistanceParam_N_Estimate
    parameterDF$h0_Estimate <- parameterDF$h0_N_Estimate
  }
  return(parameterDF)
}
#**************************************

Try the sicegar package in your browser

Any scripts or data that you put into this service are public.

sicegar documentation built on Nov. 16, 2025, 1:07 a.m.