Commit 4ef12688 authored by davidkep's avatar davidkep

Merge branch 'release/2.0.1'

parents 99d7844f 0f7be3b1
......@@ -5,4 +5,5 @@
*.Rproj
.DS_Store
.Rbuildignore
.gitignore
*.trace
Package: pense
Type: Package
Title: Penalized Elastic Net S/MM-Estimator of Regression
Version: 2.0.1
Date: 2019-10-24
Authors@R: c(
person("David", "Kepplinger", , "david.kepplinger@gmail.com",
role = c("aut", "cre")),
person("Matias", "Salibian-Barrera", role = c("aut")),
person("Gabriela", "Cohen Freue", role = "aut"),
person("Derek", "Cho", role = "ctb")
)
Copyright: See the file COPYRIGHTS for copyright details on some of the
functions and algorithms used.
Encoding: UTF-8
Biarch: true
SystemRequirements: C++11
URL: https://gitlab.math.ubc.ca/dakep/pense
BugReports: https://gitlab.math.ubc.ca/dakep/pense/issues
Description: Robust penalized elastic net S and MM estimator for linear
regression. The method is described in detail in
Cohen Freue, G. V., Kepplinger, D., Salibian-Barrera, M., and Smucler, E.
(2017) <https://gcohenfr.github.io/pdfs/PENSE_manuscript.pdf>.
Depends:
R (>= 3.4.0),
Matrix
Imports:
Rcpp,
parallel,
methods,
lifecycle
LinkingTo:
nsoptim,
Rcpp,
RcppArmadillo (>= 0.9.100)
Suggests:
testthat (>= 2.1.0)
License: MIT + file LICENSE
NeedsCompilation: yes
RoxygenNote: 6.1.1
Roxygen: list(markdown = TRUE)
MIT License
Copyright (c) 2019 David Kepplinger
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
YEAR: 2019
COPYRIGHT HOLDER: David Kepplinger
# Generated by roxygen2: do not edit by hand
S3method(print,nsoptim_metrics)
export(consistency_const)
export(elnet)
export(en_admm_options)
export(en_dal_options)
export(enpy_initial_estimates)
export(enpy_options)
export(mloc)
export(mlocscale)
export(mscale)
export(mscale_algorithm_options)
export(pense)
export(pense_admm_options)
export(pense_mm_options)
export(rho_function)
export(tau_size)
importFrom(Rcpp,evalCpp)
importFrom(lifecycle,deprecate_soft)
importFrom(methods,as)
importFrom(stats,mad)
useDynLib(pense, .registration = TRUE)
#' Options for the ENPY Algorithm
#'
#' @param max_it maximum number of PY iterations.
#' @param eps numerical tolerance to check for convergence.
#' @param en_algorithm_opts options for the LS-EN algorithm. See [en_algorithm_options] for details.
#' @param keep_psc_proportion how many observations should be kept based on the Principal Sensitivy Components.
#' @param keep_residuals_measure how to determine what observations to keep, based on their residuals.
#' If `proportion`, a fixed number of observations is kept, while if `threshold`,
#' only observations with residuals below the threshold are kept.
#' @param keep_residuals_proportion how many observations should be kept based on their residuals.
#' @param keep_residuals_threshold only observations with (standardized) residuals less than this threshold are kept.
#' @param retain_best_factor in addition to the candidates from the last iteration, also keep candidates
#' that are within this factor of the best candidate.
#'
#' @return options for the ENPY algorithm.
#' @export
enpy_options <- function (max_it = 10, eps = 1e-6, keep_psc_proportion = 0.5,
en_algorithm_opts,
keep_residuals_measure = c('threshold', 'proportion'),
keep_residuals_proportion = 0.5, keep_residuals_threshold = 2,
retain_best_factor = 1.1) {
list(max_it = .as(max_it[[1L]], 'integer'),
en_options = if (missing(en_algorithm_opts)) { NULL } else { en_algorithm_opts },
eps = .as(eps[[1L]], 'numeric'),
keep_psc_proportion = .as(keep_psc_proportion[[1L]], 'numeric'),
use_residual_threshold = match.arg(keep_residuals_measure) == 'threshold',
keep_residuals_proportion = .as(keep_residuals_proportion[[1L]], 'numeric'),
keep_residuals_threshold = .as(keep_residuals_threshold[[1L]], 'numeric'),
retain_best_factor = .as(retain_best_factor[[1L]], 'numeric'))
}
#' Options for the M-scale Estimation Algorithm
#'
#' @param max_it maximum number of iterations.
#' @param eps numerical tolerance to check for convergence.
#'
#' @return options for the M-scale estimation algorithm.
#' @export
mscale_algorithm_options <- function (max_it = 200, eps = 1e-8) {
list(max_it = .as(max_it[[1L]], 'integer'),
eps = .as(eps[[1L]], 'numeric'))
}
#' Full options for the M-scale Estimation Algorithm
#'
#' @param mscale_opts "public" control options created by [mscale_algorithm_options].
#' @param bdp the breakdown point, i.e., `delta` in the M-estimation equation.
#' @param cc the cutoff threshold for the bisquare rho function.
#'
#' @return full options for the M-scale estimation algorithm.
.full_mscale_algo_options <- function (bdp, cc, mscale_opts) {
mscale_opts$delta <- .as(bdp, 'numeric')
if (isTRUE(mscale_opts$delta < .Machine$double.eps) || isTRUE(mscale_opts$delta > 0.5)) {
stop("`bdp` is outside of 0 and 0.5")
}
mscale_opts$cc <- if (is.null(cc)) { .bisquare_consistency_const(mscale_opts$delta) } else { .as(cc, 'numeric') }
return(mscale_opts)
}
#' Control the Algorithm to Compute Penalized Elastic Net S-Estimates
#'
#' The package provides different algorithms to compute the PENSE estimate.
#' To select a specific algorithm and set its parameters, use any of the `pense_***_options` functions.
#'
#' * [pense_admm_options]: Select the ADMM algorithm.
#' * [pense_mm_options]: Select the MM-algorithm.
#'
#' @name pense_algorithm_options
NULL
#' Options for the S-Estimate Algorithm
#'
#' @param max_it maximum number of iterations.
#' @param sparse use sparse coefficients.
#' @param tau step size for the algorithm.
#' @param prox_eps numerical tolerance for computing the proximal operator of the S-loss.
#' @param prox_max_it maximum number of iterations for computing the proximal operator of the S-loss.
#' @param prox_oscillate_window moving average size to determine oscillation for computing the proximal operator of the
#' S-loss.
#' @param prox_minimum_step_size minimum step size for computing the proximal operator of the S-loss.
#' @param prox_wolfe_c1 constant to check the first Wolfe condition for computing the proximal operator of the
#' S-loss.
#' @param prox_wolfe_c2 constant to check the second Wolfe condition for computing the proximal operator of the
#' S-loss.
#' @param step_size_adj multiplicative factor to decrease the step size if Wolfe's conditions are not satisfied
#' when computing the proximal operator of the S-loss.
#' @param prox_max_small_steps maximum number of consecutive small steps when computing the proximal operator of the
#' S-loss. After this many steps with minimal step size, a large step is performed to
#' escape the current neighborhood.
#'
#' @return options for the S-Estimate algorithm.
#' @export
pense_admm_options <- function (max_it = 5000, tau, sparse = TRUE, prox_eps = 1e-8, prox_max_it = 200,
prox_oscillate_window = 4, prox_minimum_step_size = 0.01, prox_wolfe_c1 = 1e-4,
prox_wolfe_c2 = 0.9, prox_step_size_adj = 0.9, prox_max_small_steps = 10) {
list(algorithm = 'admm',
max_it = .as(max_it[[1L]], 'integer'),
prox_opts = list(
eps = .as(prox_eps[[1L]], 'numeric'),
max_it = .as(max_it[[1L]], 'integer'),
tau = if (missing(tau) || is.null(tau)) { -1 } else { .as(tau[[1L]], 'numeric') },
oscillate_window = .as(prox_oscillate_window[[1L]], 'integer'),
minimum_step_size = .as(prox_minimum_step_size[[1L]], 'numeric'),
wolfe_ss_c1 = .as(prox_wolfe_c1[[1L]], 'numeric'),
wolfe_ss_c2 = .as(prox_wolfe_c2[[1L]], 'numeric'),
ss_adjustment = .as(prox_step_size_adj[[1L]], 'numeric'),
max_minimum_step_size_steps = .as(prox_max_small_steps[[1]], 'integer')
),
sparse = isTRUE(sparse[[1L]]))
}
#' MM-algorithm to Compute Penalized Elastic Net S-Estimates
#'
#' @param max_it maximum number of iterations.
#' @param tightening how to make inner iterations more precise as the algorithm approaches a local minimum.
#' @param tightening_steps for *adaptive* tightening strategy, how often to tighten until the desired tolerance
#' is attained.
#' @param en_algorithm_opts options for the inner LS-EN algorithm. See [en_algorithm_options] for details.
#'
#' @return options for the PENSE MM-algorithm.
#' @export
pense_mm_options <- function (max_it = 500, tightening = c('adaptive', 'exponential', 'none'),
tightening_steps = 10, en_algorithm_opts) {
list(algorithm = 'mm',
max_it = .as(max_it[[1L]], 'integer'),
tightening = .tightening_id(match.arg(tightening)),
tightening_steps = .as(tightening_steps[[1L]], 'integer'),
en_options = if (missing(en_algorithm_opts)) { NULL } else { en_algorithm_opts })
}
#' Control the Algorithm to Compute (Weighted) Least-Squares Elastic Net Estimates
#'
#' The package supports multiple different algorithms to compute the EN estimate for weighted LS loss functions.
#' Each algorithm has certain characteristics that make it useful for some problems.
#' To select a specific algorithm and set its parameters, use any of the `en_***_options` functions.
#'
#' * [en_admm_options]: Select an iterative ADMM-type algorithm. There are two versions available:
#' `admm_type = "linearized"` needs *O(n p)* operations per iteration and converges linearly, while
#' `admm_type = "var-stepsize"` needs *O(n p^3)* operations per iteration but converges quadratically.
#' * [en_dal_options]: Select the iterative Dual Augmented Lagrangian (DAL) method. DAL needs O(n^3 p^2) operations
#' per iteration, but converges exponentially.
#'
#' @name en_algorithm_options
NULL
#' Options for the ADMM Elastic Net Algorithm
#'
#' @param max_it maximum number of iterations.
#' @param eps numerical tolerance to check for convergence.
#' @param sparse use sparse coefficients.
#' @param admm_type what type of the ADMM algorithm to use. If `linearized`,
#' uses a linearized version of ADMM which has runtime $O()$
#' and converges linearly.
#' If `var-stepsize`, uses a variable step-size ADMM
#' algorithm which converges quadratically for "true" EN
#' penalties (i.e., \eqn{alpha < 1}) but has runtime $O()$.
#' If `auto` (the default), chooses the type based on the
#' penalty and the problem size.
#' @param tau step size for the algorithm if using the `linearized` version
#' and the largest step size if using the `var-stepsize` version.
#' @param tau_adjustment_lower (smallest) multiplicative factor for the
#' adjustment of the step size
#' `tau = tau_adjustment * tau`
#' (only for the `var-stepsize` version).
#' @param tau_adjustment_upper (largest) multiplicative factor for the
#' adjustment of the step size
#' `tau = tau_adjustment * tau`
#' (only for the `var-stepsize` version).
#'
#' @return options for the ADMM EN algorithm.
#' @family EN algorithms
#' @export
en_admm_options <- function (max_it = 1000, eps = 1e-8, tau, sparse = FALSE,
admm_type = c('auto', 'linearized', 'var-stepsize'),
tau_lower_mult = 0.01, tau_adjustment_lower = 0.98,
tau_adjustment_upper = 0.999) {
tau <- if (missing(tau) || is.null(tau)) { -1 } else { .as(tau[[1L]], 'numeric') }
list(algorithm = 'admm',
admm_type = match.arg(admm_type),
sparse = isTRUE(sparse[[1L]]),
max_it = .as(max_it[[1L]], 'integer'),
eps = .as(eps[[1L]], 'numeric'),
tau = tau,
prox_opts = list(tau = tau),
tau_lower_mult = .as(tau_lower_mult[[1L]], 'numeric'),
tau_adjustment_lower = .as(tau_adjustment_lower[[1L]], 'numeric'),
tau_adjustment_upper = .as(tau_adjustment_upper[[1L]], 'numeric'))
}
#' Options for the DAL Elastic Net Algorithm
#'
#' @param max_it maximum number of (outer) iterations.
#' @param max_inner_it maximum number of (inner) iterations in each outer iteration.
#' @param eps numerical tolerance to check for convergence.
#' @param eta_multiplier multiplier for the barrier parameter. In each iteration, the barrier must be more restrictive
#' (i.e., the multiplier must be > 1).
#' @param eta_start_conservative conservative initial barrier parameter. This is used if the previous penalty is
#' undefined or too far away.
#' @param eta_start_aggressive aggressive initial barrier parameter. This is used if the previous penalty is close.
#' @param lambda_relchange_aggressive how close must the lambda parameter from the previous penalty term be to use
#' an aggressive initial barrier parameter (i.e., what constitues "too far").
#'
#' @return options for the DAL EN algorithm.
#' @family EN algorithms
#' @export
en_dal_options <- function (max_it = 100, max_inner_it = 100, eps = 1e-6, eta_multiplier = 2,
eta_start_conservative = 0.01, eta_start_aggressive = 1,
lambda_relchange_aggressive = 0.25) {
list(algorithm = 'dal',
sparse = TRUE,
max_it = .as(max_it[[1L]], 'integer'),
max_inner_it = .as(max_inner_it[[1L]], 'integer'),
eps = .as(eps[[1L]], 'numeric'),
eta_start_numerator_conservative = .as(eta_start_conservative[[1L]], 'numeric'),
eta_start_numerator_aggressive = .as(eta_start_aggressive[[1L]], 'numeric'),
lambda_relchange_aggressive = .as(lambda_relchange_aggressive[[1L]], 'numeric'),
eta_multiplier = .as(eta_multiplier[[1L]], 'numeric'))
}
#' Ridge optimizer using an Augmented data matrix.
#' Only available for Ridge problems (`alpha=0``) and selected automatically in this case.
en_ridge_options <- function () {
list(algorithm = 'augridge', sparse = FALSE)
}
## Check if the selected EN algorithm can handle the given penalty.
.check_en_algorithm <- function (en_algorithm_opts, alpha) {
if (en_algorithm_opts$algorithm == 'dal') {
if(!isTRUE(alpha > 0)) {
warning('The DAL algorithm can not handle a Ridge penalty. Using default
algorithm as fallback.')
return(FALSE)
}
}
return(TRUE)
}
## Choose the appropriate ADMM algorithm type based on the penalty and
## the problem size.
.choose_admm_algorithm <- function (en_algorithm_opts, alpha, x) {
if (isTRUE(en_algorithm_opts$admm_type != 'auto')) {
return(en_algorithm_opts)
}
en_algorithm_opts$admm_type <- if (isTRUE(alpha < 1 && ncol(x) < 1.4 * nrow(x))) {
'var-stepsize'
} else {
'linearized'
}
return(en_algorithm_opts)
}
.en_algorithm_id <- function (en_algorithm_opts) {
switch (en_algorithm_opts$algorithm,
admm = switch(en_algorithm_opts$admm_type, `var-stepsize` = 2L, 1L),
dal = 3L,
augridge = 4L,
1L)
}
.pense_algorithm_id <- function (opts) {
switch (opts$algorithm,
admm = 2L,
mm = 1L,
1L)
}
.tightening_id <- function (tightening) {
switch (tightening, exponential = 1L, adaptive = 2L, 0L)
}
#' A wrapper around `methods::as` which raises an error if the conversion results in NA.
#' @param ... passed on to [methods::as].
#'
#' @importFrom methods as
.as <- function (object, class, ...) {
object_var <- deparse(substitute(object))
tryCatch(methods::as(object, class, ...), warning = function (w) {
stop(sprintf('`%s` is not of type `%s`', object_var, class))
})
}
.select_en_algorithm <- function (en_options, alpha, x) {
if (!is.null(en_options)) {
# Check if the EN algorithm can handle the given `alpha`
if (en_options$algorithm == 'dal') {
if(!isTRUE(alpha > 0)) {
warning('The DAL algorithm can not handle a Ridge penalty. Using default algorithm as fallback.')
en_options <- NULL
}
}
}
if (is.null(en_options)) {
en_options <- if (alpha > 0) {
en_admm_options()
} else {
en_ridge_options()
}
}
if (en_options$algorithm == 'admm') {
# Choose a "good" ADMM algorithm for the given problem size.
if (isTRUE(en_options$admm_type == 'auto')) {
en_options$admm_type <- if (isTRUE(alpha < 1 && ncol(x) < 1.4 * nrow(x))) {
'var-stepsize'
} else {
'linearized'
}
}
}
en_options$algorithm <- .en_algorithm_id(en_options)
return(en_options)
}
#' Compute the Elastic Net Regularization Path
#'
#' Compute the EN estimator for linear regression with optional observation
#' weights and penalty loadings.
#'
#' The elastic net estimator for the linear regression model solves
#' the optimization problem
#'
#' \deqn{argmin_{\mu, \beta}
#' (1/n) \sum_i w_i (y_i - \mu - x_i' \beta)^2 +
#' \lambda \sum_j 0.5 (1 - \alpha) \beta_j^2 + \alpha l_i |\beta_j| }
#'
#' with observation weights \eqn{w_i} and penalty loadings \eqn{l_i}.
#'
#' @param alpha value for the alpha parameter, the balance between L1 and L2
#' penalization.
#' @param lambdas a vector of positive values for the lambda parameter.
#' @param penalty_loadings a vector of positive penalty loadings
#' (a.k.a. weights) for different penalization of each
#' coefficient.
#' @param weights a vector of positive observation weights.
#' @param include_intercept include an intercept in the model.
#' @param en_algorithm_opts options for the EN algorithm. See [en_algorithm_options]
#' for details.
#' @seealso [pense] for an S-estimate of regression with elastic net penalty.
#' @export
elnet <- function(x, y, alpha, lambdas, penalty_loadings, weights,
include_intercept = TRUE, en_algorithm_opts) {
optional_args <- list()
# Normalize input
y <- .as(y, 'numeric')
x_dim <- dim(x)
if (length(y) != x_dim[[1L]]) {
stop("Number of observations does not match between `x` and `y`.")
}
alpha <- .as(alpha[[1L]], 'numeric')
lambdas <- sort(.as(lambdas, 'numeric'), decreasing = TRUE)
if (alpha < 0 || alpha > 1) {
stop("`alpha` is outside 0 and 1.")
}
if (alpha < sqrt(.Machine$double.eps)) {
alpha <- 0
}
penalty_loadings <- if (!missing(penalty_loadings)) {
if(alpha == 0) {
stop("Penalty loadings are only supported for alpha > 0.")
} else if (length(penalty_loadings) != x_dim[[2L]]) {
stop("Penalty loadings are not of length `p`.")
}
.as(penalty_loadings, 'numeric')
} else {
NULL
}
weights <- if (!missing(weights)) {
if (length(weights) != x_dim[[1L]]) {
stop("Observation weights are not the same length as `y`.")
}
.as(weights, 'numeric')
} else {
NULL
}
if (missing(en_algorithm_opts)) {
en_algorithm_opts <- NULL
}
# Check EN algorithm
optional_args$en_options <- .select_en_algorithm(en_algorithm_opts, alpha, x)
# Call internal function
res <- .elnet_internal(x, y, alpha, lambdas, penalty_loadings, weights,
include_intercept, optional_args)
if (!is.null(res$metrics)) {
class(res$metrics) <- 'nsoptim_metrics'
}
return(res)
}
#' Perform some final input adjustments and call the internal C++ code.
.elnet_internal <- function(x, y, alpha, lambdas, penalty_loadings = NULL,
weights = NULL, include_intercept = TRUE,
optional_args) {
# Create penalties-list, without sorting the lambda sequence
penalties <- lapply(lambdas, function (l) { list(lambda = l, alpha = alpha) })
include_intercept <- isTRUE(include_intercept)
if (!is.null(penalty_loadings)) {
optional_args$pen_loadings <- penalty_loadings
}
if (!is.null(weights)) {
optional_args$obs_weights <- weights
}
return(.Call(C_lsen_regression, x, y, penalties, include_intercept,
optional_args))
}
#' ENPY Initial Estimates
#'
#' @param x `n` by `p` matrix of numeric predictors.
#' @param y vector of response values of length `n`.
#' @param alpha value for the alpha parameter, the balance between L1 and L2
#' penalization.
#' @param lambdas a vector of positive values for the lambda parameter.
#' @param penalty_loadings a vector of positive penalty loadings
#' (a.k.a. weights) for different penalization of each coefficient.
#' @param include_intercept include an intercept in the model.
#' @param bdp the desired breakdown point of the estimator, between 0 and 0.5.
#' @param cc consistency constant for the scale estimate. By default, the scale estimate is made consistent for the
#' given breakdown point under the Normal model.
#' @param en_algorithm_opts options for the EN algorithm. See [en_algorithm_options] for details.
#' @param mscale_opts options for the M-scale estimation. See [mscale_algorithm_options] for details.
#' @export
enpy_initial_estimates <- function (x, y, alpha, lambdas, bdp = 0.25, cc, include_intercept = TRUE,
enpy_opts = enpy_options(), mscale_opts = mscale_algorithm_options()) {
alpha <- .as(alpha[[1L]], 'numeric')
if (missing(cc)) {
cc <- NULL
}
if (alpha < 0 || alpha > 1) {
stop("`alpha` is outside 0 and 1.")
}
if (alpha < sqrt(.Machine$double.eps)) {
alpha <- 0
}
penalties <- make_penalties(alpha, lambdas)
s_loss_params <- list(mscale = .full_mscale_algo_options(bdp = bdp, cc = cc, mscale_opts = mscale_opts),
intercept = include_intercept)
# Check EN algorithm for ENPY
enpy_opts$en_options <- .select_en_algorithm(enpy_opts$en_options, alpha, x)
res <- .Call(C_penpy, x, drop(y), penalties, s_loss_params, enpy_opts, list())
lapply(res, function (res) {
if (!is.null(res$metrics)) {
class(res$metrics) <- 'nsoptim_metrics'
}
return(res)
})
}
## Make a list of penalties and ensure that `alpha` and `lambdas` are of
## correct type and order.
make_penalties <- function (alpha, lambdas) {
lapply(sort(as.numeric(lambdas), decreasing = TRUE), function (lambda) {
list(alpha = alpha, lambda = lambda)
})
}
#' @useDynLib pense, .registration = TRUE
#' @importFrom Rcpp evalCpp
NULL
#' Compute the PENSE Regularization Path
#'
#' @section Strategies:
#' The function supports several different strategies to compute, and use given, starting points for the optimization.
#' They are all applied to the entire regularization path, i.e., for every value in the descending lambda grid.
#' * *Zero-based* (`'Z'`): Use the 0-vector to start the optimization at the largest lambda value.
#' For subsequent lambda values, use the optimum at the previous lambda value as starting point.
#' * *EN-PY (shared)* (`'Es'`): compute the EN-PY estimates at each of the lambda values in `enpy_lambdas`. Each of
#' these estimates is used as starting point at every lambda value in `lambdas`.
#' * *EN-PY (individual)* (`'Ei'`): compute the EN-PY estimates at each of the lambda values in `enpy_lambdas`.
#' These estimates are used as starting points at the corresponding value in `lambdas`. Only the best `nr_tracks`
#' optima are retained. At the following lambda value, the previous optima are used as starting points.
#' At the next smaller lambda which is both in `enpy_lambdas` and in `lambdas`, an additional `nr_tracks` optima are
#' computed by starting from the EN-PY estimates, however, only the best `nr_tracks` estimates (from using the
#' previous optima and the EN-PY estimates as starting points) are retained.
#' * *Other Starting Points (shared)* (`'Os'`): use the starting points given in `other_starts` at every lambda value.
#' Only used if `other_starts` is a list of estimates.
#' * *Other Starting Points (individual)* (`'Oi'`): use the starting points in `other_starts` at the corresponding
#' lambda values. Only used if `other_starts` is a list of estimates and if at least one of these estimates includes
#' the value of the regularization parameter `lambda`. The same "warm-start" strategy as for *EN-PY (individual)*
#' is used.
#'
#' For example, the strategies string `strategies = 'ZEiOs'` uses the strategies *Zero-based*, *EN-PY (individual)*,
#' and *Other Starting Points (shared)* to compute the optima along the regularization path. The strategies
#' don't have to be in any specific order; for example `strategies = 'EiZOs'` and `strategies = 'OsEiZ'` would also
#' select the same strategies as in the previous example.
#'
#' @param x `n` by `p` matrix of numeric predictors.
#' @param y vector of response values of length `n`.
#' @param alpha value for the alpha parameter, the balance between L1 and L2
#' penalization.
#' @param lambdas a vector of positive values for the lambda parameter.
#' @param penalty_loadings a vector of positive penalty loadings
#' (a.k.a. weights) for different penalization of each
#' coefficient.
#' @param enpy_lambdas a vector of lambda values at which *cold* initial
#' estimates are computed (see [enpy] for details).
#' @param strategies a string specifying the strategies used to get starting points for the optimization.
#' See *Strategies* for a list of supported strategies. The strategies can be in any order.
#' @param other_starts a list of other initial estimates to try only at the corresponding lambda values.
#' Each list item must contain the `intercept` and `beta` coefficients.
#' If the strategy *Other Starting Points (individual)* is requested, any starting point which has
#' the value for the regularization parameter `lambda` included (alongside `intercept` and `beta`),
#' will be used only at the specified `lambda` value.
#' @param include_intercept include an intercept in the model.
#' @param bdp the desired breakdown point of the estimator, between 0 and 0.5.
#' @param cc consistency constant for the scale estimate. By default, the scale estimate is made consistent for the
#' given breakdown point under the Normal model.
#' @param eps convergence tolerance for the algorithm.
#' @param explore_solutions number of optima to track for the *EN-PY (individual)* and
#' *Other Starting Points (individual)* regularization path strategies. Also controls how many
#' starting points for *EN-PY (individual)* and *Other Starting Points (shared)* are
#' computed to the full accuracy.
#' @param explore_tol relative numerical tolerance when exploring possible solutions.
#' @param max_solutions only retain up to `max_solutions` unique optima per lambda value.
#' @param comparison_tol numeric tolerance to determine if two optima are equal The comparison is first done
#' on the absolute difference in the value of the objective function at the optima.
#' If this is less than `comparison_tol`, two optima are deemed equal if the squared difference
#' of the intercepts is less than `comparison_tol` and the squared L_2 norm of the difference
#' vector is less than `comparison_tol`.
#' @param sparse use sparse coefficient vector.
#' @param algorithm_opts options for the PENSE algorithm. See [pense_algorithm_options] for details.
#' @param mscale_opts options for the M-scale estimation. See [mscale_algorithm_options] for details.
#' @param enpy_opts options for the ENPY initial estimates, created with the
#' [enpy_options] function. See [enpy_initial_estimates] for details.
#' @param cold_lambdas a deprecated alias for `enpy_lambdas`.
#' @param ... currently not used.
#'
#' @export
#' @importFrom lifecycle deprecate_soft
pense <- function(x, y, alpha, lambdas, enpy_lambdas, penalty_loadings, strategies = 'ZEsEiOsOi', other_starts,
include_intercept = TRUE, bdp = 0.25, cc, eps = 1e-6, algorithm_opts,
explore_solutions = 10, explore_tol = 0.1, max_solutions = 10, comparison_tol = sqrt(eps),
sparse = TRUE, mscale_opts = mscale_algorithm_options(), enpy_opts = enpy_options(),
...) {
optional_args <- list()
mc <- match.call()
if ('cold_lambdas' %in% names(mc)) {
deprecate_soft('1.1', 'pense(cold_lambdas = )', with = 'pense(enpy_lambdas =)')
if (missing(enpy_lambdas)) {
enpy_lambdas <- list(...)[['cold_lambdas']]
} else {
warning("Using values in `enpy_lambdas` for EN-PY estimator, ignoring `cold_lambdas`")
}
}
if (missing(enpy_lambdas)) {
enpy_lambdas <- numeric(0L)
}
# Normalize input
y <- .as(y, 'numeric')
x_dim <- dim(x)
if (length(y) != x_dim[[1L]]) {
stop("Number of observations in `x` and `y` does not match.")
}
if (missing(cc)) {
cc <- NULL
}
alpha <- .as(alpha[[1L]], 'numeric')
if (alpha < sqrt(.Machine$double.eps)) {
alpha <- 0
}
lambdas <- sort(.as(lambdas, 'numeric'), decreasing = TRUE)
sparse <- isTRUE(sparse)
# Split the `other_starts` into individual and shared starts.
if (!missing(other_starts)) {
# All starting points can be used for the "shared" strategy
optional_args$shared_starts <- if (inherits(other_starts, 'pense_estimate_list')) {