Commit 78fdcd09 authored by davidkep's avatar davidkep

add CV for elnet and coef methods

parent 8181a910
......@@ -21,7 +21,7 @@ Description: Robust penalized elastic net S and MM estimator for linear
Cohen Freue, G. V., Kepplinger, D., Salibian-Barrera, M., and Smucler, E.
(2017) <https://projecteuclid.org/euclid.aoas/1574910036>.
Depends:
R (>= 3.4.0),
R (>= 3.5.0),
Matrix
Imports:
Rcpp,
......
# Generated by roxygen2: do not edit by hand
S3method(coef,elnet_cv_res)
S3method(coef,elnet_res)
S3method(print,nsoptim_metrics)
export(consistency_const)
export(elnet)
export(elnet_cv)
export(en_admm_options)
export(en_dal_options)
export(en_lars_options)
......@@ -23,11 +26,14 @@ export(pense_options)
export(prinsens)
export(rho_function)
export(tau_size)
importFrom(Matrix,drop)
importFrom(Matrix,sparseVector)
importFrom(Rcpp,evalCpp)
importFrom(lifecycle,deprecate_stop)
importFrom(lifecycle,deprecate_warn)
importFrom(lifecycle,deprecated)
importFrom(lifecycle,is_present)
importFrom(methods,as)
importFrom(methods,is)
importFrom(stats,mad)
useDynLib(pense, .registration = TRUE)
This diff is collapsed.
......@@ -270,7 +270,7 @@ pense <- function(x, y, alpha, nlambda = 50, lambda_min_ratio, nlambda_enpy = 10
.Call(C_pense_max_lambda, x, y, pense_options, optional_args) / max(0.01, alpha)
}
#' Generate a log-spaced grid of decreasing lambda values
## Generate a log-spaced grid of decreasing lambda values
.pense_lambda_grid <- function (x, y, alpha, nlambda, lambda_min_ratio, pense_options, penalty_loadings = NULL) {
alpha <- max(0.01, alpha)
x_dim <- dim(x)
......@@ -282,7 +282,7 @@ pense <- function(x, y, alpha, nlambda = 50, lambda_min_ratio, nlambda_enpy = 10
}
}
max_lambda <- .pense_max_lambda(x, y, alpha, pense_options, penalty_loadings)
lambdas <- rev(exp(seq(log(lambda_min_ratio * max_lambda), log(max_lambda), length.out = nlambda)))
rev(exp(seq(log(lambda_min_ratio * max_lambda), log(max_lambda), length.out = nlambda)))
}
## Perform some final input adjustments and call the internal C++ code.
......
......@@ -150,7 +150,7 @@ rho_function <- function (rho) {
##
## @param x,table see [base::match] for details.
## @param eps numerical tolerance for matching.
## @return a vector the same lenght of `x` with integers giving the position in
## @return a vector the same length as `x` with integers giving the position in
## `table` of the first match if there is a match, or `NA_integer_`
## otherwise.
.approx_match <- function(x, table,
......@@ -188,3 +188,166 @@ extract_metric <- function (metrics, attr, node) {
}
return(object)
}
## Run replicated K-fold CV with random splits
##
## @param std_data standardized full data set (standardized by `.standardize_data`)
## @param cv_k number of folds per CV split
## @param cv_rep number of CV replications.
## @param cv_est_fun function taking the standardized training set and the indices of the left-out observations and
## returns a list of estimates. The function always needs to return the same number of estimates!
## @param metric function taking a vector of prediction errors and returning the scale of the prediction error.
#' @importFrom Matrix drop
.run_replicated_cv <- function (std_data, cv_k, cv_rep, cv_est_fun, metric) {
est_fun <- match.fun(cv_est_fun)
prediction_metrics <- lapply(integer(cv_rep), function (.) {
test_segments <- split(seq_along(y), sample(rep_len(seq_len(cv_k), length(y))))
prediction_errors <- lapply(test_segments, function (test_ind) {
train_x <- x[-test_ind, , drop = FALSE]
train_y <- y[-test_ind]
test_x <- x[test_ind, , drop = FALSE]
test_y <- y[test_ind]
train_std <- std_data$standardize(train_x, train_y)
cv_ests <- est_fun(train_std, test_ind)
matrix(unlist(lapply(cv_ests, function (est) {
unstd_est <- std_data$unstandardize_coef(est)
test_y - drop(test_x %*% unstd_est$beta) - unstd_est$intercept
}), use.names = FALSE, recursive = FALSE), ncol = length(cv_ests))
})
prediction_errors_mat <- do.call(rbind, prediction_errors)
apply(prediction_errors_mat, 2, metric)
})
matrix(unlist(prediction_metrics, recursive = FALSE, use.names = FALSE), ncol = cv_rep)
}
## Standardize data
##
## @param x predictor matrix. Can also be a list with components `x` and `y`, in which case `y` is ignored.
## @param y response vector.
## @param standardize standardize or not.
## @param robust use robust standardization.
## @param location_rho,location_cc,mscale_opts,... further arguments to `mlocscale()` and `mloc()`
## @return a list with the following entries:
#' @importFrom Matrix drop
#' @importFrom methods is
.standardize_data <- function (x, y, standardize = TRUE, robust = TRUE, location_rho = 'bisquare', location_cc = 4.5,
mscale_opts = mscale_algorithm_options(), ...) {
if (is.list(x) && !is.null(x$x) && !is.null(x$y)) {
y <- x$y
x <- x$x
}
ret_list <- list(scale_x = 1, mux = 0, muy = 0, x = x, y = y)
## standardize data
if (isTRUE(standardize)) {
if (!isTRUE(robust)) {
ret_list$scale_x <- apply(x, 2, sd)
ret_list$mux <- colMeans(x)
ret_list$muy <- mean(y)
} else {
locscale <- apply(x, 2, function (xj) {
mlocscale(xj, location_rho = location_rho, location_cc = location_cc, ...)
})
ret_list$mux <- locscale[1L, ]
ret_list$scale_x <- locscale[2L, ]
ret_list$muy <- mloc(y, rho = location_rho, cc = location_cc)
}
if (!isTRUE(all(ret_list$scale_x > 0))) {
stop("Standardization failed. One or more variables in x have a scale of 0.")
}
ret_list$x <- scale(x, center = ret_list$mux, scale = ret_list$scale_x)
ret_list$y <- y - ret_list$muy
}
ret_list$standardize <- function(x, y) {
if (is.list(x) && !is.null(x$x) && !is.null(x$y)) {
y <- x$y
x <- x$x
}
.standardize_data(x, y, standardize = standardize, robust = robust, location_rho = location_rho,
location_cc = location_cc, mscale_opts = mscale_opts, ... = ...)
}
ret_list$standardize_coefs <- function(coef_obj) {
if (!isTRUE(standardize) || is.null(coef_obj)) {
return(coef_obj)
}
if (is.null(coef_obj$intercept)) {
coef_obj$intercept <- 0
}
coef_obj$intercept <- coef_obj$intercept - ret_list$muy + drop(ret_list$mux %*% coef_obj$beta)
coef_obj$beta <- coef_obj$beta * ret_list$scale_x
return(coef_obj)
}
ret_list$unstandardize_coefs <- function(coef_obj) {
if (!isTRUE(standardize) || is.null(coef_obj)) {
return(coef_obj)
}
if (is.null(coef_obj$intercept)) {
coef_obj$intercept <- 0
}
if (is(coef_obj$beta, 'dsparseVector')) {
coef_obj$beta@x <- coef_obj$beta@x / ret_list$scale_x[coef_obj$beta@i]
} else {
coef_obj$beta <- coef_obj$beta / ret_list$scale_x
}
coef_obj$intercept <- coef_obj$intercept + ret_list$muy - drop(ret_list$mux %*% coef_obj$beta)
return(coef_obj)
}
return(ret_list)
}
.cv_mape <- function (r) {
median(abs(r))
}
.cv_rmspe <- function (r) {
sqrt(mean(r^2))
}
#' @importFrom methods is
#' @importFrom Matrix sparseVector
.concat_coefs <- function (coef, call) {
if (is(coef$beta, 'dsparseVector')) {
sparseVector(c(coef$intercept, coef$beta@x), i = c(1L, coef$beta@i + 1L), length = coef$beta@length + 1L)
} else {
coefvec <- c(coef$intercept, coef$beta)
varnames <- tryCatch(colnames(eval(call$x)), error = function (e) NULL)
if (!is.null(varnames)) {
names(coefvec) <- c('(Intercept)', varnames)
}
return(coefvec)
}
}
.cv_se_selection <- function (cvm, cvsd, se_fact) {
type <- rep.int(factor('none', levels = c('none', 'min', 'se_fact')), length(cvm))
best <- which.min(cvm)
candidates <- which(cvm <= cvm[[best]] + se_fact * cvsd[[best]])
candidates <- candidates[candidates <= best] # only consider sparser solutions
# "ignore" solutions after which the prediction performance comes back down
best_1se <- if (any(diff(candidates) > 1)) {
min(candidates[seq_len(max(which(diff(candidates) > 1)))])
} else {
min(candidates)
}
type[[best]] <- 'min'
type[[best_1se]] <- 'se_fact'
return(type)
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/elnet.R
\name{coef.elnet_cv_res}
\alias{coef.elnet_cv_res}
\title{Extract Model Coefficients}
\usage{
\method{coef}{elnet_cv_res}(
object,
lambda = c("min", "se"),
se_mult = 2,
exact = deprecated(),
sparse = deprecated(),
...
)
}
\arguments{
\item{object}{LS-EN estimate fitted by \code{\link[=elnet]{elnet()}} to extract coefficients from.}
\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{exact}{deprecated. Always gives a warning if \code{lambda} is not part of the fitted sequence and coefficients
are interpolated.}
\item{sparse}{deprecated. Always uses the \code{sparse} argument in the call to \code{\link[=elnet]{elnet()}}.}
\item{...}{currently not used.}
}
\value{
either a numeric vector or a sparse vector of type \link[Matrix:sparseVector-class]{dsparseVector}
of size \eqn{p + 1}, depending on the \code{sparse} argument in the call to \code{\link[=elnet]{elnet()}}.
}
\description{
Extract coefficients from the LS-EN estimate with hyper-parameters chosen by cross-validation.
}
\details{
If \code{lambda = "se"} and \code{object} has fitted estimates for every penalization level in the sequence, extract the
coefficients of the most parsimonious model with prediction performance statistically indistinguishable from the best
model. This is determined to be the model with prediction performance within \code{se_mult * cv_se} from the best model.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/elnet.R
\name{coef.elnet_res}
\alias{coef.elnet_res}
\title{Extract Model Coefficients}
\usage{
\method{coef}{elnet_res}(object, lambda, exact = deprecated(), sparse = deprecated(), ...)
}
\arguments{
\item{object}{LS-EN estimate fitted by \code{\link[=elnet]{elnet()}} to extract coefficients from.}
\item{lambda}{a single value of the penalty parameter.}
\item{exact}{deprecated. Always gives a warning if \code{lambda} is not part of the fitted sequence and coefficients
are interpolated.}
\item{sparse}{deprecated. Always uses the \code{sparse} argument in the call to \code{\link[=elnet]{elnet()}}.}
\item{...}{currently not used.}
}
\value{
either a numeric vector or a sparse vector of type \link[Matrix:sparseVector-class]{dsparseVector}
of size \eqn{p + 1}, depending on the \code{sparse} argument in the call to \code{\link[=elnet]{elnet()}}.
}
\description{
Extract Model Coefficients
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/elnet.R
\name{.elnet_internal}
\alias{.elnet_internal}
\title{Perform some final input adjustments and call the internal C++ code.}
\usage{
.elnet_internal(
x,
y,
alpha,
lambdas,
penalty_loadings = NULL,
weights = NULL,
include_intercept = TRUE,
optional_args
)
}
\description{
Perform some final input adjustments and call the internal C++ code.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/pense_regression.R
\name{.pense_lambda_grid}
\alias{.pense_lambda_grid}
\title{Generate a log-spaced grid of decreasing lambda values}
\usage{
.pense_lambda_grid(
x,
y,
alpha,
nlambda,
lambda_min_ratio,
pense_options,
penalty_loadings = NULL
)
}
\description{
Generate a log-spaced grid of decreasing lambda values
}
......@@ -2,34 +2,51 @@
% Please edit documentation in R/elnet.R
\name{elnet}
\alias{elnet}
\title{Compute the Elastic Net Regularization Path}
\title{Compute the Least Squares Elastic Net Regularization Path}
\usage{
elnet(
x,
y,
alpha,
lambdas,
nlambda = 100,
lambda_min_ratio,
lambda,
penalty_loadings,
weights,
include_intercept = TRUE,
intercept = TRUE,
en_algorithm_opts,
sparse = FALSE,
eps = 1e-06
eps = 1e-06,
standardize = TRUE,
correction = deprecated(),
xtest = deprecated(),
options = deprecated()
)
}
\arguments{
\item{x}{\code{n} by \code{p} matrix of numeric predictors.}
\item{y}{vector of response values of length \code{n}.}
\item{alpha}{value for the alpha parameter, the balance between L1 and L2
penalization.}
\item{lambdas}{a vector of positive values for the lambda parameter.}
\item{nlambda}{number of penalization levels.}
\item{lambda_min_ratio}{Smallest value of the penalization level as a fraction of the largest level (i.e., the
smallest value for which all coefficients are zero). The default depends on the sample
size relative to the number of variables and \code{alpha}. If more observations than variables
are available, the default is \code{1e-3 * alpha}, otherwise \code{1e-2 * alpha}.}
\item{lambda}{optional user-supplied sequence of penalization levels. If given and not \code{NULL}, \code{nlambda} and
\code{lambda_min_ratio} are ignored.}
\item{penalty_loadings}{a vector of positive penalty loadings
(a.k.a. weights) for different penalization of each
\item{penalty_loadings}{a vector of positive penalty loadings (a.k.a. weights) for different penalization of each
coefficient.}
\item{weights}{a vector of positive observation weights.}
\item{include_intercept}{include an intercept in the model.}
\item{intercept}{include an intercept in the model.}
\item{en_algorithm_opts}{options for the EN algorithm. See \link{en_algorithm_options}
for details.}
......@@ -37,9 +54,32 @@ for details.}
\item{sparse}{use sparse coefficient vectors.}
\item{eps}{numerical tolerance.}
\item{standardize}{standardize variables to have unit variance. Coefficients are always returned in original scale.}
\item{correction}{deprecated. Instead, extract coefficients using \code{\link[=coef.elnet_est]{coef.elnet_est()}}.}
\item{xtest}{deprecated. Instead, predict values with \code{\link[=predict.elnet_est]{predict.elnet_est()}}.}
\item{options}{deprecated. Use \code{en_algorithm_opts} instead.}
}
\value{
\describe{
\item{\code{call}}{the original call.}
\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{status}}{if \verb{> 0} the algorithm had an issue with computing the estimate}
\item{\code{status_message}}{optional status message from the algorithm}
}
}
}
}
\description{
Compute the EN estimator for linear regression with optional observation
Compute least squares EN estimates for linear regression with optional observation
weights and penalty loadings.
}
\details{
......@@ -47,11 +87,13 @@ 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| }
(1/2n) \sum_i w_i (y_i - \mu - x_i' \beta)^2 +
\lambda \sum_j 0.5 (1 - \alpha) \beta_j^2 + \alpha l_j |\beta_j| }
with observation weights \eqn{w_i} and penalty loadings \eqn{l_i}.
with observation weights \eqn{w_i} and penalty loadings \eqn{l_j}.
}
\seealso{
\link{pense} for an S-estimate of regression with elastic net penalty.
\code{\link[=elnet_cv]{elnet_cv()}} for cross-validating prediction performance of the estimates.
\code{\link[=pense]{pense()}} for an S-estimate of regression with elastic net penalty.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/elnet.R
\name{elnet_cv}
\alias{elnet_cv}
\title{Cross-validation for Least-Squares Elastic Net Estimates}
\usage{
elnet_cv(
x,
y,
lambda,
cv_k = 100,
cv_repl = 1,
cv_metric = c("rmspe", "tau_size", "mape"),
fit_all = TRUE,
...
)
}
\arguments{
\item{x}{\code{n} by \code{p} matrix of numeric predictors, as for \code{\link[=elnet]{elnet()}}.}
\item{y}{vector of response values of length \code{n}, as for \code{\link[=elnet]{elnet()}}.}
\item{lambda}{optional user-supplied sequence of penalization levels. For details see \code{\link[=elnet]{elnet()}}.}
\item{cv_k}{number of folds per cross-validation.}
\item{cv_repl}{number of cross-validation replications.}
\item{cv_metric}{either a string specifying the performance metric to use, or a function to evaluate prediction
errors in a single CV replication. If a function, it is called with a single numeric vector of
prediction errors and must return a scalar number.}
\item{fit_all}{If \code{TRUE}, fit the EN estimate for all penalization levels. Otherwise, only at penalization
level with smallest average CV performance.}
\item{...}{further arguments to \code{\link[=elnet]{elnet()}}.}
}
\value{
a list with components:
\describe{
\item{\code{cv_performance}}{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{cv_performance}.}
\item{\code{call}}{the original call.}
\item{\code{estimates}}{the estimates fitted on the full data. Same format as returned by \code{\link[=elnet]{elnet()}}.}
}
}
\description{
Perform (repeated) K-fold cross-validation for \code{\link[=elnet]{elnet()}}.
}
\details{
The built-in CV metrics are
\describe{
\item{\code{"rmspe"}}{Root mean squared prediction error (default).}
\item{\code{"mape"}}{Median absolute prediction error.}
\item{\code{"tau_size"}}{\eqn{\tau}-size of the prediction error, computed by \code{\link[=tau_size]{tau_size()}}.}
}
}
\seealso{
\code{\link[=pense_cv]{pense_cv()}} for cross-validation for S-estimates of regression with elastic net penalty.
}
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