...
 
Commits (4)
......@@ -34,3 +34,5 @@ Suggests:
testthat (>= 2.1.0)
License: MIT + file LICENSE
NeedsCompilation: yes
RoxygenNote: 7.1.0
Roxygen: list(markdown = TRUE, load = "source")
# Generated by roxygen2: do not edit by hand
S3method(coef,post_lasso)
S3method(coef,thrive)
export(post_lasso)
export(post_lasso_cv)
export(thrive)
importFrom(Rcpp,evalCpp)
useDynLib(thrive, .registration = TRUE)
#' Extract Model Coefficients
#'
#' @param object a post-LASSO estimate to extract coefficients from.
#' @param s the number of variables in the model (for `what = "plasso"`) or the value of the
#' penalty parameter (for `what = "lasso"`).
#' @param sparse return sparse vector or dense (base R) vectors.
#' @param ... currently not used.
#' @return a list of coefficients. Each coefficient is a list with two components: `intercept` and `beta`.
#' If `sparse = FALSE`, `beta` is a numeric vector of size \eqn{p}, otherwise it is a sparse vector.
#'
#' @export
coef.post_lasso <- function (object, s, what = c("plasso", "lasso"), sparse = TRUE, ...) {
what <- match.arg(what)
res <- switch(what, plasso = .plasso_coefs(object, s), lasso = .lasso_coefs(object, s))
if (!isTRUE(sparse)) {
return(as.numeric(res))
}
return(res)
}
#' Extract Model Coefficients
#'
#' @param object a thrive estimate to extract coefficients from.
#' @param nnz the number of exposures with non-zero coefficient.
#' @param sparse return a sparse vector or a dense (base R) numeric vector
#' @return a list of coefficients. Each coefficient is a list with two components: `intercept` and `beta`.
#' If `sparse = FALSE`, `beta` is a numeric vector of size \eqn{p}, otherwise it is a sparse vector.
#'
#' @export
coef.thrive <- function (object, nnz, sparse = TRUE, ...) {
nvars <- sapply(object$estimates, function (x) length(x$beta@i) )
lapply(nnz, function (nk) {
ret_ind <- which(nk == nvars)
if (length(ret_ind) == 0L) {
warning(sprintf("requested number of variables (%d) not available", nk))
return(list())
}
if (isFALSE(sparse)) {
object$estimates[[ret_ind]]$beta <- as.numeric(object$estimates[[ret_ind]]$beta)
}
return(object$estimates[[ret_ind]])
})
}
.plasso_coefs <- function (object, nkeep) {
nvars <- sapply(object$plasso, function (x) length(x$beta@i) )
unlist(lapply(nkeep, function (nk) {
ret_ind <- match(nk, nvars)
if (length(ret_ind) == 0L || anyNA(ret_ind)) {
warning(sprintf("requested number of variables (%d) not available", nk))
return(list())
}
return(object$plasso[ret_ind])
}), recursive = FALSE, use.names = FALSE)
}
.lasso_coefs <- function(object, lambda) {
lambda_seq <- unlist(lapply(object$lasso, `[[`, 'lambda'))
if (length(lambda_seq) == 0L) {
stop("No LASSO solutions computed.")
}
lapply(lambda, function (l) {
if (l >= max(lambda_seq)) {
return(object$lasso[[length(object$lasso)]])
} else if (l <= min(lambda_seq)) {
return(object$lasso[[1L]])
}
upper_lambda_ind <- min(which(l <= lambda_seq))
lower_lambda_ind <- max(which(l >= lambda_seq))
if (upper_lambda_ind == lower_lambda_ind) {
return(object$lasso[[upper_lambda_ind]])
}
interp_beta <- ((l - lambda_seq[[lower_lambda_ind]]) * object$lasso[[upper_lambda_ind]]$beta +
(lambda_seq[[upper_lambda_ind]] - l) * object$lasso[[lower_lambda_ind]]$beta) /
(lambda_seq[[upper_lambda_ind]] - lambda_seq[[lower_lambda_ind]])
interp_int <- ((l - lambda_seq[[lower_lambda_ind]]) * object$lasso[[upper_lambda_ind]]$intercept +
(lambda_seq[[upper_lambda_ind]] - l) * object$lasso[[lower_lambda_ind]]$intercept) /
(lambda_seq[[upper_lambda_ind]] - lambda_seq[[lower_lambda_ind]])
return(list(intercept = interp_int, beta = interp_beta, lambda = l))
})
}
#' Compute (Post-)LASSO Estimates of Linear Regression
#'
#' Estimate the LASSO and Post-LASSO (PLASSO) regression coefficients.
#' Post-LASSO estimates are computed for all models of size in `post_keep`, whereas LASSO estimates are computed
#' for all penalization levels in `lambda`.
#'
#' Post-LASSO estimates are ordinary least square estimates of the coefficients for variables selected by LASSO.
#'
#' The minimization problem is formulated as
#' \deqn{\frac{1}{2 N} RSS + \lambda \left(
#' \frac{(1 - \alpha)} {2} \| \beta \|_2^2 + \alpha \| \beta \|_1
#' \right)}{
#' (1/2N) RSS + \lambda * (
#' (1 - \alpha) / 2 * L2(\beta)^2 + \alpha * L1(\beta)
#' )}
#'
#' @param x data matrix with variables.
#' @param y response vector.
#' @param post_keep sizes of models for which post-LASSO estimates are to be computed. If missing, post-LASSO estimates
#' are computed for all possible models, i.e., models of size 1 ... `ncol(x)`.
#' If `NULL`, no post-LASSO estimates are computed.
#' @param nlambda number of penalization levels for which LASSO estimates are computed.
#' @param 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. If more observations than variables
#' are available, the default is `1e-3`, otherwise `1e-2`.
#' @param lambda optional user-supplied sequence of penalization levels. If given and not `NULL`, `nlambda` and
#' `lambda_min_ratio` are ignored.
#' @param intercept should an intercept be estimated?
#' @return
#' \item{lambda}{vector of lambda values.}
#' \item{lasso}{list of LASSO regression coefficients.}
#' \item{plasso}{list of PLASSO regression coefficients.}
#'
#' @importFrom Rcpp evalCpp
#' @useDynLib thrive, .registration = TRUE
#' @export
post_lasso <- function(x, y, post_keep, nlambda = 0, lambda, lambda_min_ratio, intercept = TRUE) {
cl <- match.call()
args_call <- cl
args_call[[1L]] <- quote(thrive:::.post_lasso_args)
args <- eval.parent(args_call)
plassores <- .Call(C_fastplasso, args$x, args$y, args$intercept, args$post_keep, args$lambda)
plassores$has_intercept <- args$intercept
plassores$call <- cl
class(plassores) <- 'post_lasso'
return(plassores)
}
#' @rdname post_lasso
#' @description [post_lasso_cv()] automatically cross-validates (post)-LASSO estimates.
#'
#' @param cv_k number of cross-validation folds.
#' @param cv_repl number of replications of cross-validation.
#' @param 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.
#' @export
post_lasso_cv <- function (x, y, post_keep, cv_k, nlambda = 0, lambda, lambda_min_ratio, intercept = TRUE, cv_repl = 1,
cv_metric = c('rmspe', 'mape')) {
cl <- match.call()
args_call <- cl
args_call[[1L]] <- quote(thrive:::.post_lasso_args)
args <- eval.parent(args_call)
x_dim <- dim(args$x)
cv_k <- as.integer(cv_k[[1L]])
cv_repl <- as.integer(cv_repl[[1L]])
if (!isTRUE(cv_k > 1L)) {
stop("`cv_k` must be greater than 1.")
}
if (!isTRUE(cv_repl > 0L)) {
stop("`cv_repl` must be greater than 0.")
}
cv_metric <- if (is.character(cv_metric)) {
cv_measure_str <- match.arg(cv_metric)
switch(cv_measure_str, mape = .cv_mape, rmspe = .cv_rmspe, tau_size = tau_size)
} else {
cv_measure_str <- 'user_fun'
match.fun(cv_metric)
}
# Fit to the full data
fit_full <- .Call(C_fastplasso, args$x, args$y, args$intercept, args$post_keep, args$lambda)
nact_full <- vapply(fit_full$plasso, FUN.VALUE = integer(1L), function (coef) length(coef$beta@i))
nact_full_present <- unique(nact_full)
args$post_keep <- rev(nact_full_present)
cv_base_seed <- sample.int(max(cv_repl, .Machine$integer.max - cv_repl), 1L)
cv_perf <- lapply(seq_len(cv_repl), function (repl) {
# Create CV segments
set.seed(cv_base_seed + repl)
cv_segments <- split(seq_len(x_dim[[1L]]), sample(rep_len(seq_len(cv_k), x_dim[[1L]])))
pred_errs <- lapply(cv_segments, function(cv_ind) {
x_train <- args$x[-cv_ind, , drop = FALSE]
y_train <- args$y[-cv_ind]
x_test <- args$x[cv_ind, , drop = FALSE]
y_test <- args$y[cv_ind]
train_fit <- .Call(C_fastplasso, x_train, y_train, args$intercept, args$post_keep, args$lambda)
prederr_lasso <- if (length(train_fit$lasso) > 0L) {
vapply(train_fit$lasso, FUN.VALUE = numeric(length(y_test)), function (coef) {
y_test - coef$intercept - drop(x_test %*% coef$beta)
})
} else {
matrix(nrow = 0L, ncol = 0L)
}
prederr_plasso <- if (length(train_fit$plasso) > 0L) {
plasso_nact <- vapply(train_fit$plasso, FUN.VALUE = integer(1L), function (coef) {
length(coef$beta@i)
})
plasso_prederr <- lapply(train_fit$plasso, function (coef) {
y_test - coef$intercept - as.numeric(x_test %*% coef$beta)
})
if (any(duplicated(plasso_nact))) {
plasso_prederr <- vapply(split(plasso_prederr, plasso_nact), FUN.VALUE = numeric(length(y_test)),
function (pes) {
cvm <- vapply(pes, FUN.VALUE = numeric(1L), cv_metric)
pes[[which.min(cvm)]]
}, USE.NAMES = FALSE)
plasso_nact <- unique(plasso_nact)
} else {
plasso_prederr <- matrix(unlist(plasso_prederr), nrow = length(y_test))
}
plasso_prederr[ , match(nact_full_present, plasso_nact)]
} else {
numeric(0L)
}
return(list(lasso = prederr_lasso, plasso = prederr_plasso))
})
list(lasso = apply(do.call(rbind, lapply(pred_errs, `[[`, 'lasso')), 2, cv_metric),
plasso = apply(do.call(rbind, lapply(pred_errs, `[[`, 'plasso')), 2, cv_metric))
})
fit_full$cv_perf_lasso <- .make_cv_performance_df(cv_perf, 'lasso', cv_repl = cv_repl, lambda = rev(args$lambda))
fit_full$cv_perf_plasso <- .make_cv_performance_df(cv_perf, 'plasso', cv_repl = cv_repl, nnz = rev(args$post_keep))
fit_full$call <- cl
class(fit_full) <- 'post_lasso_cv'
return(fit_full)
}
.make_cv_performance_df <- function (cv_performance, what, cv_repl, ...) {
cv_objs <- matrix(unlist(lapply(cv_performance, `[[`, what)), ncol = cv_repl)
data.frame(..., cvm = rowMeans(cv_objs), cvse = apply(cv_objs, 1, sd))
}
.cv_mape <- function (r) {
median(abs(r))
}
.cv_rmspe <- function (r) {
sqrt(mean(r^2))
}
.post_lasso_args <- function (x, y, post_keep, nlambda = 0, lambda, lambda_min_ratio, intercept = TRUE, ...) {
args <- list()
y <- drop(y)
x_dim <- dim(x)
y_dim <- dim(y)
yl <- length(y)
if (is.null(yl) || (!is.null(y_dim) && length(y_dim) != 1L) || !is.numeric(y)) {
stop("`y` must be a numeric vector.")
}
if (is.null(x_dim) || length(x_dim) != 2L || !is.numeric(x) || x_dim[1L] != yl) {
stop("`x` must be a numeric matrix with the same number of observations as `y`.")
}
if (anyNA(x) || anyNA(y)) {
stop("Missing values in `x` or `y` are not supported.")
}
args$post_keep <- if (missing(post_keep)) {
seq_len(x_dim[[2L]])
} else if (is.null(post_keep)) {
integer(0L)
} else if (is.numeric(post_keep)) {
sort(as.integer(post_keep))
} else {
stop("`post_keep` must be an integer vector.")
}
args$lambda <- if (missing(lambda) || is.null(lambda)) {
if (isTRUE(nlambda > 0L)) {
# Generate lambda sequence automatically.
if (missing(lambda_min_ratio) || is.null(lambda_min_ratio)) {
lambda_min_ratio <- if (x_dim[[1L]] > x_dim[[2L]]) { 1e-3 } else { 1e-2 }
}
max_lambda <- ((x_dim[[1L]] - 1) / x_dim[[1L]]) * max(abs(cov(x, y)))
rev(exp(seq(log(lambda_min_ratio * max_lambda), log(max_lambda), length.out = nlambda)))
} else {
numeric(0L)
}
} else if (is.numeric(lambda)) {
sort(lambda, decreasing = TRUE)
} else {
stop("`lambda` must be a numeric vector.")
}
args$intercept <- !isFALSE(intercept)
args$x <- x
args$y <- y
return(args)
}
#' Compute Thrice-Robust Instrumental Variable Estimates
#'
#' Compute thrice-robust instrumental variable estimates.
#' The exposures `x` are first de-confounded using instruments `z`, and the response `y` is then regressed on
#' the de-confounded exposures.
#' Estimates for de-confounding the exposures are computed via post-LASSO and the reguarlization parameter
#' is selected via cross-validation.
#'
#'
#' @param x matrix of exposures. Should be properly standardized beforehand.
#' @param y vector of response values.
#' @param z matrix of potential instruments for all exposures. Should be properly standardized beforehand.
#' @param model_size_span span of the grid of model sizes to consider in the second stage. See details.
#' @param model_sizes_x model sizes to consider for the exposures. If missing, consider all possible model sizes
#' for each exposure (i.e., models of size 1 ... `ncol(z)`). If an integer vector,
#' the same model sizes are considered for every exposure. If a list the same length as `x`
#' columns, different model sizes are considered for each exposure.
#' @param model_size_y model sizes to consider for the second stage regression, i.e., for regressing the
#' response `y` on the de-confounded matrix of exposures.
#' @param intercept should an intercept be estimated?
#' @param cv_k number of cross-validation folds.
#' @param cv_repl number of replications of cross-validation.
#' @param 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.
#' @return
#' \item{plasso}{list of 2nd-stage PLASSO regression coefficients.}
#' \item{x_fitted}{the fitted exposures.}
#' \item{cv_perf_plasso}{cross-validated prediction performance of the estimates based on `x_fitted`.}
#' \item{cv_perf_grid}{cross-validated prediction performance of the estimates based on all the considered fitted
#' exposures.}
#' \item{z_coefs}{coefficients for regressing each exposure on the instruments.}
#' @export
thrive <- function (x, y, z, cv_k, model_size_span = 1, model_sizes_x, model_size_y,
intercept = TRUE, cv_repl = 1, cv_metric = c('rmspe', 'mape')) {
y <- drop(y)
x_dim <- dim(x)
z_dim <- dim(z)
y_dim <- dim(y)
yl <- length(y)
if (is.null(yl) || (!is.null(y_dim) && length(y_dim) != 1L) || !is.numeric(y)) {
stop("`y` must be a numeric vector.")
}
if (is.null(x_dim) || length(x_dim) != 2L || !is.numeric(x) || x_dim[1L] != yl) {
stop("`x` must be a numeric matrix with the same number of observations as `y`.")
}
if (is.null(x_dim) || length(x_dim) != 2L || !is.numeric(x) || x_dim[1L] != yl) {
stop("`z` must be a numeric matrix with the same number of observations as `y`.")
}
if (anyNA(x) || anyNA(y) || anyNA(z)) {
stop("Missing values in `x`, `y`, or `z` are not supported.")
}
cv_metric <- if (is.character(cv_metric)) {
cv_measure_str <- match.arg(cv_metric)
switch(cv_measure_str, mape = .cv_mape, rmspe = .cv_rmspe, tau_size = tau_size)
} else {
cv_measure_str <- 'user_fun'
match.fun(cv_metric)
}
model_size_y <- if (missing(model_size_y) || is.null(model_size_y)) {
seq_len(x_dim[[2L]] + 1L) - 1L
} else if (is.numeric(model_size_y)) {
sort(as.integer(model_size_y))
} else {
stop("`model_size_y` must be an integer vector.")
}
model_sizes_x <- if (missing(model_sizes_x) || is.null(model_sizes_x)) {
rep.int(list(seq_len(z_dim[[2L]] + 1L) - 1L), x_dim[[2L]])
} else if (is.list(model_sizes_x)) {
if (length(model_sizes_x) != x_dim[[2L]]) {
stop("If `model_sizes_x` is a list, it must contain one integer vector per column in `x`.")
}
lapply(model_sizes_x, function (msx) {
if (!is.numeric(msx)) {
stop("If `model_sizes_x` is a list, it must contain one integer vector per column in `x`.")
}
sort(as.integer(msx))
})
} else if (is.numeric(model_sizes_x)) {
rep.int(list(sort(as.integer(model_size_y))), x_dim[[2L]])
} else {
stop("`model_size_y` must be an integer vector or a list of integer vectors.")
}
intercept <- !isFALSE(intercept)
# First stage, regress each exposure on all instruments
stage_1 <- lapply(seq_len(x_dim[[2L]]), function (j) {
pr <- post_lasso_cv(x = z, y = x[, j, drop = TRUE], post_keep = model_sizes_x[[j]], lambda = numeric(0L),
intercept = intercept, cv_k = cv_k, cv_repl = cv_repl, cv_metric = cv_metric)
best_model <- which.min(pr$cv_perf_plasso$cvm)
model_grid <- .hard_threshold(best_model + c(-rev(seq_len(model_size_span)), 0L, seq_len(model_size_span)),
1L, length(pr$cv_perf_plasso$cvm))
coefs <- pr$plasso[model_grid]
fitted <- lapply(coefs, function (coef) {
as.numeric(z %*% coef$beta) + coef$intercept
})
list(coefs = coefs, fitted = fitted)
})
stage_2 <- lapply(seq_len(2L * model_size_span + 1L), function (k) {
x_fitted <- vapply(stage_1, FUN.VALUE = y, function (s1) { s1$fitted[[k]] })
pr <- post_lasso_cv(x = x_fitted, y = y, post_keep = model_size_y, lambda = numeric(0L), intercept = intercept,
cv_k = cv_k, cv_repl = cv_repl, cv_metric = cv_metric)
pr$x_fitted <- x_fitted
pr$cv_perf_plasso$grid_step <- k - model_size_span - 1L
pr$z_coefs <- lapply(stage_1, function (s1) { s1$coefs[[k]] })
return(pr)
})
stage_2_min_cv <- vapply(stage_2, FUN.VALUE = numeric(1L), function (s2) {
min(s2$cv_perf_plasso$cvm, na.rm = TRUE)
})
best_stage_2 <- stage_2[[which.min(stage_2_min_cv)]]
best_stage_2$cv_perf_lasso <- NULL
best_stage_2$lasso <- NULL
best_stage_2$cv_perf_grid <- do.call(rbind, lapply(stage_2, `[[`, 'cv_perf_plasso'))
class(best_stage_2) <- 'thrive'
names(best_stage_2)[[which(names(best_stage_2) == 'plasso')]] <- 'estimates'
return(best_stage_2)
}
.hard_threshold <- function (x, .min, .max) {
pmax(.min, pmin(x, .max))
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/coef-methods.R
\name{coef.post_lasso}
\alias{coef.post_lasso}
\title{Extract Model Coefficients}
\usage{
\method{coef}{post_lasso}(object, s, what = c("plasso", "lasso"), sparse = TRUE, ...)
}
\arguments{
\item{object}{a post-LASSO estimate to extract coefficients from.}
\item{s}{the number of variables in the model (for \code{what = "plasso"}) or the value of the
penalty parameter (for \code{what = "lasso"}).}
\item{sparse}{return sparse vector or dense (base R) vectors.}
\item{...}{currently not used.}
}
\value{
a list of coefficients. Each coefficient is a list with two components: \code{intercept} and \code{beta}.
If \code{sparse = FALSE}, \code{beta} is a numeric vector of size \eqn{p}, otherwise it is a sparse vector.
}
\description{
Extract Model Coefficients
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/coef-methods.R
\name{coef.thrive}
\alias{coef.thrive}
\title{Extract Model Coefficients}
\usage{
\method{coef}{thrive}(object, nnz, sparse = TRUE, ...)
}
\arguments{
\item{object}{a thrive estimate to extract coefficients from.}
\item{nnz}{the number of exposures with non-zero coefficient.}
\item{sparse}{return a sparse vector or a dense (base R) numeric vector}
}
\value{
a list of coefficients. Each coefficient is a list with two components: \code{intercept} and \code{beta}.
If \code{sparse = FALSE}, \code{beta} is a numeric vector of size \eqn{p}, otherwise it is a sparse vector.
}
\description{
Extract Model Coefficients
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/post_lasso.R
\name{post_lasso}
\alias{post_lasso}
\alias{post_lasso_cv}
\title{Compute (Post-)LASSO Estimates of Linear Regression}
\usage{
post_lasso(
x,
y,
post_keep,
nlambda = 0,
lambda,
lambda_min_ratio,
intercept = TRUE
)
post_lasso_cv(
x,
y,
post_keep,
cv_k,
nlambda = 0,
lambda,
lambda_min_ratio,
intercept = TRUE,
cv_repl = 1,
cv_metric = c("rmspe", "mape")
)
}
\arguments{
\item{x}{data matrix with variables.}
\item{y}{response vector.}
\item{post_keep}{sizes of models for which post-LASSO estimates are to be computed. If missing, post-LASSO estimates
are computed for all possible models, i.e., models of size 1 ... \code{ncol(x)}.
If \code{NULL}, no post-LASSO estimates are computed.}
\item{nlambda}{number of penalization levels for which LASSO estimates are computed.}
\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{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. If more observations than variables
are available, the default is \code{1e-3}, otherwise \code{1e-2}.}
\item{intercept}{should an intercept be estimated?}
\item{cv_k}{number of cross-validation folds.}
\item{cv_repl}{number of replications of cross-validation.}
\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.}
}
\value{
\item{lambda}{vector of lambda values.}
\item{lasso}{list of LASSO regression coefficients.}
\item{plasso}{list of PLASSO regression coefficients.}
}
\description{
Estimate the LASSO and Post-LASSO (PLASSO) regression coefficients.
Post-LASSO estimates are computed for all models of size in \code{post_keep}, whereas LASSO estimates are computed
for all penalization levels in \code{lambda}.
\code{\link[=post_lasso_cv]{post_lasso_cv()}} automatically cross-validates (post)-LASSO estimates.
}
\details{
Post-LASSO estimates are ordinary least square estimates of the coefficients for variables selected by LASSO.
The minimization problem is formulated as
\deqn{\frac{1}{2 N} RSS + \lambda \left(
\frac{(1 - \alpha)} {2} \| \beta \|_2^2 + \alpha \| \beta \|_1
\right)}{
(1/2N) RSS + \lambda * (
(1 - \alpha) / 2 * L2(\beta)^2 + \alpha * L1(\beta)
)}
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/thrive.R
\name{thrive}
\alias{thrive}
\title{Compute Thrice-Robust Instrumental Variable Estimates}
\usage{
thrive(
x,
y,
z,
cv_k,
model_size_span = 1,
model_sizes_x,
model_size_y,
intercept = TRUE,
cv_repl = 1,
cv_metric = c("rmspe", "mape")
)
}
\arguments{
\item{x}{matrix of exposures. Should be properly standardized beforehand.}
\item{y}{vector of response values.}
\item{z}{matrix of potential instruments for all exposures. Should be properly standardized beforehand.}
\item{cv_k}{number of cross-validation folds.}
\item{model_size_span}{span of the grid of model sizes to consider in the second stage. See details.}
\item{model_sizes_x}{model sizes to consider for the exposures. If missing, consider all possible model sizes
for each exposure (i.e., models of size 1 ... \code{ncol(z)}). If an integer vector,
the same model sizes are considered for every exposure. If a list the same length as \code{x}
columns, different model sizes are considered for each exposure.}
\item{model_size_y}{model sizes to consider for the second stage regression, i.e., for regressing the
response \code{y} on the de-confounded matrix of exposures.}
\item{intercept}{should an intercept be estimated?}
\item{cv_repl}{number of replications of cross-validation.}
\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.}
}
\value{
\item{plasso}{list of 2nd-stage PLASSO regression coefficients.}
\item{x_fitted}{the fitted exposures.}
\item{cv_perf_plasso}{cross-validated prediction performance of the estimates based on \code{x_fitted}.}
\item{cv_perf_grid}{cross-validated prediction performance of the estimates based on all the considered fitted
exposures.}
\item{z_coefs}{coefficients for regressing each exposure on the instruments.}
}
\description{
Compute thrice-robust instrumental variable estimates.
The exposures \code{x} are first de-confounded using instruments \code{z}, and the response \code{y} is then regressed on
the de-confounded exposures.
Estimates for de-confounding the exposures are computed via post-LASSO and the reguarlization parameter
is selected via cross-validation.
}
## Use the R_HOME indirection to support installations of multiple R version
CXX_STD = CXX11
PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
# PKG_CXXFLAGS= -D__STDC_LIMIT_MACROS
## Debug!
PKG_OPTFLAGS = -O0 -g
PKG_CXXFLAGS= -pipe -Wall -Wextra -pedantic -Wno-unused-parameter -DDEBUG -D__STDC_LIMIT_MACROS
CXX_STD = CXX11
PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
## Debug!
PKG_OPTFLAGS = -O0 -g
PKG_CXXFLAGS= -pipe -Wall -Wextra -pedantic -Wno-unused-parameter -DDEBUG -D__STDC_LIMIT_MACROS @RESTRICT_CXXFLAGS@
CXX_STD = CXX11
## Use the R_HOME indirection to support installations of multiple R version
PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
PKG_CFLAGS= -D__STDC_LIMIT_MACROS
PKG_CXXFLAGS= -D__STDC_LIMIT_MACROS
//
// Rinterface.cpp
// thrive
//
// Created by David Kepplinger on 2016-02-03.
// Copyright © 2016 David Kepplinger. All rights reserved.
//
#include <stdexcept>
#include <memory>
#include <forward_list>
#include <R_ext/Rdynload.h>
#include "./config.h"
#include "./nsoptim.hpp"
using namespace nsoptim;
using namespace arma;
using Coefficients = nsoptim::RegressionCoefficients<arma::sp_vec>;
using CoefsPtrList = std::forward_list< std::unique_ptr<Coefficients> >;
extern "C" void R_init_fastplasso(DllInfo *dll);
namespace {
/**
*
*
* @param x numeric The transpose of the numeric X matrix (size `nvar` x `nobs`)
* @param y numeric The numeric y vector (size `nobs`)
* @param intercept bool Should an intercept be included?
* @param keep integer What steps (nr. of active predictors) to keep. If empty (length=0), keep all steps.
* @param lambda numeric What steps (in terms of the regularizatin parameter lambda) to keep,
* additionally to the ones kept for `keep`. If empty, no additional
* steps are retained.
*
* @return List Returns a list with two elements:
* lasso: list of coefficients with lasso estimates
* plasso: list of coefficients with plasso estimates
*/
SEXP C_fastplasso(SEXP x, SEXP y, SEXP intercept, SEXP keep, SEXP lambda);
struct MatrixDimensions {
int n_rows;
int n_cols;
};
MatrixDimensions GetMatrixDimensions(SEXP matrix) noexcept;
std::unique_ptr<const PredictorResponseData> MakePredictorResponseData(SEXP x, SEXP y);
std::unique_ptr<const vec> MakeVectorView(SEXP numeric_vector) noexcept;
std::unique_ptr<const uvec> MakeIntegerVectorView(SEXP integer_vector) noexcept;
template<typename InputIterator>
SEXP WrapList(const CoefsPtrList& list, const std::string& other_slot, InputIterator other) noexcept;
SEXP WrapList(const CoefsPtrList& list) noexcept;
template <typename T>
SEXP WrapSparseVector(const arma::SpCol<T>& svec);
const R_CallMethodDef exportedCallMethods[] = {
{"C_fastplasso", (DL_FUNC) &C_fastplasso, 5},
{NULL, NULL, 0}
};
} // namespace
extern "C" void R_init_thrive(DllInfo *dll) {
R_registerRoutines(dll, NULL, exportedCallMethods, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
R_forceSymbols(dll, TRUE);
}
namespace {
SEXP C_fastplasso(SEXP x, SEXP y, SEXP intercept, SEXP Rkeep, SEXP Rlambda) {
LsRegressionLoss loss(MakePredictorResponseData(x, y), Rcpp::as<bool>(intercept));
auto nkeep_vec = MakeIntegerVectorView(Rkeep);
auto lambda_vec = MakeVectorView(Rlambda);
Rcpp::List ret_list;
BEGIN_RCPP
AugmentedLarsOptimizer optimizer(loss);
CoefsPtrList keep_coefs;
CoefsPtrList lambda_coefs;
std::forward_list<double> lambda_coefs_lambda_seq;
const auto nkeep_end = nkeep_vec->end();
const auto lambda_end = lambda_vec->end();
auto nkeep_it = nkeep_vec->begin();
auto lambda_it = lambda_vec->begin();
const arma::uword max_nkeep = std::max(loss.data().n_pred(), loss.data().n_obs()) + 1;
while (nkeep_it != nkeep_end || lambda_it != lambda_end) {
const arma::uword next_nkeep = (nkeep_it != nkeep_end) ? *nkeep_it : max_nkeep;
const double next_lambda = (lambda_it != lambda_end) ? *lambda_it : 0;
auto optim = optimizer.Optimize(next_lambda, next_nkeep);
if (optim.interpolated_coefs) {
// the lambda threshold is reached
lambda_coefs.push_front(std::move(optim.interpolated_coefs));
lambda_coefs_lambda_seq.push_front(optim.lambda);
if (lambda_it != lambda_end) {
++lambda_it;
}
}
if (optim.coefs) {
// the nkeep threshold is reached
keep_coefs.push_front(std::move(optim.coefs));
}
if (optim.nnz >= next_nkeep && nkeep_it != nkeep_end) {
++nkeep_it;
}
}
ret_list = Rcpp::List::create(Rcpp::Named("lasso") = WrapList(lambda_coefs, "lambda",
lambda_coefs_lambda_seq.begin()),
Rcpp::Named("plasso") = WrapList(keep_coefs));
VOID_END_RCPP
return Rcpp::wrap(ret_list);
}
MatrixDimensions GetMatrixDimensions(SEXP matrix) noexcept {
SEXP r_dims;
PROTECT(r_dims = Rf_getAttrib(matrix, R_DimSymbol));
int const * const dims = INTEGER(r_dims);
MatrixDimensions mat_dims {dims[0], dims[1]};
UNPROTECT(1);
return mat_dims;
}
std::unique_ptr<const PredictorResponseData> MakePredictorResponseData(SEXP x, SEXP y) {
const int nobs_y = Rf_length(y);
const MatrixDimensions x_dims = GetMatrixDimensions(x);
// Check that x and y have the same nr. of observations
if (nobs_y != x_dims.n_rows) {
throw std::invalid_argument("y and x have differing number of observations");
}
// Check that x and y have the required data type
if (TYPEOF(x) != REALSXP || TYPEOF(y) != REALSXP) {
throw std::invalid_argument("y and x must be numeric");
}
return std::unique_ptr<const PredictorResponseData>(
new PredictorResponseData(mat(REAL(x), x_dims.n_rows, x_dims.n_cols, false, true),
vec(REAL(y), nobs_y, false, true)));
}
std::unique_ptr<const vec> MakeVectorView(SEXP numeric_vector) noexcept {
if (TYPEOF(numeric_vector) != REALSXP) {
return std::unique_ptr<const vec>(new vec());
}
return std::unique_ptr<const vec>(new vec(REAL(numeric_vector), Rf_length(numeric_vector), false, true));
}
std::unique_ptr<const uvec> MakeIntegerVectorView(SEXP int_vector) noexcept {
if (TYPEOF(int_vector) != INTSXP) {
return std::unique_ptr<const uvec>(new uvec());
}
return std::unique_ptr<const uvec>(new uvec(reinterpret_cast<arma::uword*>(INTEGER(int_vector)),
Rf_length(int_vector), false, true));
}
template<typename InputIterator>
SEXP WrapList(const CoefsPtrList& list, const std::string& other_slot, InputIterator other) noexcept {
Rcpp::List r_list;
for (auto&& element : list) {
if (element) {
auto coef_list = Rcpp::List::create(Rcpp::Named("intercept") = element->intercept,
Rcpp::Named("beta") = WrapSparseVector(element->beta));
coef_list.push_back(*other++, other_slot);
r_list.push_back(coef_list);
} else {
r_list.push_back(R_NilValue);
++other;
}
}
return r_list;
}
SEXP WrapList(const CoefsPtrList& list) noexcept {
Rcpp::List r_list;
for (auto&& element : list) {
if (element) {
auto coef_list = Rcpp::List::create(Rcpp::Named("intercept") = element->intercept,
Rcpp::Named("beta") = WrapSparseVector(element->beta));
r_list.push_back(coef_list);
} else {
r_list.push_back(R_NilValue);
}
}
return r_list;
}
// SEXP WrapList(const std::forward_list<double>& list) noexcept {
// Rcpp::List r_list;
// for (auto&& element : list) {
// r_list.push_back(Rcpp::wrap(element));
// }
// return r_list;
// }
//! Wrap a sparse vector into a Matrix::sparseVector object
template <typename T>
SEXP WrapSparseVector(const arma::SpCol<T>& svec) {
const int RTYPE = Rcpp::traits::r_sexptype_traits<T>::rtype;
// important: update internal state of SpMat object
svec.sync();
auto length = Rcpp::IntegerVector::create(svec.n_elem);
// copy the data into R objects
const Rcpp::Vector<RTYPE> values(svec.values, svec.values + svec.n_nonzero);
Rcpp::IntegerVector rowind(svec.row_indices, svec.row_indices + svec.n_nonzero);
// the sparseVector uses 1-based row indices.
for (arma::uword i = 0; i < svec.n_nonzero; ++i) {
rowind[i] += 1;
}
Rcpp::S4 r_sparse_vector("dsparseVector");
r_sparse_vector.slot("length") = length;
r_sparse_vector.slot("i") = rowind;
r_sparse_vector.slot("x") = values;
return r_sparse_vector;
}
} // namespace
//
// config.hpp
// thrive
//
// Created by David Kepplinger on 2020-06-18.
// Copyright © 2020 David Kepplinger. All rights reserved.
//
#ifndef THRIVE_CONFIG_H
#define THRIVE_CONFIG_H
#ifndef DEBUG
# ifndef ARMA_NO_DEBUG
# define ARMA_NO_DEBUG
# endif
#endif
#define NUMERIC_EPS 1e-16
#endif // THRIVE_CONFIG_H
//
// nsoptim.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_HPP_
#define NSOPTIM_HPP_
#if defined(Rcpp_hpp)
#error "The file 'Rcpp.h' should not be included. Please correct to include only 'nsoptim.hpp'."
#endif
#if defined(RcppArmadillo__RcppArmadillo__h)
#error "The file 'RcppArmadillo.h' should not be included. Please correct to include only 'nsoptim.hpp'."
#endif
#include "nsoptim/config.hpp"
#include "nsoptim/armadillo.hpp"
#include "nsoptim/container.hpp"
#include "nsoptim/objective.hpp"
#include "nsoptim/optimizer.hpp"
#endif // NSOPTIM_HPP_
//
// armadillo.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_ARMADILLO_HPP_
#define NSOPTIM_ARMADILLO_HPP_
#define ARMA_USE_CXX11 1
#define ARMA_DONT_USE_OPENMP 1
#ifdef __clang__
# pragma clang diagnostic push
# pragma clang diagnostic ignored "-Weverything"
#endif
#include <RcppArmadillo.h>
#ifdef __clang__
# pragma clang diagnostic pop
#endif
#endif // NSOPTIM_ARMADILLO_HPP_
//
// armadillo.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_CONFIG_HPP_
#define NSOPTIM_CONFIG_HPP_
#define NSOPTIM_METRICS_LEVEL 1
#ifdef NSOPTIM_METRICS_DETAILED
# undef NSOPTIM_METRICS_LEVEL
# define NSOPTIM_METRICS_LEVEL 2
#endif
#ifdef NSOPTIM_METRICS_ENABLED
# undef NSOPTIM_METRICS_LEVEL
# define NSOPTIM_METRICS_LEVEL 1
#endif
#ifdef NSOPTIM_METRICS_DISABLED
# undef NSOPTIM_METRICS_LEVEL
# define NSOPTIM_METRICS_LEVEL 0
#endif
#endif // NSOPTIM_CONFIG_HPP_
//
// container.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_CONTAINER_HPP_
#define NSOPTIM_CONTAINER_HPP_
#include "./container/data.hpp"
#include "./container/regression_coefficients.hpp"
#endif // NSOPTIM_CONTAINER_HPP_
//
// data.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_CONTAINER_DATA_HPP_
#define NSOPTIM_CONTAINER_DATA_HPP_
#include <iostream>
#include "../armadillo.hpp"
namespace nsoptim {
//! Simple structure holding a matrix with predictor data *x* (with `n_obs` rows and `n_pred` colums)
//! and the response vector *y* (with `n_obs` rows)
class PredictorResponseData {
public:
//! Initialize the predictor-response data with empty x and y
PredictorResponseData() noexcept : x_(), y_(), n_obs_(0), n_pred_(0) {}
//! Initialize predictor-response data with the given x and y.
//! @note the given data will be copied!
//!
//! @param other_x predictor matrix to copy.
//! @param other_y response vector to copy.
PredictorResponseData(const arma::mat& other_x, const arma::vec& other_y) noexcept
: x_(other_x), y_(other_y), n_obs_(other_x.n_rows), n_pred_(other_x.n_cols) {}
//! Initialize predictor-response data with the given x and y.
//! @note the given data will be moved to this container!
//!
//! @param other_x predictor matrix to move.
//! @param other_y response vector to move.
PredictorResponseData(arma::mat&& other_x, arma::vec&& other_y) noexcept
: x_(std::move(other_x)), y_(std::move(other_y)), n_obs_(x_.n_rows), n_pred_(x_.n_cols) {}
//! Copy the given predictor-response data, but pointing to the same underlying data!
//!
//! @param other predictor-response data to copy.
PredictorResponseData(const PredictorResponseData& other) = default;
PredictorResponseData& operator=(const PredictorResponseData& other) = default;
//! Get a data set with the observations at the requested indices.
//!
//! @param indices the indicies of the observations to get.
//! @return the subset of the data with the requested observations.
PredictorResponseData Observations(const arma::uvec& indices) const {
return PredictorResponseData(x_.rows(indices), y_.rows(indices));
}
//! Get a data set with the given observation removed.
//!
//! @param index the index of the observation to remove.
//! @return the subset of the data with the observation removed.
PredictorResponseData RemoveObservation(const arma::uword index) const {
return PredictorResponseData(arma::join_vert(x_.head_rows(index), x_.tail_rows(n_obs_ - index - 1)),
arma::join_vert(y_.head(index), y_.tail(n_obs_ - index - 1)));
}
//! Get a data set with the first `n_obs` observations of the data.
//!
//! @param n_obs number of observations to extract.
//! @return the subset of the data with the requested observations.
PredictorResponseData HeadRows(const arma::uword n_obs) const {
return PredictorResponseData(x_.head_rows(n_obs), y_.head_rows(n_obs));
}
//! Get a data set with the last `n_obs` observations of the data.
//!
//! @param n_obs number of observations to extract.
//! @return the subset of the data with the requested rows.
PredictorResponseData TailRows(const arma::uword n_obs) const {
return PredictorResponseData(x_.tail_rows(n_obs), y_.tail_rows(n_obs));
}
//! Get a constant reference to the predictor matrix.
//! Only valid as long as the PredictorResponseData object is in scope.
//!
//! @return constant reference to the predictor matrix
const arma::mat& cx() const noexcept {
return x_;
}
//! Get a constant reference to the response vector.
//! Only valid as long as the PredictorResponseData object is in scope.
//!
//! @return constant reference to the response vector
const arma::vec& cy() const noexcept {
return y_;
}
//! Get non-const references to the data
//! Get a reference to the predictor matrix.
//! Only valid as long as the PredictorResponseData object is in scope.
//!
//! @return reference to the predictor matrix
arma::mat& x() noexcept {
return x_;
}
//! Get a reference to the response vector.
//! Only valid as long as the PredictorResponseData object is in scope.
//!
//! @return reference to the response vector
arma::vec& y() noexcept {
return y_;
}
//! Get the number of observations in this data set.
arma::uword n_obs() const noexcept {
return n_obs_;
}
//! Get the number of observations in this data set.
arma::uword n_pred() const noexcept {
return n_pred_;
}
private:
arma::mat x_;
arma::vec y_;
arma::uword n_obs_; //< The number of observations in the data.
arma::uword n_pred_; //< The number of variables in the data.
};
} // namespace nsoptim
#endif // NSOPTIM_CONTAINER_DATA_HPP_
//
// regression_coefficients.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_CONTAINER_REGRESSION_COEFFICIENTS_HPP_
#define NSOPTIM_CONTAINER_REGRESSION_COEFFICIENTS_HPP_
#include <memory>
#include <type_traits>
#include "../armadillo.hpp"
namespace nsoptim {
//! Regression coefficients in the form of a scalar intercept and a (sparse) vector of slope coefficients.
//! The slope coefficients must be either of type `arma::vec` or `arma::sp_vec`.
template <class T>
class RegressionCoefficients {
static_assert(std::is_same<T, arma::vec>::value || std::is_same<T, arma::sp_vec>::value,
"T must be a (sparse) vector.");
public:
//! The dimensions of regression coefficients is the number if *slope* parameters.
//! The intercept is not included in this number!
using Dimensions = int;
//! The type of the regression slope vector.
using SlopeCoefficient = T;
//! Initialize an empty coefficient vector and an intercept of 0.
RegressionCoefficients() noexcept : intercept(0), beta() {}
//! Initialize the coefficients as 0-vector
//!
//! @param n_pred the number of predictors, i.e., number of slope coefficients.
explicit RegressionCoefficients(const arma::uword n_pred) noexcept
: RegressionCoefficients(n_pred, typename std::is_same<T, arma::sp_vec>::type{}) {}
//! Copy the slope coefficients from the given vector. The intercept will be set to 0.
//!
//! @param beta the value for the slope coefficients.
explicit RegressionCoefficients(const SlopeCoefficient& beta) noexcept : intercept(0), beta(beta) {}
//! Copy the coefficients from the given intercept and slope.
//!
//! @param intercept the value for the intercept.
//! @param beta the value for the slope coefficients.
RegressionCoefficients(const double intercept, const SlopeCoefficient& beta) noexcept
: intercept(intercept), beta(beta) {}
//! Copy constructor
//!
//! @param coef the RegressionCoefficients object to copy.
RegressionCoefficients(const RegressionCoefficients& coef) = default;
//! Move constructor.
//!
//! @param coef the RegressionCoefficients object to move.
RegressionCoefficients(RegressionCoefficients&& coef) = default;
//! Copy assignment
//!
//! @param coef the RegressionCoefficients object to copy.
RegressionCoefficients& operator=(const RegressionCoefficients& other) = default;
//! Move assignment
//!
//! @param coef the RegressionCoefficients object to move.
RegressionCoefficients& operator=(RegressionCoefficients&& other) = default;
//! Reset the coefficients to an empty vector and an intercept of 0.
inline void Reset() noexcept {
intercept = 0;
beta.reset();
}
//! Print the coefficients to the given output stream
//!
//! @param user_stream the stream to print the coefficients to.
inline void Print(std::ostream &user_stream) const {
user_stream << "Intercept: " << intercept << "\n";
beta.t().print(user_stream, "Beta:");
}
//! Print the coefficients to stdout
inline void Print() const {
std::cout << "Intercept: " << intercept << "\n";
beta.t().print("Beta:");
}
double intercept; //< The intercept coefficient
SlopeCoefficient beta; //< The slope coefficients
private:
// Initialize the regression coefficients with all zeros if the slope is a dense vector.
RegressionCoefficients(const arma::uword n_pred, std::false_type) noexcept
: intercept(0), beta(n_pred, arma::fill::zeros) {}
// Initialize the regression coefficients with all zeros if the slope is a sparse vector.
RegressionCoefficients(const arma::uword n_pred, std::true_type) noexcept
: intercept(0), beta(n_pred) {}
};
} // namespace nsoptim
#endif // NSOPTIM_CONTAINER_REGRESSION_COEFFICIENTS_HPP_
//
// objective.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_OBJECTIVE_HPP_
#define NSOPTIM_OBJECTIVE_HPP_
// Loss functions
#include "./objective/ls_regression_loss.hpp"
#endif // NSOPTIM_OBJECTIVE_HPP_
//
// convex.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_OBJECTIVE_CONVEX_HPP_
#define NSOPTIM_OBJECTIVE_CONVEX_HPP_
namespace nsoptim {
//! CRTP helper class for convex functions which returns the object itself as convex surrogate.
template <typename Function>
class ConvexFunction {
public:
using ConvexSurrogateType = Function;
template<typename T>
Function& GetConvexSurrogate(const T&) {
return static_cast<Function&>(*this);
}
};
} // namespace nsoptim
#endif // NSOPTIM_OBJECTIVE_CONVEX_HPP_
//
// loss.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_OBJECTIVE_LOSS_HPP_
#define NSOPTIM_OBJECTIVE_LOSS_HPP_
namespace nsoptim {
//! Boilerplate base class for all loss functions.
//!
//! Loss functions must at least the following methods:
//! `data()` to give read access to the internal data, and
//! `operator(where)` to evaluate the loss at the given coefficients values.
//! `ZeroCoefficients()` to obtain the 0-coefficient value.
//!
//! Loss functions can optionally also implement the following methods:
//! `Difference(a, b)` to evaluate the difference of two coefficient values.
//!
//! Loss functions should be easy and quick to copy and move. The main purpose is not to provide functionality but
//! context.
template<class Data>
class LossFunction {
public:
using DataType = Data;
//! Access the data the loss operates on.
//!
//! @return the data the loss operates on.
//! const Data& data() const;
//! Evaluate the loss function.
//!
//! @param where where to evaluate the loss function.
//! @return the loss evaluated at the given coefficients.
//! double operator()(const Coefficients& where) const;
//! Get the zero coefficients for this loss type.
//!
//! @return zero coefficients.
// Coefficients ZeroCoefficients() const;
//! Get the difference between two sets of coefficients.
//!
//! @param x a set of regression coefficients.
//! @param y the other set of regression coefficients.
//! @return the relative difference between `x` and `y`.
// double Difference(const Coefficients& x, const Coefficients& y) const;
};
} // namespace nsoptim
#endif // NSOPTIM_OBJECTIVE_LOSS_HPP_
//
// ls_regression_loss.hpp
// nsoptim
//
// Created by David Kepplinger on 2018-11-30.
// Copyright © 2018 David Kepplinger. All rights reserved.
//
#ifndef NSOPTIM_OBJECTIVE_LS_REGRESSION_LOSS_HPP_
#define NSOPTIM_OBJECTIVE_LS_REGRESSION_LOSS_HPP_
#include <algorithm>
#include "../armadillo.hpp"
#include "../objective/convex.hpp"
#include "../container/data.hpp"
#include "../container/regression_coefficients.hpp"
#include "../objective/loss.hpp"
namespace nsoptim {
namespace ls_regression_loss {
//! Compute the minimum of the 1- and infinity norm of the weighted matrix :math:`W X`, where :math:`W` is a diagonal
//! matrix with entries `sqrt_weights`.
inline double TwoNormUpper(const arma::mat& x, const arma::vec& sqrt_weights) {
double norm_1 = 0, norm_inf = 0;
// Compute the inf-norm of the weighted matrix
for (arma::uword i = 0; i < x.n_rows; ++i) {
const double tmp = sqrt_weights[i] * arma::norm(x.row(i), 1);
if (tmp > norm_inf) {
norm_inf = tmp;
}
}
// Compute the 1-norm of the weighted matrix
for (arma::uword j = 0; j < x.n_cols; ++j) {
const double tmp = arma::norm(sqrt_weights % x.col(j), 1);
if (tmp > norm_1) {
norm_1 = tmp;
}
}
if (norm_1 < norm_inf) {
return norm_1;
}
return norm_inf;
}
} // namespace ls_regression_loss
//! A regression loss function implementing the un-weighted least-squares loss defined as
//! 1/(2n) * sum_{i = 1}^n (y_i - \hat{\mu} - x_i' . \beta)^2
class LsRegressionLoss : public LossFunction<PredictorResponseData>, public ConvexFunction<LsRegressionLoss> {
public:
template<typename>
using GradientType = RegressionCoefficients<arma::vec>;
//! Tag this loss function as an LS-type loss.
struct is_ls_regression_loss_tag {};
explicit LsRegressionLoss(std::shared_ptr<const PredictorResponseData> data, const bool include_intercept = true)
: include_intercept_(include_intercept), data_(data), pred_norm_(-1) {}
LsRegressionLoss(const LsRegressionLoss& other) = default;
LsRegressionLoss(LsRegressionLoss&& other) = default;
LsRegressionLoss& operator=(const LsRegressionLoss& other) = default;
LsRegressionLoss& operator=(LsRegressionLoss&& other) = default;
~LsRegressionLoss() = default;
//! Get the zero coefficients for this loss type.
//!
//! @return zero coefficients.
template<typename T, typename = typename std::enable_if<
std::is_same<T, RegressionCoefficients<typename T::SlopeCoefficient>>::value, void>::type>
T ZeroCoefficients() const noexcept {
return RegressionCoefficients<typename T::SlopeCoefficient>(data_->n_pred());
}
//! Evaluate the LS loss function.
//!
//! @param where point where to evaluate the loss function.
//! @return loss evaluated at `where`.
template<typename T>
double operator()(const RegressionCoefficients<T>& where) const {
return Evaluate(where);
}
//! Evaluate the LS loss function.
//!
//! @param where point where to evaluate the loss function.
//! @return loss evaluated at `where`.
template<typename T>
double Evaluate(const RegressionCoefficients<T>& where) const {
const arma::vec residuals = data_->cy() - data_->cx() * where.beta;
if (include_intercept_) {
return 0.5 * arma::mean(arma::square(residuals - where.intercept));
}
return 0.5 * arma::mean(arma::square(residuals));
}
//! Get the difference between two sets of regression coefficients.
//!
//! For the LS loss, the relative difference is an approximation to the 2-norm of the matrix-vector product
//! ||X . (beta1 - beta2) + ( mu1 - mu2 )||_2 < |mu1 - mu2| sqrt(n) + ||X|| ||beta1 - beta2||_2
//!
//! @param x a set of regression coefficients.
//! @param y the other set of regression coefficients.
//! @return the relative difference between `x` and `y`.
template<typename T>
double Difference(const RegressionCoefficients<T>& x, const RegressionCoefficients<T>& y) const {
const double pred_norm = (pred_norm_ < 0) ? std::min(arma::norm(data_->cx(), "inf"), arma::norm(data_->cx(), 1)) :
pred_norm_;
return std::sqrt(data_->n_obs()) * std::abs(x.intercept - y.intercept) +
pred_norm * arma::norm(x.beta - y.beta, 2);
}
//! Get the difference between two sets of regression coefficients.
//!
//! For the LS loss, the relative difference is an approximation to the 2-norm of the matrix-vector product
//! ||X . (beta1 - beta2) + ( mu1 - mu2 )||_2 < |mu1 - mu2| sqrt(n) + ||X|| ||beta1 - beta2||_2
//!
//! @param x a set of regression coefficients.
//! @param y the other set of regression coefficients.
//! @return the relative difference between `x` and `y`.
template<typename T>
double Difference(const RegressionCoefficients<T>& x, const RegressionCoefficients<T>& y) {
if (pred_norm_ < 0) {
pred_norm_ = std::min(arma::norm(data_->cx(), "inf"), arma::norm(data_->cx(), 1));
}
return std::sqrt(data_->n_obs()) * std::abs(x.intercept - y.intercept) +
pred_norm_ * arma::norm(x.beta - y.beta, 2);
}
//! Evaluate the gradient of the weighted LS loss at the given coefficient value.
//!
//! The gradient of the LS loss is given by
//! intercept = -1/n * sum_{i = 1}^n w_i (y_i - \hat{\mu} - x_i' . \beta)
//! beta = -1/n sum_{i = 1}^n w_i (y_i - \hat{\mu} - x_i' . \beta) x_i
//!
//! @param coefs the point where the gradient should be evaluated.
//! @param gradient_intercept a pointer where the gradient of the intercept should be stored at.
//! @param gradient_beta a pointer where the gradient of the slope should be stored at.
template<typename T>
GradientType<T> Gradient(const RegressionCoefficients<T>& where) const {
const arma::vec neg_residuals = data_->cx() * where.beta - data_->cy();
if (include_intercept_) {
return GradientType<T>(arma::mean(neg_residuals) + where.intercept,
arma::mean(data_->cx().each_col() % (neg_residuals + where.intercept), 0));
}
return GradientType<T>(0, arma::mean(data_->cx().each_col() % neg_residuals, 0));
}
//! Access the data used by this weighted LS loss function.
//!
//! @return a constant reference to this loss' data.
const PredictorResponseData& data() const noexcept {
return *data_;
}
//! Check if the intercept term is included in this weighted LS loss function.
//!
//! @return `true` if the intercept term is included, `false` otherwise.
bool IncludeIntercept() const noexcept {
return include_intercept_;