Commit afd345a0 authored by davidkep's avatar davidkep

add methods to make pense fits starting points

parent b6387a64
# Generated by roxygen2: do not edit by hand
S3method(as_starting_point,enpy_starting_points)
S3method(as_starting_point,pense_cvfit)
S3method(as_starting_point,pense_fit)
S3method(as_starting_point,starting_point)
S3method(c,starting_point)
S3method(c,starting_points)
......
......@@ -61,7 +61,7 @@ coef.pense_fit <- function (object, lambda, sparse = NULL, exact = deprecated(),
#' @param object PENSE with cross-validated hyper-parameters to extract coefficients from.
#' @param lambda either a string specifying which penalty level to use or a a single numeric value of the penalty
#' parameter. See details.
#' @param se_mult If `lambda = "se"`, the multiplier to the standard error.
#' @param se_mult If `lambda = "se"`, the multiple of standard errors to tolerate.
#' @param sparse should coefficients be returned as sparse or dense vectors? Defaults to the `sparse` argument
#' supplied to [pense_cv()]. Can also be set to `sparse = 'matrix'`, in which case a sparse matrix
#' is returned instead of a sparse vector.
......
......@@ -213,6 +213,7 @@ en_ridge_options <- function () {
}
#' @importFrom rlang warn
## also adds an element `sparse` to the returned list, which is the required sparsity parameter!
.select_en_algorithm <- function (en_options, alpha, sparse, eps) {
if (!is.null(en_options)) {
# Check if the EN algorithm can handle the given `alpha`
......
......@@ -41,12 +41,12 @@
#' \item{`lambda`}{the sequence of penalization parameters.}
#' \item{`estimates`}{a list of estimates. Each estimate contains the following information:
#' \describe{
#' \item{`intercept`}{intercept estimate}
#' \item{`beta`}{beta (slope) estimate}
#' \item{`lambda`}{penalization level at which the estimate is computed}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate}
#' \item{`status`}{optional status message from the algorithm}
#' \item{`intercept`}{intercept estimate.}
#' \item{`beta`}{beta (slope) estimate.}
#' \item{`lambda`}{penalization level at which the estimate is computed.}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed.}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate.}
#' \item{`status`}{optional status message from the algorithm.}
#' }
#' }
#' \item{`call`}{the original call.}
......@@ -123,7 +123,8 @@ elnet <- function(x, y, alpha, nlambda = 100, lambda_min_ratio, lambda, penalty_
#'
#' @return a list with components:
#' \describe{
#' \item{`cvres`}{data frame of average cross-validated performance}
#' \item{`lambda`}{the sequence of penalization levels.}
#' \item{`cvres`}{data frame of average cross-validated performance.}
#' \item{`cv_replications`}{matrix of cross-validated performance metrics, one column per replication.
#' Rows are in the same order as in `cvres`.}
#' \item{`call`}{the original call.}
......@@ -188,7 +189,8 @@ elnet_cv <- function (x, y, lambda, cv_k, cv_repl = 1, cv_metric = c('rmspe', 't
args$restore_coef_length(args$std_data$unstandardize_coefs(est))
})
return(structure(list(call = call, cvres = cv_perf_df, cv_replications = cv_perf, cv_measure = cv_measure_str,
estimates = .metrics_attrib(fit$estimates, fit$metrics)), class = c('pense_en', 'pense_cvfit')))
lambda = args$lambda, estimates = .metrics_attrib(fit$estimates, fit$metrics)),
class = c('pense_en', 'pense_cvfit')))
}
## Perform some final input adjustments and call the internal C++ code.
......
......@@ -54,6 +54,7 @@
#' Note that if a the starting point is *specific* to a penalization level, this penalization
#' level is added to the grid of penalization levels (either the manually specified grid in `lambda`
#' or the automatically generated grid of size `nlambda`).
#' If `standardize = TRUE`, the starting points are also scaled.
#' @param standardize logical flag to standardize the `x` variables prior to fitting the PENSE estimates.
#' Coefficients are always returned on the original scale. This can fail for variables with a large
#' proportion of a single value (e.g., zero-inflated data). In this case, either compute with
......@@ -84,16 +85,16 @@
#'
#' @return a list-like object with the following items
#' \describe{
#' \item{`lambda`}{the sequence of penalization levels}
#' \item{`lambda`}{the sequence of penalization levels.}
#' \item{`estimates`}{a list of estimates. Each estimate contains the following information:
#' \describe{
#' \item{`intercept`}{intercept estimate}
#' \item{`beta`}{beta (slope) estimate}
#' \item{`lambda`}{penalization level at which the estimate is computed}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed}
#' \item{`objf_value`}{value of the objective function at the solution}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate}
#' \item{`status`}{optional status message from the algorithm}
#' \item{`intercept`}{intercept estimate.}
#' \item{`beta`}{beta (slope) estimate.}
#' \item{`lambda`}{penalization level at which the estimate is computed.}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed.}
#' \item{`objf_value`}{value of the objective function at the solution.}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate.}
#' \item{`status`}{optional status message from the algorithm.}
#' }
#' }
#' \item{`call`}{the original call.}
......@@ -175,7 +176,8 @@ pense <- function(x, y, alpha, nlambda = 50, nlambda_enpy = 10, lambda, lambda_m
#'
#' @return a list with components:
#' \describe{
#' \item{`cvres`}{data frame of average cross-validated performance}
#' \item{`lambda`}{the sequence of penalization levels.}
#' \item{`cvres`}{data frame of average cross-validated performance.}
#' \item{`cv_replications`}{matrix of cross-validated performance metrics, one column per replication.
#' Rows are in the same order as in `cvres`.}
#' \item{`call`}{the original call.}
......@@ -243,8 +245,9 @@ pense_cv <- function(x, y, standardize = TRUE, lambda, cv_k, cv_repl = 1,
fit$estimates <- lapply(fit$estimates, function (ests) {
args$restore_coef_length(args$std_data$unstandardize_coefs(ests[[1L]]))
})
return(structure(list(call = call, cvres = cv_perf_df, cv_replications = cv_perf, cv_measure = cv_measure_str,
estimates = .metrics_attrib(fit$estimates, fit$metrics)), class = c('pense', 'pense_cvfit')))
return(structure(list(call = call, lambda = args$lambda, cvres = cv_perf_df, cv_replications = cv_perf,
cv_measure = cv_measure_str, estimates = .metrics_attrib(fit$estimates, fit$metrics)),
class = c('pense', 'pense_cvfit')))
}
#' @description `adapense_cv()` is a convenience wrapper to compute adaptive PENSE estimates.
......@@ -478,7 +481,6 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
nlambda_enpy <- length(enpy_lambda)
}
sparse <- isTRUE(sparse)
pense_opts <- list(algo_opts = algorithm_opts,
strategy_0 = isTRUE(add_zero_based),
strategy_enpy_individual = isTRUE(enpy_specific) && (nlambda_enpy > 0L),
......@@ -493,6 +495,7 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
nr_tracks = .as(explore_solutions[[1L]], 'integer'),
max_optima = .as(max_solutions[[1L]], 'integer'),
num_threads = max(1L, .as(ncores[[1L]], 'integer')),
sparse = isTRUE(sparse),
mscale = .full_mscale_algo_options(bdp = bdp, mscale_opts = mscale_opts))
if (pense_opts$explore_tol < pense_opts$eps) {
......@@ -503,16 +506,17 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
}
# Check EN algorithm for ENPY
enpy_opts$en_options <- .select_en_algorithm(enpy_opts$en_options, alpha, sparse, eps)
enpy_opts$en_options <- .select_en_algorithm(enpy_opts$en_options, alpha, pense_opts$sparse, eps)
pense_opts$sparse <- enpy_opts$en_options$sparse
# If using the MM algorithm, ensure that the EN options are set.
if (pense_opts$algorithm == 1L) {
pense_opts$algo_opts$en_options <- .select_en_algorithm(pense_opts$algo_opts$en_options, alpha, sparse, eps)
if (!isTRUE(enpy_opts$en_options$sparse == pense_opts$algo_opts$en_options$sparse)) {
pense_opts$algo_opts$en_options <- .select_en_algorithm(pense_opts$algo_opts$en_options, alpha, pense_opts$sparse,
eps)
if (!isTRUE(pense_opts$sparse == pense_opts$algo_opts$en_options$sparse)) {
abort("The `sparse` option for the EN-PY algorithm and the MM algorithm for PENSE disagree.")
}
}
sparse <- enpy_opts$en_options$sparse
# Set the number of cores for the ENPY options
if (pense_opts$num_threads > 1L && !isTRUE(.k_multithreading_support)) {
......@@ -533,7 +537,7 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
# Check penalty loadings
if (!missing(penalty_loadings) && !is.null(penalty_loadings)) {
checked_pls <- .prepare_penalty_loadings(penalty_loadings, x, alpha, sparse = sparse)
checked_pls <- .prepare_penalty_loadings(penalty_loadings, x, alpha, sparse = pense_opts$sparse)
penalty_loadings <- checked_pls$loadings
restore_coef_length <- checked_pls$restore_fun
x <- checked_pls$trimmed_x
......@@ -546,8 +550,9 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
pense_opts$intercept <- TRUE
warn("All values in `penalty_loadings` are infinite. Only computing the intercept.")
std_data <- .standardize_data(matrix(runif(x_dim[[1L]]), ncol = 1L), y, intercept = TRUE,
sparse = sparse, standardize = standardize, robust = TRUE, mscale_opts = mscale_opts,
bdp = pense_opts$mscale$delta, scale_cc = pense_opts$mscale$cc)
sparse = pense_opts$sparse, standardize = standardize, robust = TRUE,
mscale_opts = mscale_opts, bdp = pense_opts$mscale$delta,
scale_cc = pense_opts$mscale$cc)
# Compute only the 0-based solution.
pense_opts$strategy_enpy_individual <- FALSE
pense_opts$strategy_enpy_shared <- FALSE
......@@ -560,7 +565,7 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
}
std_data <- .standardize_data(x, y, intercept = pense_opts$intercept, standardize = standardize, robust = TRUE,
sparse = sparse, mscale_opts = mscale_opts, bdp = pense_opts$mscale$delta,
sparse = pense_opts$sparse, mscale_opts = mscale_opts, bdp = pense_opts$mscale$delta,
scale_cc = pense_opts$mscale$cc)
# Scale penalty loadings appropriately
......@@ -590,7 +595,9 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
}
# Ensure the `beta` coefficients in `other_starts` agree with the desired vector class (sparse vs. dense)
other_starts <- .sparsify_other_starts(other_starts, sparse)
# and standardize them.
other_starts <- lapply(.sparsify_other_starts(other_starts, pense_opts$sparse),
std_data$standardize_coefs)
# Identify which other starts are shared and which are specific.
other_starts_shared <- unlist(lapply(other_starts, is, 'shared_starting_point'), recursive = FALSE,
......@@ -604,7 +611,7 @@ adapense_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...)
}
if (any(other_starts_specific)) {
pense_opts$strategy_other_individual <- TRUE
ind_starts <- .make_initest_list(other_starts[other_starts_specific], lambda, sparse = sparse)
ind_starts <- .make_initest_list(other_starts[other_starts_specific], lambda, sparse = pense_opts$sparse)
optional_args$individual_starts <- ind_starts$starting_points
lambda <- ind_starts$extended_lambda
}
......
......@@ -44,16 +44,16 @@
#'
#' @return a list-like object with the following items
#' \describe{
#' \item{`lambda`}{the sequence of penalization levels}
#' \item{`lambda`}{the sequence of penalization levels.}
#' \item{`estimates`}{a list of estimates. Each estimate contains the following information:
#' \describe{
#' \item{`intercept`}{intercept estimate}
#' \item{`beta`}{beta (slope) estimate}
#' \item{`lambda`}{penalization level at which the estimate is computed}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed}
#' \item{`objf_value`}{value of the objective function at the solution}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate}
#' \item{`status`}{optional status message from the algorithm}
#' \item{`intercept`}{intercept estimate.}
#' \item{`beta`}{beta (slope) estimate.}
#' \item{`lambda`}{penalization level at which the estimate is computed.}
#' \item{`alpha`}{*alpha* hyper-parameter at which the estimate is computed.}
#' \item{`objf_value`}{value of the objective function at the solution.}
#' \item{`statuscode`}{if `> 0` the algorithm experienced issues when computing the estimate.}
#' \item{`status`}{optional status message from the algorithm.}
#' }
#' }
#' \item{`call`}{the original call.}
......@@ -122,7 +122,8 @@ regmest <- function(x, y, alpha, nlambda = 50, lambda, lambda_min_ratio, scale,
#'
#' @return a list with components:
#' \describe{
#' \item{`cvres`}{data frame of average cross-validated performance}
#' \item{`lambda`}{the sequence of penalization levels.}
#' \item{`cvres`}{data frame of average cross-validated performance.}
#' \item{`cv_replications`}{matrix of cross-validated performance metrics, one column per replication.
#' Rows are in the same order as in `cvres`.}
#' \item{`call`}{the original call.}
......@@ -189,7 +190,8 @@ regmest_cv <- function(x, y, standardize = TRUE, lambda, cv_k, cv_repl = 1,
})
return(structure(list(call = call, cvres = cv_perf_df, cv_replications = cv_perf, cv_measure = cv_measure_str,
estimates = .metrics_attrib(fit$estimates, fit$metrics)), class = c('mest', 'pense_cvfit')))
lambda = args$lambda, estimates = .metrics_attrib(fit$estimates, fit$metrics)),
class = c('mest', 'pense_cvfit')))
}
#' @description `adamest_cv()` is a convenience wrapper to compute adaptive elastic-net M-estimates.
......@@ -238,51 +240,6 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
return(adamstep)
}
## Make a list of initial estimates
##
## @return a list with 2 components:
## `extended_lambda` an extended grid of penalization levels to contain both the given lambda values
## plus the lambda values in `other_starts`
## `starting_points` a list the same length as `extended_lambda` with a list of initial estimates
## for each value in `extended_lambda`.
#' @importFrom Matrix sparseVector
#' @importClassesFrom Matrix dsparseVector
#' @importFrom rlang warn
.make_initest_list <- function (other_starts, lambda, sparse) {
if (length(other_starts) == 0L) {
return(list(extended_lambda = lambda, starting_points = rep.int(list(list()), length(lambda))))
}
# Check for impostor starting points without lambda
init_est_lambda <- unlist(lapply(other_starts, `[[`, 'lambda'), use.names = FALSE, recursive = FALSE)
if (length(init_est_lambda) != length(other_starts)) {
abort("Some starting points in `other_starts` are marked as \"specific\" but do not have a `lambda` component.")
}
init_est_lambda <- .as(init_est_lambda, 'numeric')
init_est_inds <- .approx_match(init_est_lambda, lambda)
new_initest_lambda <- which(is.na(init_est_inds))
if (length(new_initest_lambda) > 0L) {
# Some starts are for unknown lambda. Add lambdas to the grid!
lambda <- sort(c(lambda, unique(init_est_lambda[new_initest_lambda])), decreasing = TRUE)
init_est_inds <- .approx_match(init_est_lambda, lambda)
}
starting_points <- lapply(seq_along(lambda), function (i) {
matches <- which(i == init_est_inds)
if (length(matches) > 0L) {
return(other_starts[matches])
} else {
return(list())
}
})
return(list(extended_lambda = lambda, starting_points = starting_points))
}
## Get the smallest lambda such that the PENSE estimate gives the empty model.
.regmest_max_lambda <- function (x, y, alpha, scale, mest_options, penalty_loadings = NULL) {
optional_args <- list()
......@@ -350,12 +307,12 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
mscale_opts <- .full_mscale_algo_options(bdp = mscale_bdp, mscale_opts = mscale_opts)
sparse <- isTRUE(sparse)
mest_opts <- list(algo_opts = algorithm_opts,
cc = .as(cc[[1L]], 'numeric'),
strategy_0 = isTRUE(add_zero_based),
algorithm = .regmest_algorithm_id(algorithm_opts),
intercept = !isFALSE(intercept),
sparse = isTRUE(sparse),
eps = .as(eps[[1L]], 'numeric'),
max_it = mscale_opts$max_it, # max. iterations for finding the largest lambda
comparison_tol = .as(comparison_tol[[1L]], 'numeric'),
......@@ -377,8 +334,8 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
# If using the MM algorithm, ensure that the EN options are set.
if (mest_opts$algorithm == 1L) {
mest_opts$algo_opts$en_options <- .select_en_algorithm(mest_opts$algo_opts$en_options, alpha, sparse, eps)
sparse <- mest_opts$algo_opts$en_options$sparse
mest_opts$algo_opts$en_options <- .select_en_algorithm(mest_opts$algo_opts$en_options, alpha, mest_opts$sparse, eps)
mest_opts$sparse <- mest_opts$algo_opts$en_options$sparse
}
# Set the number of cores for the ENPY options
......@@ -399,7 +356,7 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
# Check penalty loadings
if (!missing(penalty_loadings) && !is.null(penalty_loadings)) {
checked_pls <- .prepare_penalty_loadings(penalty_loadings, x, alpha, sparse = sparse)
checked_pls <- .prepare_penalty_loadings(penalty_loadings, x, alpha, sparse = mest_opts$sparse)
penalty_loadings <- checked_pls$loadings
restore_coef_length <- checked_pls$restore_fun
x <- checked_pls$trimmed_x
......@@ -411,7 +368,7 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
if (ncol(x) == 0L) {
warn("All values in `penalty_loadings` are infinite. Only computing the intercept.")
mest_opts$intercept <- TRUE
std_data <- .standardize_data(matrix(runif(x_dim[[1L]]), ncol = 1L), y, intercept = TRUE, sparse = sparse,
std_data <- .standardize_data(matrix(runif(x_dim[[1L]]), ncol = 1L), y, intercept = TRUE, sparse = mest_opts$sparse,
standardize = standardize, robust = TRUE, mscale_opts = mscale_opts,
bdp = mscale_opts$delta, scale_cc = mscale_opts$cc)
# Compute only the 0-based solution.
......@@ -423,7 +380,7 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
}
std_data <- .standardize_data(x, y, intercept = mest_opts$intercept, standardize = standardize, robust = TRUE,
sparse = sparse, mscale_opts = mscale_opts, bdp = mscale_opts$delta,
sparse = mest_opts$sparse, mscale_opts = mscale_opts, bdp = mscale_opts$delta,
scale_cc = mscale_opts$cc)
# Scale penalty loadings appropriately
......@@ -444,7 +401,9 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
nlambda <- length(lambda)
if (!missing(starting_points)) {
if (is(starting_points, 'starting_point')) {
if (is(starting_points, 'pense_fit') || is(starting_points, 'pense_cvfit')) {
starting_points <- as_starting_point(starting_points)
} else if (is(starting_points, 'starting_point')) {
starting_points <- structure(list(starting_points), class = 'starting_points')
} else if (!is(starting_points, 'starting_points')) {
abort(paste("`starting_points` must be a list of starting points created by",
......@@ -452,7 +411,9 @@ adamest_cv <- function (x, y, alpha, alpha_preliminary = 0, exponent = 1, ...) {
}
# Ensure the `beta` coefficients in `other_starts` agree with the desired vector class (sparse vs. dense)
optional_args$shared_starts <- .sparsify_other_starts(starting_points, sparse)
# and standardize them.
optional_args$shared_starts <- lapply(.sparsify_other_starts(starting_points, mest_opts$sparse),
std_data$standardize_coefs)
}
if (any(lambda < .Machine$double.eps)) {
......
......@@ -385,6 +385,7 @@ extract_metric <- function (metrics, attr, node) {
} else if (!isTRUE(sparse) && !is.numeric(est$beta)) {
est$beta <- .as(est$beta, 'numeric')
}
class(est) <- NULL
return(est)
})
}
......@@ -226,6 +226,89 @@ as_starting_point.enpy_starting_points <- function (object, specific = FALSE, ..
}
}
#' @rdname starting_point
#' @param lambda optional penalization level(s) to extract from `object`. Penalization levels not present in `object`
#' are ignored with a warning.
#' @export
#' @importFrom rlang warn abort
as_starting_point.pense_fit <- function (object, specific = FALSE, lambda, ...) {
lambda_indices <- if (missing(lambda)) {
seq_along(object$estimates)
} else if (is.numeric(lambda)) {
lambda_indices <- .approx_match(lambda, object$lambda)
bad_lambda_indices <- which(is.na(lambda_indices))
if (length(bad_lambda_indices) == length(lambda)) {
abort("No penalization level requested with `lambda` is available in `object`.")
} else if (length(bad_lambda_indices) > 0L) {
warn("Some penalization levels requested with `lambda` are not available in `object`.")
lambda_indices[-bad_lambda_indices]
} else {
lambda_indices
}
} else {
abort("`lambda` must be missing or a numeric vector.")
}
class_list <- if (isTRUE(specific)) {
c('specific_starting_point', 'starting_point')
} else {
c('shared_starting_point', 'starting_point')
}
structure(lapply(object$estimates[lambda_indices], structure, class = class_list), class = 'starting_points')
}
#' @rdname starting_point
#' @param lambda optionally either a string specifying which penalty level to use (`"min"` or `"se"`) or a numeric
#' vector of the penalty levels to extract from `object`. Penalization levels not present in `object`
#' are ignored with a warning. If `NULL`, all estimates in `object` are extracted.
#' @param se_mult If `lambda = "se"`, the multiple of standard errors to tolerate.
#'
#' @details
#' When creating starting points from cross-validated fits, it is possible to extract only the estimate with best
#' CV performance (`lambda = "min"`), or the estimate with CV performance statistically indistinguishable from the
#' best performance (`lambda = "se"`). This is determined to be the estimate with prediction performance
#' within `se_mult * cv_se` from the best model.
#'
#' @export
#' @importFrom rlang warn abort
as_starting_point.pense_cvfit <- function (object, specific = FALSE, lambda = c('min', 'se'), se_mult = 2, ...) {
lambda_indices <- if (is.character(lambda)) {
lambda <- match.arg(lambda)
se_mult <- if (lambda == 'min') {
0
} else {
.as(se_mult[[1L]], 'numeric')
}
se_selection <- .cv_se_selection(object$cvres$cvavg, object$cvres$cvse, se_mult)
which(se_selection == 'se_fact')
} else if (is.numeric(lambda)) {
lambda_indices <- .approx_match(lambda, object$lambda)
bad_lambda_indices <- which(is.na(lambda_indices))
if (length(bad_lambda_indices) == length(lambda)) {
abort("No penalization level requested with `lambda` is available in `object`.")
} else if (length(bad_lambda_indices) > 0L) {
warn("Some penalization levels requested with `lambda` are not available in `object`.")
lambda_indices[-bad_lambda_indices]
} else {
lambda_indices
}
} else if (is.null(lambda)) {
seq_along(object$estimates)
} else {
abort("`lambda` must be missing, a numeric vector, or either one of \"min\" or \"se\".")
}
class_list <- if (isTRUE(specific)) {
c('specific_starting_point', 'starting_point')
} else {
c('shared_starting_point', 'starting_point')
}
structure(lapply(object$estimates[lambda_indices], structure, class = class_list), class = 'starting_points')
}
#' @export
#' @keywords internal
#' @importFrom rlang abort
......
......@@ -20,7 +20,7 @@
\item{lambda}{either a string specifying which penalty level to use or a a single numeric value of the penalty
parameter. See details.}
\item{se_mult}{If \code{lambda = "se"}, the multiplier to the standard error.}
\item{se_mult}{If \code{lambda = "se"}, the multiple of standard errors to tolerate.}
\item{sparse}{should coefficients be returned as sparse or dense vectors? Defaults to the \code{sparse} argument
supplied to \code{\link[=pense_cv]{pense_cv()}}. Can also be set to \code{sparse = 'matrix'}, in which case a sparse matrix
......
......@@ -71,12 +71,12 @@ a list-like object with the following items
\item{\code{lambda}}{the sequence of penalization parameters.}
\item{\code{estimates}}{a list of estimates. Each estimate contains the following information:
\describe{
\item{\code{intercept}}{intercept estimate}
\item{\code{beta}}{beta (slope) estimate}
\item{\code{lambda}}{penalization level at which the estimate is computed}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate}
\item{\code{status}}{optional status message from the algorithm}
\item{\code{intercept}}{intercept estimate.}
\item{\code{beta}}{beta (slope) estimate.}
\item{\code{lambda}}{penalization level at which the estimate is computed.}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed.}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate.}
\item{\code{status}}{optional status message from the algorithm.}
}
}
\item{\code{call}}{the original call.}
......
......@@ -38,7 +38,8 @@ level with smallest average CV performance.}
\value{
a list with components:
\describe{
\item{\code{cvres}}{data frame of average cross-validated performance}
\item{\code{lambda}}{the sequence of penalization levels.}
\item{\code{cvres}}{data frame of average cross-validated performance.}
\item{\code{cv_replications}}{matrix of cross-validated performance metrics, one column per replication.
Rows are in the same order as in \code{cvres}.}
\item{\code{call}}{the original call.}
......
......@@ -76,7 +76,8 @@ If the output of \code{\link[=enpy_initial_estimates]{enpy_initial_estimates()}}
among all penalization levels.
Note that if a the starting point is \emph{specific} to a penalization level, this penalization
level is added to the grid of penalization levels (either the manually specified grid in \code{lambda}
or the automatically generated grid of size \code{nlambda}).}
or the automatically generated grid of size \code{nlambda}).
If \code{standardize = TRUE}, the starting points are also scaled.}
\item{eps}{numerical tolerance.}
......@@ -115,16 +116,16 @@ proportion of a single value (e.g., zero-inflated data). In this case, either co
\value{
a list-like object with the following items
\describe{
\item{\code{lambda}}{the sequence of penalization levels}
\item{\code{lambda}}{the sequence of penalization levels.}
\item{\code{estimates}}{a list of estimates. Each estimate contains the following information:
\describe{
\item{\code{intercept}}{intercept estimate}
\item{\code{beta}}{beta (slope) estimate}
\item{\code{lambda}}{penalization level at which the estimate is computed}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed}
\item{\code{objf_value}}{value of the objective function at the solution}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate}
\item{\code{status}}{optional status message from the algorithm}
\item{\code{intercept}}{intercept estimate.}
\item{\code{beta}}{beta (slope) estimate.}
\item{\code{lambda}}{penalization level at which the estimate is computed.}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed.}
\item{\code{objf_value}}{value of the objective function at the solution.}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate.}
\item{\code{status}}{optional status message from the algorithm.}
}
}
\item{\code{call}}{the original call.}
......
......@@ -53,7 +53,8 @@ level with smallest average CV performance.}
\value{
a list with components:
\describe{
\item{\code{cvres}}{data frame of average cross-validated performance}
\item{\code{lambda}}{the sequence of penalization levels.}
\item{\code{cvres}}{data frame of average cross-validated performance.}
\item{\code{cv_replications}}{matrix of cross-validated performance metrics, one column per replication.
Rows are in the same order as in \code{cvres}.}
\item{\code{call}}{the original call.}
......
......@@ -93,16 +93,16 @@ proportion of a single value (e.g., zero-inflated data). In this case, either co
\value{
a list-like object with the following items
\describe{
\item{\code{lambda}}{the sequence of penalization levels}
\item{\code{lambda}}{the sequence of penalization levels.}
\item{\code{estimates}}{a list of estimates. Each estimate contains the following information:
\describe{
\item{\code{intercept}}{intercept estimate}
\item{\code{beta}}{beta (slope) estimate}
\item{\code{lambda}}{penalization level at which the estimate is computed}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed}
\item{\code{objf_value}}{value of the objective function at the solution}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate}
\item{\code{status}}{optional status message from the algorithm}
\item{\code{intercept}}{intercept estimate.}
\item{\code{beta}}{beta (slope) estimate.}
\item{\code{lambda}}{penalization level at which the estimate is computed.}
\item{\code{alpha}}{\emph{alpha} hyper-parameter at which the estimate is computed.}
\item{\code{objf_value}}{value of the objective function at the solution.}
\item{\code{statuscode}}{if \verb{> 0} the algorithm experienced issues when computing the estimate.}
\item{\code{status}}{optional status message from the algorithm.}
}
}
\item{\code{call}}{the original call.}
......
......@@ -53,7 +53,8 @@ level with smallest average CV performance.}
\value{
a list with components:
\describe{
\item{\code{cvres}}{data frame of average cross-validated performance}
\item{\code{lambda}}{the sequence of penalization levels.}
\item{\code{cvres}}{data frame of average cross-validated performance.}
\item{\code{cv_replications}}{matrix of cross-validated performance metrics, one column per replication.
Rows are in the same order as in \code{cvres}.}
\item{\code{call}}{the original call.}
......
......@@ -4,6 +4,8 @@
\alias{starting_point}
\alias{as_starting_point}
\alias{as_starting_point.enpy_starting_points}
\alias{as_starting_point.pense_fit}
\alias{as_starting_point.pense_cvfit}
\title{Create Starting Points for the PENSE Algorithm}
\usage{
starting_point(beta, intercept, lambda)
......@@ -11,6 +13,16 @@ starting_point(beta, intercept, lambda)
as_starting_point(object, specific = FALSE, ...)
\method{as_starting_point}{enpy_starting_points}(object, specific = FALSE, ...)
\method{as_starting_point}{pense_fit}(object, specific = FALSE, lambda, ...)
\method{as_starting_point}{pense_cvfit}(
object,
specific = FALSE,
lambda = c("min", "se"),
se_mult = 2,
...
)
}
\arguments{
\item{beta}{beta coefficients at the starting point. Can be a numeric vector, a sparse vector of class
......@@ -19,8 +31,11 @@ as_starting_point(object, specific = FALSE, ...)
\item{intercept}{intercept coefficient at the starting point.}
\item{lambda}{optional penalization level. If missing, the starting point is used at \emph{every} penalization level,
otherwise, only at this penalization level.}
\item{lambda}{optionally either a string specifying which penalty level to use (\code{"min"} or \code{"se"}) or a numeric
vector of the penalty levels to extract from \code{object}. Penalization levels not present in \code{object}
are ignored with a warning. If \code{NULL}, all estimates in \code{object} are extracted.}
\item{se_mult}{If \code{lambda = "se"}, the multiple of standard errors to tolerate.}
}
\value{
an object of type \code{starting_points} to be used as starting point for \code{\link[=pense]{pense()}}.
......@@ -42,4 +57,9 @@ See \code{\link[=pense]{pense()}} for details.
Starting points computed via \code{\link[=enpy_initial_estimates]{enpy_initial_estimates()}} are by default \emph{shared} starting points but can
be transformed to \emph{specific} starting points via \code{enpy_starting_point(..., specific = TRUE)}.
When creating starting points from cross-validated fits, it is possible to extract only the estimate with best
CV performance (\code{lambda = "min"}), or the estimate with CV performance statistically indistinguishable from the
best performance (\code{lambda = "se"}). This is determined to be the estimate with prediction performance
within \code{se_mult * cv_se} from the best model.
}
......@@ -52,9 +52,10 @@ constexpr double kDefaultBisquareMscaleCc = 2.937015;
//! Default breakdown point for the M-scale equation.
constexpr double kDefaultMscaleDelta = 0.25;
constexpr EnAlgorithm kDefaultEnAlgorithm = EnAlgorithm::kLinearizedAdmm;
constexpr EnAlgorithm kDefaultEnAlgorithm = EnAlgorithm::kLars;
constexpr PenseAlgorithm kDefaultPenseAlgorithm = PenseAlgorithm::kMm;
constexpr bool kDefaultUseSparse = true;
constexpr MestEnAlgorithm kDefaultMestAlgorithm = MestEnAlgorithm::kMm;
constexpr bool kDefaultUseSparse = false;
} // namespace pense
......
......@@ -82,16 +82,17 @@ template<typename PenaltyFunction>
SEXP MMAlgorithmDispatch(SEXP x, SEXP y, SEXP scale, SEXP penalties, SEXP r_mest_opts,
const Rcpp::List& optional_args) {
using pense::GetFallback;
using Rcpp::List;
using SparseCoefs = nsoptim::RegressionCoefficients<arma::sp_vec>;
using DenseCoefs = nsoptim::RegressionCoefficients<arma::vec>;
using pense::r_interface::MakeOptimizer;
using MLoss = pense::MLoss<pense::RhoBisquare>;
using SurrogateLoss = typename MLoss::ConvexSurrogateType;
const auto mest_opts = as<Rcpp::List>(r_mest_opts);
const auto mm_options = GetFallback(mest_opts, "algo_opts", Rcpp::List());
const auto en_options = GetFallback(mm_options, "en_options", Rcpp::List());
const bool use_sparse_coefs = GetFallback(mm_options, "sparse", pense::kDefaultUseSparse);
const auto mest_opts = as<List>(r_mest_opts);
const auto mm_options = GetFallback(mest_opts, "algo_opts", List());
const auto en_options = GetFallback(mm_options, "en_options", List());
const bool use_sparse_coefs = GetFallback(mest_opts, "sparse", pense::kDefaultUseSparse);
switch (GetFallback(en_options, "algorithm", pense::kDefaultEnAlgorithm)) {
case pense::EnAlgorithm::kDal: {
......
......@@ -366,7 +366,7 @@ SEXP PenseMMDispatch(SEXP x, SEXP y, SEXP penalties, SEXP enpy_inds, const Rcpp:
using SurrogateLoss = typename SLoss::ConvexSurrogateType;
const auto mm_options = GetFallback(pense_opts, "algo_opts", Rcpp::List());
const auto en_options = GetFallback(mm_options, "en_options", Rcpp::List());
const bool use_sparse_coefs = GetFallback(mm_options, "sparse", pense::kDefaultUseSparse);
const bool use_sparse_coefs = GetFallback(pense_opts, "sparse", pense::kDefaultUseSparse);
switch (GetFallback(en_options, "algorithm", pense::kDefaultEnAlgorithm)) {
case pense::EnAlgorithm::kDal: {
......@@ -412,14 +412,14 @@ template<typename PenaltyFunction>
SEXP PenseAdmmDispatch(SEXP x, SEXP y, SEXP penalties, SEXP enpy_inds, const Rcpp::List& pense_opts,
SEXP enpy_opts, const Rcpp::List& optional_args) {
const auto admm_options = GetFallback(pense_opts, "algo_opts", Rcpp::List());
const bool use_sparse_coefs = GetFallback(admm_options, "sparse", pense::kDefaultUseSparse);
const bool use_sparse_coefs = GetFallback(pense_opts, "sparse", pense::kDefaultUseSparse);
if (use_sparse_coefs) {
using Optimizer = LinearizedAdmmOptimizer<pense::PenseProximalOperator, PenaltyFunction, SparseCoefs>;
return PenseRegressionImpl(MakeOptimizer<Optimizer>(admm_options), x, y, penalties, enpy_inds,
pense_opts, enpy_opts, optional_args);
} else {
using Optimizer = LinearizedAdmmOptimizer<pense::PenseProximalOperator, PenaltyFunction, SparseCoefs>;
using Optimizer = LinearizedAdmmOptimizer<pense::PenseProximalOperator, PenaltyFunction, DenseCoefs>;
return PenseRegressionImpl(MakeOptimizer<Optimizer>(admm_options), x, y, penalties, enpy_inds,
pense_opts, enpy_opts, optional_args);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment