cgaim: Constrained Groupwise Additive Index Models

Documented in cgaim.control

################################################################################
#
#  Setting up parameters to control for the estimation
#
################################################################################

#' Parameters controlling the CGAIM fit
#' 
#' Internal function setting the parameters to control the CGAIM fit. Sets default values and check parameters passed by the user. It is called internally by \code{\link{cgaim}} and should not be called directly by the user.
#' 
#' @param max.iter The maximum number of iteration allowed. 
#' @param tol The tolerance for convergence. The algorithm stops when the convergence criterion falls below \code{tol}.
#' @param halving Logical turning on and off halving for bad steps. See details.
#' @param min.step.len Numeric between 0 and 1 giving the minimum descent step length to be considered in case of bad step. See details.
#' @param convergence_criterion Character indicating the convergence criterion for the algorithm. See details.
#' @param trace If TRUE, prints the convergence criterion and alpha coefficients at each step.
#' @param alpha.start An optional vector or list of starting alpha values. If \code{NULL}, starting values are generated internally. See the \code{init.type} argument.
#' @param init.type The type of initialization to perform if no initial value is provided. If \code{init.type = "regression"} (the default), starting values are generated by regressing the index design matrix on the response. If \code{init.type = "random"}, feasible starting values are generated randomly.
#' @param norm.type The type of norm used to normalize index coefficients vectors. See \code{\link[base]{norm}} for available norms. Default to L1 norm meaning that, for each index, absolute values of coefficients sum to 1.
#' @param check.Cmat Logical indicating whether to check for redundant constraints and remove them from \code{Cmat}.
#' @param solver The quadratic programming solver to use. One of \code{"osqp"} (the default), \code{"quadprog"} or \code{"coneproj"}.
#' @param ctol Tolerance value on constraints. See details.
#' @param qp_pars A named list of parameters to be passed to the \code{solver} function. 
#' @param sample_pars A named list of parameters to be passed to \code{\link[limSolve]{xsample}} when randomly generating initial alpha coefficients.
#' @param sm_method Character specifying which method to use for constrained smoothing. Either \code{\link[scam]{scam}} (the default), \code{\link[cgam]{cgam}} or \code{\link[scar]{scar}}.
#' @param sm_pars Named list to pass specific parameters to the smoothing function of \code{sm_method}. See help pages of corresponding functions.
#' 
#' @details The model is fitted by an iterative sequential quadratic programming algorithm. The algorithm iterates between updating the index alpha coefficients and updating the smoothing of indices and covariates. It stops when the criterion given in \code{convergence_criterion} is below \code{tol}, or when \code{max.iter} is reached. Convergence criteria include \code{rss} for which the algorithm stops when the relative decrease in residual sum of squares, \code{(rss_new - rss_old) / rss_old} is below \code{tol}, \code{alpha} for which the algorithm stops when the largest update \code{max(abs(alpha_new - alpha_old) / abs(alpha_old))} is below \code{tol}, and \code{offset} when the scalar product between the RSS and current direction (measuring orthogonality) is below \code{tol} (EXPERIMENTAL, use at your own risk).
#' 
#' By default, when the RSS fails to decrease during a step (a "bad" step), the step length is iteratively halved until the RSS decreases. The minimum step length allowed is controlled by \code{min.step.len} as the proportion of the original step length. This is a common behaviour in non-linear least squares and is implement in \code{\link[stats]{nls}} for instance, but can be turned off by setting \code{halving = FALSE}, in which case the algorithm stops for any bad step. 
#' 
#' The alpha updating step consists in estimating an update vector in a descent direction by a constrained regression of index derivatives on the current residuals of the model. This is fitted through a quadratic program, ensuring the updated coefficients respect the constraints at each step of the algorithm. Initial values can either be provided by the user through the argument \code{alpha.start} or be internally generated. The latter is controlled by the argument \code{init.type} allowing to initialize the weights either by regressing the index variables on the response (\code{init.type = "regression"}) ensuring feasible starting values (the default), or by randomly generating feasible values (\code{init.type = "random"}). In the latter case, random generation is performed by the function \code{\link[limSolve]{xsample}} which can be controlled by the parameter \code{sample_pars}. When random initial values are chosen, it is recommended to fit the algorithm several time and keep the best fit, to avoid falling into a local minimum. 
#' 
#' At the moment, three solvers are available to perform quadratic programming, which can be controlled by the argument \code{solver}. By default the function \code{\link[osqp]{solve_osqp}} (\code{solver = "osqp"}) is used. Alternatively the more established but slower function \code{\link[quadprog]{solve.QP}} (\code{solver = "quadprog"}) as well as \code{\link[coneproj]{qprog}} (\code{solver = "coneproj"}) functions can be used. Although default parameters are internally set for these function, they can entirely be controlled through the argument \code{qp_pars}. See their specific help pages for details.
#' 
#' In some cases, minimal numerical imprecision in the repeated call to quadratic program, along with the normalization of alpha coefficients, can lead to unfeasible alphas at convergence. To avoid this, these imprecision are compensated by adding a small tolerance \code{ctol} to the constraints, defaulting to 0.001. If no tolerance is wanted, it can be set to 0.
#' 
#' By default, the package automatically checks that \code{Cmat} is irreducible, i.e. that no constraint is redundant. A constraint is redundant if it can be expressed as a non-negative linear combination of other constraints. If \code{check.Cmat = TRUE}, such constraints are removed with a warning.
#' 
#' @returns A named list containing all arguments to be used in \code{\link{cgaim}}.
#' 
#' @seealso These parameters control the fitting of \code{\link{cgaim}}.
#' 
#' @references 
#' Bates, D.M., Watts, D.G., 1981. A Relative Off set Orthogonality Convergence Criterion for Nonlinear least Squares. Technometrics 23, 179–183.
#' 
#' Bates, D.M., Watts, D.G., 1988. Nonlinear Regression Analysis and Its Applications, Wiley Series in Probability and Statistics. Wiley.
#' 
#' @export
#' 
cgaim.control <- function(
  # For the overall algorithm  
  max.iter = 50, tol = 1e-3, halving = TRUE, min.step.len = 0.1, 
  convergence_criterion = "rss", trace = FALSE,
  # For alpha estimation
  alpha.start = NULL, init.type = "regression", norm.type = "1", check.Cmat = TRUE, 
  solver = "osqp", ctol = 1e-3, qp_pars = list(), sample_pars = list(),
  # For smoothing
  sm_method = "scam", sm_pars = list())
{
  # Catch passed and default parameters in list
  pars <- as.list(match.call())[-1]
  defpars <- as.list(formals(cgaim.control))
  # Merge them
  pars <- utils::modifyList(defpars, pars)
  # Check specific parameters
  pars$convergence_criterion <- match.arg(pars$convergence_criterion,
    c("rss", "alpha", "offset"))
  pars$init.type <- match.arg(pars$init.type, c("regression", "random"))
  pars$solver <- match.arg(pars$solver, c("quadprog", "osqp"))
  pars$sm_method <- match.arg(pars$sm_method, c("scam", "cgam", "scar"))
  # Default parameters for the QP solver
  pars$qp_pars <- do.call(sprintf("def.%s", pars$solver), list(pars$qp_pars))
  # Evaluate lists
  lang <- sapply(pars, is.language)
  pars[lang] <- lapply(pars[lang], eval)
  # Return
  pars
}