R/curveFit_utils.R

In peakPantheR: Peak Picking and Annotation of High Resolution Experiments

#' Check if object is of class peakPantheR_curveFit
#' 
#' Check if object is of class peakPantheR_curveFit
#' @param x object to test
#' @return (bool) TRUE or FALSE
is.peakPantheR_curveFit <- function(x){inherits(x, "peakPantheR_curveFit")}


#' Curve fitting using minpack.lm
#' 
#' Fit different curve models using minpack. Fitting parameters can be passed or guessed.
#'
#' @param x (numeric) x values (e.g. retention time)
#' @param y (numeric) y observed values (e.g. spectra intensity)
#' @param curveModel (str) name of the curve model to fit (currently \code{skewedGaussian})
#' @param params (list or str) either 'guess' for automated parametrisation or list of initial parameters ($init_params), lower parameter bounds ($lower_bounds) and upper parameter bounds ($upper_bounds)
#'
#' @return A 'peakPantheR_curveFit': a list of fitted curve parameters, \code{fitStatus} from \code{nls.lm$info} and curve shape name \code{curveModel}. \code{fitStatus=0} unsuccessful completion: improper input parameters, \code{fitStatus=1} successful completion: first convergence test is successful, \code{fitStatus=2} successful completion: second convergence test is successful, \code{fitStatus=3} successful completion: both convergence test are successful, \code{fitStatus=4} questionable completion: third convergence test is successful but should be carefully examined (maximizers and saddle points might satisfy), \code{fitStatus=5} unsuccessful completion: excessive number of function evaluations/iterations
#' 
#' @examples
#' \dontrun{
#' ## x is retention time, y corresponding intensity
#' input_x  <- c(3362.102, 3363.667, 3365.232, 3366.797, 3368.362, 3369.927, 3371.492, 3373.057,
#'              3374.622, 3376.187, 3377.752, 3379.317, 3380.882, 3382.447, 3384.012, 3385.577,
#'              3387.142, 3388.707, 3390.272, 3391.837, 3393.402, 3394.966, 3396.531, 3398.096,
#'              3399.661, 3401.226, 3402.791, 3404.356, 3405.921, 3407.486, 3409.051)
#' input_y  <- c(51048, 81568, 138288, 233920, 376448, 557288, 753216, 938048, 1091840, 1196992,
#'               1261056, 1308992, 1362752, 1406592, 1431360, 1432896, 1407808, 1345344, 1268480,
#'               1198592, 1126848, 1036544, 937600, 849792, 771456, 692416, 614528, 546088, 492752,
#'               446464, 400632)
#' 
#' ## Fit
#' fitted_curve <- fitCurve(input_x, input_y, curveModel='skewedGaussian', params='guess')
#' 
#' ## Returns the optimal fitting parameters
#' fitted_curve
#' #
#' # $amplitude
#' # [1] 275371.1
#' # 
#' # $center
#' # [1] 3382.577
#' # 
#' # $sigma
#' # [1] 0.07904697
#' # 
#' # $gamma
#' # [1] 0.001147647
#' # 
#' # $fitStatus
#' # [1] 2
#' # 
#' # $curveModel
#' # [1] "skewedGaussian"
#' #
#' # attr(,"class")
#' # [1] "peakPantheR_curveFit"
#' }
fitCurve <- function(x, y, curveModel='skewedGaussian', params='guess') {
  
  ## Check inputs
  # x and y length
  if (length(x) != length(y)) {
    stop('Error: length of "x" and "y" must match!')
  }
  # known curveModel
  known_curveModel <- c('skewedGaussian')
  if (!(curveModel %in% known_curveModel)) {
    stop(paste('Error: "curveModel" must be one of:', known_curveModel))
  }
  # params
  if (!(typeof(params) %in% c('list', 'character'))) {
    stop('Error: "params" must be a list or "guess"')
  }
  useGuess = TRUE
  if (any(params != "guess")) {
    useGuess = FALSE
    # check init_params, lower and upper bounds are defined
    if (!all(c("init_params", "lower_bounds", "upper_bounds") %in% names(params))) {
      stop('Error: "params must be a list of "init_params", "lower_bounds" and "upper_bounds"')
    }
    # init_params is list
    if (typeof(params$init_params) != 'list') {
      stop('Error: "params$init_params" must be a list of parameters')
    }
    # lower_bounds is list
    if (typeof(params$lower_bounds) != 'list') {
      stop('Error: "params$lower_bounds" must be a list of parameters')
    }
    # upper_bounds is list
    if (typeof(params$upper_bounds) != 'list') {
      stop('Error: "params$upper_bounds" must be a list of parameters')
    }
  }
  
  ## Init
  fittedCurve <- list()
  
  ## Run fitting
  # skewed gaussian
  if (curveModel == 'skewedGaussian') {
    
    # Guess parameters and bounds
    if (useGuess) {
      new_params   <- skewedGaussian_guess(x, y)
    } else {
      new_params   <- params
    }
    
    # ensure order of init params and bounds (init is a list, lower and upper are ordered numeric vectors)
    init  <- list(amplitude = new_params$init_params$amplitude,
                  center = new_params$init_params$center,
                  sigma = new_params$init_params$sigma,
                  gamma = new_params$init_params$gamma)
    lower <- unlist( c(new_params$lower_bounds['amplitude'], new_params$lower_bounds['center'], new_params$lower_bounds['sigma'], new_params$lower_bounds['gamma']) )
    upper <- unlist( c(new_params$upper_bounds['amplitude'], new_params$upper_bounds['center'], new_params$upper_bounds['sigma'], new_params$upper_bounds['gamma']) )
    
    # perform fit
    resultFit  <- minpack.lm::nls.lm(par=init, lower=lower, upper=upper, fn=skewedGaussian_minpack.lm_objectiveFun, observed=y, xx=x)
    
    # prepare output
    fittedCurve             <- resultFit$par
    fittedCurve$fitStatus   <- resultFit$info
    fittedCurve$curveModel  <- curveModel
    class(fittedCurve)      <- 'peakPantheR_curveFit'
  }
  # for future curve shapes
  #} else if () {
  #}
  
  return(fittedCurve)
}


#' Predict curve values
#'
#' Evaluate fitted curve values at \code{x} data points
#'
#' @param fittedCurve (peakPantheR_curveFit) A 'peakPantheR_curveFit': a list of curve fitting parameters, curve shape model \code{curveModel} and nls.lm fit status \code{fitStatus}.
#' @param x (numeric) values at which to evaluate the fitted curve
#' 
#' @return  fitted curve values at x
#' 
#' @examples
#' \dontrun{
#' ## Input a fitted curve
#' fittedCurve <- list(amplitude=275371.1, center=3382.577, sigma=0.07904697, gamma=0.001147647,
#'                     fitStatus=2, curveModel="skewedGaussian")
#' class(fittedCurve)  <- 'peakPantheR_curveFit'
#' input_x     <- c(3290, 3300, 3310, 3320, 3330, 3340, 3350, 3360, 3370, 3380, 3390, 3400, 3410)
#'
#' ## Predict y at each input_x
#' pred_y      <- predictCurve(fittedCurve, input_x)
#' pred_y
#' # [1] 2.347729e-08 1.282668e-05 3.475590e-03 4.676579e-01 3.129420e+01 1.043341e+03 1.736915e+04
#' # [8] 1.447754e+05 6.061808e+05 1.280037e+06 1.369651e+06 7.467333e+05 2.087477e+05
#' }
predictCurve  <- function(fittedCurve, x) {

  # Check input
  if (!is.peakPantheR_curveFit(fittedCurve)) {
    stop('Error: "fittedCurve" must be a peakPantheR_curveFit!')
  }
  known_curveModel <- c('skewedGaussian')
  if (!(fittedCurve$curveModel %in% known_curveModel)) {
    stop(paste('Error: "fittedCurve$curveModel" must be one of:', known_curveModel))
  }
  
  # Select correct model
  if (fittedCurve$curveModel == 'skewedGaussian') {
    yy <- skewedGaussian_minpack.lm(params=fittedCurve, xx=x)
  } # for future curve shapes
  #} else if () {
  #}
  
  return(yy)
}



## --------------------------------------------------------------------------------------------------
##        Skewed Gaussian
## --------------------------------------------------------------------------------------------------

#' Gaussian Error function
#'
#' Implementation of the gaussian error function
#' 
#' @param x (numeric) value at which to evaluate the gaussian error function
#' 
#' @return Value of the gaussian error function evaluated at x
skew_erf   <- function(x){
  return(2 * stats::pnorm(x * sqrt(2)) - 1)
} 


#' Implementation of the Skewed Gaussian peak shape for use with minpack.lm
#'
#' Implementation of the Skewed Gaussian peak shape for use with minpack.lm
#'
#' @param params (list) skewed gaussian parameters (\code{params$gamma}, \code{params$center}, \code{params$sigma}, \code{params$amplitude})
#' @param xx (numeric) values at which to evalute the skewed gaussian
#' 
#' @return value of the skewed gaussian evaluated at xx
skewedGaussian_minpack.lm   <- function(params, xx) {
  erf_term  <- 1 + skew_erf((params$gamma * (xx - params$center)) / params$sigma * sqrt(2))
  yy        <- (params$amplitude / (params$sigma * sqrt(2* pi))) * exp(-(xx - params$center)^2 / 2*params$sigma^2) * erf_term
  
  return(yy)
}


#' Skewed Gaussian minpack.lm objective function
#'
#' Skewed Gaussian minpack.lm objective function, calculates residuals using the skewed gaussian Peak Shape
#'
#' @param params (list) skewed gaussian parameters (\code{params$gamma}, \code{params$center}, \code{params$sigma}, \code{params$amplitude})
#' @param observed (numeric) observed y value at xx
#' @param xx (numeric) value at which to evalute the skewed gaussian
#' 
#' @return difference between observed and expected skewed gaussian value evaluated at xx
skewedGaussian_minpack.lm_objectiveFun  <- function(params, observed, xx) {
  return(observed - skewedGaussian_minpack.lm(params, xx))
}


#' Guess function for initial skewed gaussian parameters and bounds
#'
#' Guess function for initial skewed gaussian parameters and bounds, at the moment only checks the x position
#'
#' @param x (numeric) x values (e.g. retention time)
#' @param y (numeric) y observed values (e.g. spectra intensity)
#' 
#' @return A list of guessed starting parameters \code{list()$init_params}, lower \code{list()$lower_bounds} and upper bounds \code{list()$upper_bounds} (\code{$gamma}, \code{$center}, \code{$sigma}, \code{$amplitude})
skewedGaussian_guess <- function(x, y) {
  # set center as x position of max y value (e.g. highest spectra intensity)
  center_guess  <- x[which.max(y)]
  # init_param
  init_params   <- list(amplitude = 10000000, center = center_guess,     sigma = 1, gamma = 1)
  # lower_bounds
  lower_bounds  <- list(amplitude = 0,        center = center_guess - 3, sigma = 0, gamma = -0.1)
  # upper_bounds
  upper_bounds  <- list(amplitude = 1e9,      center = center_guess + 3, sigma = 5, gamma = 5)
  
  return(list(init_params=init_params, lower_bounds=lower_bounds, upper_bounds=upper_bounds))
}