\documentclass[nojss]{jss}
%% -- LaTeX packages and custom commands ---------------------------------------

%% recommended packages
\usepackage{orcidlink,thumbpdf,lmodern}
\usepackage[utf8]{inputenc}
%% another package (only for this demo article)
\usepackage{tikz,bbm,amsmath,amssymb}
\usepackage{booktabs}
\usepackage{tabularx}

%% new custom commands
\newcommand{\class}[1]{`\code{#1}'}
\newcommand{\fct}[1]{\code{#1()}}


%% -- Article metainformation (author, title, ...) -----------------------------

%% - \author{} with primary affiliation (and optionally ORCID link)
%% - \Plainauthor{} without affiliations
%% - Separate authors by \And or \AND (in \author) or by comma (in \Plainauthor).
%% - \AND starts a new line, \And does not.
\author{
Laura Vana-G\"ur~\orcidlink{0000-0002-9613-7604}\\TU Wien
\And Roman Parzer~\orcidlink{0000-0003-0893-3190}\\TU Wien
\And Peter Filzmoser~\orcidlink{0000-0002-8014-4682}\\TU Wien
}
\Plainauthor{Laura Vana G\"ur, Roman Parzer, Peter Filzmoser}

%% - \title{} in title case
%% - \Plaintitle{} without LaTeX markup (if any)
%% - \Shorttitle{} with LaTeX markup (if any), used as running title
\title{\pkg{spareg}: Sparse Projected Averaged Regression in \proglang{R}}
\Plaintitle{spareg: Sparse Projected Averaged Regression in R}
\Shorttitle{Sparse Projected Averaged Regression in \proglang{R}}

%% - \Abstract{} almost as usual
\Abstract{
The \pkg{spareg} package for \proglang{R} builds ensembles of predictive
generalized linear models for high-dimensional data. It employs an algorithm
that combines variable screening and random projection techniques in each model
of the ensemble to address the computational challenges of large predictor sets.
By implementing this modeling approach in an accessible framework,
\pkg{spareg} enables users to apply methods that have shown competitive
predictive performance against state-of-the-art alternatives, while at the same
time keeping computational costs low.

Designed with extensibility in mind, \pkg{spareg} implements the screening and
random projection techniques, as well as the generalized linear models
employed in the ensemble as \proglang{S}3 classes with user-friendly
constructor functions. This modular design allows users to seamlessly
integrate and develop new procedures, making the package a versatile tool
for high-dimensional applications.
}

\Keywords{random projection, variable screening, ensemble learning, \proglang{R}}
\Plainkeywords{random projection, variable screening, ensemble learning, R}


\Address{
Laura Vana-G\"ur\\
Computational Statistics (CSTAT)\\
Institute of Statistics and Mathematical Methods in Economics\\
TU Wien\\
Wiedner Hauptstra{\ss}e 8-10\\
1040 Vienna, Austria\\
E-mail: \email{laura.vana.guer@tuwien.ac.at}


Roman Parzer\\
Computational Statistics (CSTAT)\\
Institute of Statistics and Mathematical Methods in Economics\\
TU Wien\\
Wiedner Hauptstra{\ss}e 8-10\\
1040 Vienna, Austria\\
E-mail: \email{romanparzer1@gmail.com}


Peter Filzmoser\\
Computational Statistics (CSTAT)\\
Institute of Statistics and Mathematical Methods in Economics\\
TU Wien\\
Wiedner Hauptstra{\ss}e 8-10\\
1040 Vienna, Austria\\
E-mail: \email{peter.filzmoser@tuwien.ac.at}
}
<<r setup, echo=FALSE>>=
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
library("ggplot2")
library("spareg")
is_paper <- FALSE
tmpdir <- tempdir()
@
\newif\ifPositive
<<condition_vignette, echo = FALSE, results = tex>>=
if (is_paper) {
  cat("\\Positivetrue")
} else {
  cat("\\Positivefalse")
}
@
% \VignetteIndexEntry{spareg: Sparse Projected Averaged Regression in R}
% \VignetteDepends{robustbase,cellWise,VariableScreening,R.matlab,doParallel,doRNG,ggpubr}
% \VignetteKeyword{random projection}
% \VignetteKeyword{screening}
% \VignetteKeyword{ensemble learning}
\begin{document}

\section{Introduction} \label{sec-intro}

The \pkg{spareg} package for \proglang{R} \citep{RLanguage} introduces a
computationally efficient framework for high-dimensional generalized linear
modeling, uniquely combining variable screening, random projections, and
ensemble learning to handle settings where the number of predictors $p$ far
exceeds the sample size $n$ i.e.,~$p>n$ or even $p\gg n$.
%%

More specifically, \pkg{spareg} builds an ensemble of generalized linear models
(GLMs) where, in each model of the ensemble,
i)~variables are first screened based on a measure of the utility of
each predictor for the response,
ii)~the selected variables are then projected to a lower dimensional space
(smaller than $p$) using a random projection matrix,
iii)~a GLM is estimated using the projected predictors.
Additional sparsity in the coefficients of the original variables
can be introduced through a thresholding parameter, which, together with the
number of models in the
ensemble, can be chosen using a validation set or via cross-validation.
Finally, predictions  are obtained by averaging over the marginal models in
the ensemble.

The motivation of such an algorithm, which performs what we call a
\emph{sparse projected averaged regression} (SPAR) for both discrete and
continuous data in the GLM framework, lies in its computational
efficiency.
Random projection is a computationally efficient method which linearly maps a
set of points in high dimensions into a much lower-dimensional space, and
random projection matrices have been traditionally generated from suitable
distributions in a data-agnostic fashion.
By projecting the original predictors to a dimension lower than $p$, estimation
of the models on the reduced predictors can be done using standard methods.
%
However, random projection can suffer from noise accumulation for very large
$p$, as too many irrelevant predictors are being considered for
prediction purposes \citep{Dunson2020TargRandProj}. Therefore, the screening step
is advisable in order to eliminate the influence of irrelevant variables before
performing the random projection, while also reducing computational costs.
%
On the other hand, the ensemble approach mitigates the
variability introduced by random sampling by averaging the results from multiple iterations \citep{Thanei2017RPforHDR}.

Different variants of the algorithm have been shown to perform
well in terms of prediction power on a variety of data sets.
For example, \cite{Dunson2020TargRandProj} employ the algorithm
with the data-agnostic, sparse random
projection of \cite{ACHLIOPTAS2003JL} combined with
probabilistic marginal correlation screening.
Furthermore, \cite{parzer2024glms} show that, when paired with a carefully constructed
data-driven random projection, the algorithm performs superiorly in terms of
predictions and variable ranking in settings with different degrees of
sparsity in the coefficients and with correlated predictors.

% \textcolor{red}{TODO: clarify that data-driven is good in certain settings, others in other settings etc. check dunson paper on theoretical results.}
Package \pkg{spareg} \citep{pkg:spareg}  implements this flexible framework,
offering a diverse selection of screening coefficients, random projection
matrices, and (un)penalized marginal GLMs, and
is available from the Comprehensive \proglang{R} Archive Network (CRAN) at
\url{https://CRAN.R-project.org/package=spareg}.
A key feature of the software is its extensibility:  users can tailor the
choice of screening, projection, and modeling techniques to their data's
specific characteristics.
For instance, the data-driven random projection proposed by \cite{parzer2024glms}
excels in settings with correlated predictors and unknown sparsity levels, and
it preserves information about the original coefficients, enhancing
interpretability. In contrast, data-agnostic random projections
\citep[e.g.,][]{ACHLIOPTAS2003JL} may lose structural information, making them
less suitable when interpretability is a priority. Similarly, certain screening
methods, such as marginal correlation screening, perform well when predictors
are weakly correlated, while others, like ridge-based screening, are more robust
in the presence of multicollinearity (more details are provided in
Section~\ref{sec-models}).
The package provides user-friendly constructor
functions for \proglang{S}3 classes, enabling users to integrate custom screening,
projection, or marginal GLM procedures into the ensemble framework,
as well as methods such as \fct{plot}, \fct{predict},
\fct{coef}, \fct{print} which allow users to more easily interact with the
model output and analyze the results.

% other software
In \proglang{R}, several
other packages focus on building ensembles where the dimensionality of the
predictors is reduced.
Most notably, package \pkg{RPEnsemble}  \citep{RPEnsembleR}
implements the procedure in \cite{cannings2017random}, where
``carefully selected'' random projections are used for projecting the predictors
before they are employed in a classifier such as $k$~nearest neighbor, linear or
quadratic discriminant analysis.
On the other hand, package \pkg{RaSEn} \citep{pkg:RaSEn} implements an
algorithm for ensemble classification and regression problems, where random
subsets of predictors are generated and the optimal one is chosen to train a
weak learner on the basis of some criterion.
Similarly, package \pkg{randomGLM} \citep{pkg:randomGLM} implements
an ensemble predictor based on bootstrap aggregation (bagging) of GLMs
whose covariates are first randomly selected and then further reduced
using forward regression according to the Akaike information criterion (AIC).
For \proglang{Python} \citep{Python}, a similar workflow could be set up
for building ensembles using the
package \pkg{pycobra} \citep{pkg:pycobra} by integrating \pkg{scikit-learn}
models \citep{pedregosa2011scikit}. Random projection tools
are implemented in \pkg{sklearn.random\_projection} while
GLM estimation on the projected predictors (both classical and penalized)
 can be done using \pkg{sklearn.linear\_model} module. Alternatively,
 ensembles can be built by using native stacking.
Since the recent integration of the \proglang{H2O} machine learning models
\citep{h2o} in \proglang{Stata} \citep{Stata}, it is also possible to build
ensembles of machine learning models (including GLMs) via stacking, boosting or
bagging, but random subspace methods or similar approaches to
dimensionality reduction in GLMs are not implemented at the time of writing.
%

Unlike existing tools, \pkg{spareg} integrates screening, random projection and
ensemble learning into a cohesive workflow, enabling estimation and
interpretation of GLMs in challenging high-dimensional scenarios. Its modular
design allows users to seamlessly extend the package with custom methods.
Furthermore, \pkg{spareg} is the first package to implement this specific
algorithm for GLMs, providing a
comprehensive toolbox for model fitting, prediction, and visualization.
Through the implemented data-driven random projection which ensures
interpretability of variable effects, the package is
particularly advantageous for practical applications in fields such as
genomics, imaging, and finance. By offering this integrated functionality,
the package fills a critical gap in the current software landscape, where no
other tool provides the same capabilities for GLMs in high-dimensional settings.

%
The rest of the paper is organized as follows:
Section~\ref{sec-models} provides an overview of the methods and the
methodological details of the implemented algorithm. The package is described
in Section~\ref{sec-software} while Section~\ref{sec-extensibility} exemplifies
how a new screening procedure, a new random projection and a new marginal model
can be integrated in the package.
\ifPositive
Section~\ref{sec-illustrations} contains two examples of employing the package
on real data sets.
\fi
Section~\ref{sec-conclusion} concludes.

\section{Methods} \label{sec-models}

Throughout the section we assume to observe
high-dimensional data $\{(\boldsymbol{x}_i,y_i)\}_{i=1}^n$,
where $\boldsymbol{x}_i\in\mathbb{R}^p$ is a predictor vector and
$y_i\in\mathbb{R}$ is the response, with $p\gg n$. The predictor vectors are
collected in the rows of the predictor matrix $X\in \mathbb R^{n\times p}$.

\subsection{Variable screening} \label{sec-varscreen}
%
In high-dimensional modeling, the goal of variable screening is to reduce the
predictor set by selecting a small subset of variables with a strong \emph{utility}
to the response variable. This initial selection enables more efficient
downstream analyses by discarding less relevant predictors early in the
modeling process, thus reducing computational costs and potential noise
accumulation stemming from irrelevant
variables \citep[see e.g.,][]{Dunson2020TargRandProj}.

The field of variable screening is an active area of research, with numerous
methods developed to address different data settings. While a comprehensive
review is beyond the scope of this paper, this section provides a concise and
selective overview of key approaches for variable screening in high-dimensional
settings.
%
A classic approach is sure independence screening (SIS),
proposed by \cite{Fan2007SISforUHD}. This method uses the vector of marginal
empirical correlations
$\hat{\boldsymbol\omega}=(\hat\omega_1,\ldots ,\hat\omega_p)^\top\in\mathbb{R}^p$,
where $\hat\omega_j=\text{Cor}(\boldsymbol X_{j},\boldsymbol y)$. Here,
$\boldsymbol y$ is the $(n\times 1)$ vector of responses and
$\boldsymbol X_{j}$ is the $j$-th column of the matrix of predictors.
SIS screens predictors in a linear regression setting by selecting the variable
set $\mathcal{A}_\delta$, which contains
all variables for which $|\hat\omega_j|>\delta$, with $\delta$ being a suitable
threshold.
%
Under certain technical conditions,  this screening coefficient vector has the
\emph{sure screening property}, i.e.,~the set of truly active variables
is included in $\mathcal{A}_{\delta_n}$ with probability converging to one as
$n\to \infty$.
Extensions to SIS include modifications for GLMs \citep{Fan2010sisglms}, where
screening is performed based on the log-likelihood $\ell(.)$ or the slope
coefficient of the GLM containing only $\boldsymbol X_{j}$ as a predictor:
$\hat\omega_j:= \text{argmin}_{\beta_j\in\mathbb{R}}\text{min}_{{\beta_0}\in
\mathbb{R}}\sum_{i=1}^n -\ell(\beta_j,\beta_0;y_i,x_{ij})$, where
$x_{ij}$ is the $j$-th entry of the vector $\boldsymbol x_i$.

However, both mentioned approaches face limitations related to
the required technical conditions which can rule out practically possible
scenarios where an important variable is marginally uncorrelated to the response
due to their multicollinearity. To tackle these issues,
\cite{fan2009ultrahigh} propose to
use an iterative procedure where SIS is applied subsequently on the residuals
of the model estimated in a previous step.
Additionally, in a linear regression setting, \cite{cho2012high} propose using the tilted correlation, i.e.,~the
correlation of a tilted version of $\boldsymbol X_{j}$ with $\boldsymbol y$ where the effect of other
variables is reduced. \cite{Wang2015HOLP}
propose to take into account the correlation among the predictors by
employing the  high-dimensional ordinary least squares projection (HOLP) $\hat{\boldsymbol w}=X^\top (XX^\top)^{-1}\boldsymbol y$, which is a ridge estimator where the penalty term converges to zero.  For GLMs, joint feature screening via a sparsity-restricted maximum likelihood estimator--under an $\ell_0$ norm--has been proposed in \cite{SMLE2014}.

In order to tackle potential model misspecification, a rich stream of literature
focuses on developing semi- or non-parametric alternatives to SIS. For linear
regression, approaches include using the ranked correlation \citep{zhu2011model},
(conditional) distance correlation \citep{li2012feature,wang2015conditional} or
quantile correlation \citep{ma2016robust}.
For GLMs, \cite{fan2011nonparametric} extend \cite{Fan2010sisglms} by fitting a
generalized additive model with B-splines. Further extensions for discrete
(or categorical) outcomes include the fused Kolmogorov filter
\citep{mai2013kolmogorov}, the mean conditional variance, i.e.,~the expectation
in $\boldsymbol X_{j}$ of the variance in the response of the empirical distribution function of each predictor conditional on the response \citep{cui2015model}.
\cite{ke2023sufficient} proposes a model-free method where the contribution of
each individual predictor is quantified marginally and conditionally in the
presence of the control variables as well as the other candidates by
reproducing-kernel-based $R^2$ and partial
$R^2$ statistics.

The \proglang{R} landscape for variable screening techniques
is very rich. An overview of some notable packages on the Comprehensive
\proglang{R} Archive Network (CRAN) includes the following packages.
Package \pkg{SIS} \citep{SISR}, which implements
the (iterative) sure independence screening procedure and its extensions, as
detailed in \cite{Fan2007SISforUHD},
\cite{Fan2010sisglms}, \cite{fan2010high}. This package
also provides functionality for estimating a penalized generalized linear model
or a cox regression model for the variables selected by the screening procedure.
Package \pkg{VariableScreening} \citep{pkg:VariableScreening} offers screening
methods for independent and identically distributed (iid) data,
varying-coefficient models, and longitudinal data, and includes techniques such
as sure independent ranking and screening (SIRS), which ranks the predictors by
their correlation with the rank-ordered response, or distance correlation sure
independence screening (DC-SIS), a non-parametric extension of the correlation
coefficient.
Package \pkg{MFSIS} \citep{pkg:MFSIS} provides a collection of model-free
screening techniques including SIRS, DC-SIS, the fused Kolmogorov filter
\citep{mai2015fusedkolmogorov}, the projection correlation method using knock-off
features \citep{liu2020knockoff}, among others. Additional packages implement
specific procedures but their review is beyond the scope of the current paper.

Package \pkg{spareg} allows the integration of such (advanced) screening
techniques using a flexible framework, which in turn enables users to apply
various screening methods tailored to their data characteristics in the
algorithm generating the ensemble. This flexibility allows users to evaluate
different strategies, ensuring that the most effective approach is chosen for
the specific application at hand. Moreover, it incorporates probabilistic
screening strategies, which can be
particularly useful in ensembles, as they enhance the diversity of predictors
across ensemble models. Instead of relying on a fixed threshold $\delta$,
predictors are sampled with probabilities
proportional to their screening coefficients
\citep[see][]{Dunson2020TargRandProj,parzer2024glms}.
Note that this is different from the random subspace sampling employed
in, e.g.,~random forests.
%
%
\subsection{Random projection tools} \label{sec-rps}
%
Package \pkg{spareg} has been designed to allow the incorporation of various
random projection techniques, enabling users to tailor the procedure to their
specific data needs. Below, we provide background information on different
random projection techniques and an overview of existing software implementing
random projections.

The random projection method relies on the Johnson-Lindenstrauss (JL) lemma
\citep{JohnsonLindenstrauss1984}, which asserts that for each set of
points in $p$-dimensional Euclidean space collected in the rows of
$X\in \mathbb{R}^{n\times p}$ there exists a linear map
$\Phi\in \mathbb{R}^{m \times p}$ such that all pairwise distances are
approximately preserved within a factor of $(1\pm\epsilon)$ for
$m\geq m_l=\mathcal O(\epsilon^{-2}\log(n))$.
Computationally, an attractive feature of the method for high-dimensional
settings is that the bound does not depend on $p$.

The goal is to choose a random map $\Phi$ that satisfies the JL lemma with high
probability given that it fulfills certain technical conditions. The literature
focuses on constructing such matrices either by sampling them from some
``appropriate'' distribution,
by inducing sparsity in the matrix and/or by employing specific fast constructs
which lead to efficient matrix-vector multiplications.
It turns out that the conditions are generally satisfied by
nearly all sub-Gaussian distributions \citep{matouvsek2008variants}.
A common choice is the standard normal
distribution $\Phi_{ij} \overset{iid}{\sim} N(0,1)$ \citep{FRANKL1988JLSphere}
or a sparser version where
$\Phi_{ij}\overset{iid}{\sim} N(0,1/\sqrt{\psi})$
with probability $\psi$
and $0$ otherwise \citep{matouvsek2008variants}.
Another computationally simpler option  is the Rademacher distribution where
$\Phi_{ij} =  \pm 1/\sqrt{\psi}$ with probability $\psi/2$ and zero otherwise for $0<\psi\leq 1$,
where \cite{ACHLIOPTAS2003JL} shows results for $\psi=1$ and $\psi=1/3$ while
\cite{LiHastie2006VerySparseRP} recommend using $\psi=1/\sqrt{p}$ to obtain very
sparse matrices.
%
Further approaches include using the Haar measure to generate random orthogonal
matrices \citep{cannings2017random} or
a non-sub-Gaussian distribution like the standard Cauchy,
proposed by \cite{LiHastie2006VerySparseRP}, for preserving approximate $\ell_1$
distances in settings where the data is high-dimensional, non-sparse, and heavy-tailed.
Structured matrices, which allow for more efficient multiplication, have also been proposed
\citep[see e.g.,][]{ailon2009fast,Clarkson2013LowRankApprox}.

The conventional random projections mentioned above are
data-agnostic. However, recent work has proposed incorporating information
from the data either to select the ``best'' data-agnostic random projection or
to directly inform the random projection procedure.
\cite{cannings2017random} rely on the former approach and
build an ensemble classifier where
the random projection matrix is chosen by
selecting the one that minimizes the test error of the
classification problem among a set of data-agnostic random projections.
%
Another data-driven approach to random projection for  regression has been
proposed by  \cite{ryder2019asymmetric}, who introduce a data-informed random
projection using an asymmetric transformation of the predictor matrix without
using information of the response.
%

On the other hand,  \cite{parzer2024glms} propose
to use a random projection matrix for GLMs which directly incorporates
information about the relationship between the predictors and
the response, rather than a projection matrix
which satisfies the JL lemma. In the linear regression case
\cite{parzer2024sparse} provide  a theoretical bound on the expected reduction
in prediction error when using a
projection which incorporates information
about the true regression coefficients compared to a conventional random projection.
In particular, they rely on the sparse embedding (CW) matrix of
\cite{Clarkson2013LowRankApprox}. The CW random projection matrix $\Phi$ is
constructed as
%
\begin{align}\label{eq:rp_cw}
\Phi=BD\in \mathbb{R}^{m\times p},
\end{align}
%
where $B$ is a
$(p\times p)$  binary matrix, where for each column~$j$ an index is uniformly
sampled from $\{1,\ldots,m\}$ and the corresponding entry is set to one, and
$D$ is a $(p\times p)$ diagonal matrix,
with entries $d_j \sim \text{Unif}(\{-1, 1\})$.
\cite{parzer2024sparse}  provide theoretical
results when replacing the diagonal entries of $D$ with the true regression
coefficients.
Motivated by this result, they propose to construct such a projection matrix
where the random diagonal elements are replaced in practice by
the HOLP estimates, i.e., the ridge coefficients with a penalty converging to
zero. The constructed random projection has the advantage of approximately
capturing the true regression coefficients in the span of the random projection,
i.e.,~it ensures that the true regression coefficients can be
recovered
approximately after the projection.
In \cite{parzer2024glms}, the authors extend this idea to the GLM framework by
investigating how similar ridge coefficients can be used to construct a
suitable data-informed projection matrix.
Through simulations, they demonstrate that for non-Gaussian families, allowing the penalty to converge to zero is not always optimal. Instead, they propose selecting the smallest penalty value such that the deviance ratio
(the proportion of the null deviance explained) remains below a specified threshold. Specifically, they suggest using a threshold of 0.8 for non-Gaussian families and 0.99 for Gaussian families.
%


Several packages in \proglang{R} provide functionality for random
projections. For instance, package \pkg{RandPro}
\citep{RandProR,SIDDHARTH2020100629} allows a Gaussian random matrix
(i.e.,~where entries are simulated from a standard normal), a sparse matrix
\citep{ACHLIOPTAS2003JL,LiHastie2006VerySparseRP}
or a matrix generated using the equal probability distribution with the
elements $\{-1,1\}$, to be applied to the predictor matrix before
employing one of $k$~nearest neighbors, support vector machine or naive
Bayes classifier on the projected features.
Package \pkg{SPCAvRP} \citep{SPCAvRPR}
implements sparse
principal component analysis, based on the aggregation of eigenvector
information from ``carefully-selected'' axis-aligned random projections of the
sample covariance matrix. Additionally,
package \pkg{RPEnsembleR} \citep{RPEnsembleR} implements a similar idea when
building the ensemble of classifiers: for each classifier in the ensemble,
a collection of
(Gaussian, axis-aligned projections, or Haar) random projection matrices is
generated, and the one that minimizes a risk measure
for classification on a test set is selected.
For \proglang{Python} \citep{Python} the \pkg{sklearn.random\_projection}
module \citep{pedregosa2011scikit} implements two types of unstructured random
matrices, namely Gaussian random matrix and sparse random matrix.
%
% <!-- \cite{Clarkson2013LowRankApproxShort} propose a sparse embedding matrix ${\Phi=BD}$, where $B\in\{0,1\}^{m \times p}$ is random binary matrix and $D$ is a $p\times p$ diagonal matrix with $(D_{ii}+1)/2\sim \text{Unif}\{0,1\}$, and prove that the dimension $m$ is bounded by -->
% <!-- a polynomial in $r\varepsilon^{-1}$ for $0 <\varepsilon< 1$ and $r$ being the rank of $X$. While this is generally larger than that of FJLT, the sparse embedding matrix requires less time to compute $\Phi X$ compared to other subspace embeddings. -->
%
%
%
%
% <!-- see more info at https://www.cs.waikato.ac.nz/~bobd/ECML_Tutorial/ECML_handouts.pdf. -->
%
% <!-- See also https://web.math.princeton.edu/~amits/publications/LeanWalsh_published.pdf. -->
%
% <!-- https://arxiv.org/pdf/1710.03163.pdf -->
%
% <!-- https://cseweb.ucsd.edu/~dasgupta/papers/randomf.pdf -->
%
% <!-- https://people.math.ethz.ch/~nicolai/mv/notes6.pdf -->
%
% <!-- https://web.math.princeton.edu/~amits/publications/LeanWalsh_published.pdf -->
%
% <!-- https://pastel.hal.science/tel-01481912/document -->
%
%
\subsection{Generalized linear models}
%
After we perform an
initial screening step followed
by a projection step in each marginal model,
we assume that the reduced and projected set of
predictors $\boldsymbol{z}_i=\Phi\boldsymbol{x}_i \in \mathbb{R}^m$ together
with the response arise from a GLM with
the response having conditional density from a (reproductive) exponential
dispersion family of the form
\begin{align*}\label{eqn:y_density}
f(y_i|\theta_i,\phi) = \exp\Bigl\{\frac{y_i\theta_i- b(\theta_i)}{a(\phi)} + c(y_i,\phi)\Bigr\},
\quad
g(\E[y_i|\boldsymbol{z}_i]) = \gamma_0 + (\Phi\boldsymbol{x}_i)^\top\boldsymbol{\gamma}=:\eta_i,
\end{align*}
where $\theta_i$ is the natural parameter, $\phi$ is a dispersion
parameter,  $a(.)>0$ and $c(.)$ are specific
real-valued functions determining different families, and
$b(.)$ is the log-partition function normalizing the density to integrate to one.
If $\phi$ is known, we obtain densities in the natural exponential family for
our responses. The responses are related to the $m$-dimensional reduced and
projected predictors through the conditional mean, i.e.,~the conditional mean
of $y_i$ given ${\boldsymbol{z}}_i$
depends on a linear combination of the predictors through (invertible) link
function $g(.)$, where $\gamma_0\in\mathbb{R}$ is the intercept and
$\boldsymbol{\gamma}\in\mathbb{R}^m$ is a vector of regression coefficients
for the $m$ projected predictors. We can see that the original coefficients of the
predictors $\boldsymbol{x}_i$ can be obtained as: $\boldsymbol{\beta}=\Phi^\top\boldsymbol{\gamma}$.


Estimates for $\gamma_0\in\mathbb{R}$ and $\boldsymbol{\gamma}\in\mathbb{R}^m$
can be obtained by directly maximizing the log-likelihood below, as generally one chooses $m < n$:
\begin{align*}
  \ell(\beta_0,\boldsymbol{\beta}) = \sum_{i=1}^n \frac{y_i\theta_i(\beta_0,\boldsymbol{\beta},\boldsymbol{x}_i)- b(\theta_i(\beta_0,\boldsymbol{\beta},\boldsymbol{x}_i))}{a(\phi)} + c(y_i,\phi),
\end{align*}
where $\theta_i(\beta_0,\boldsymbol{\beta},\boldsymbol{x}_i) =
(b')^{-1}(g^{-1}(\eta_i))$.
However, even if $m < n$
it might still be desirable to add a (e.g.,~small $\ell_2$) penalty to the
likelihood of the marginal models--for the binomial family this can alleviate
problems related to separation while keeping the bias low. Moreover, if outlier
influence should be reduced,
$M$- or trimmed estimators could be employed for obtaining  robust
estimates of the regression coefficients
\citep[see e.g.,][]{cantoni2001robust}.

%
% <!-- Therefore, we propose to employ marginal models which are estimated by either  -->
% <!-- maximizing the log-likelihood or a penalized version thereof: -->
% <!-- $$ -->
% <!-- \text{argmin}_{{\gamma}\in\mathbb{R}^m}\min_{\gamma_0\in\mathbb{R}}  \sum_{i=1}^n -\ell(\gamma_0, \gamma;y_i,\boldsymbol{z}_i) + \frac{\varepsilon}{2}\sum_{j=1}^m{\gamma}_j^2, \, \varepsilon \geq 0. -->
% <!-- $$ -->
% <!-- However, the framework can in principle accommodate for other penalty functions.  -->
%
\subsection{SPAR algorithm} \label{sec-algo}

We present the general algorithm for sparse projected averaged regression (SPAR)
implemented in package \pkg{spareg}. A diagram of the algorithm is shown in
Figure~\ref{fig:algo}.
%%
\ifPositive
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=\textwidth}
\includegraphics{spareg-flowchart4th-v2}
\caption{Flowchart of the SPAR algorithm.\label{fig:algo}}
\end{figure}
\fi
%%
\begin{enumerate}
%%
\item  Choose family with corresponding log-likelihood $\ell(.)$ and link $g(.)$.
%%
\item  Standardize the $(n\times p)$ matrix of predictors $X$ for all families
and the vector of responses $\boldsymbol y$ for the Gaussian family by subtracting the
sample mean and dividing by the sample standard deviation of each variable.
%%
\item  Calculate screening coefficients $\hat{\boldsymbol \omega}$.
%%
\item  For $k=1,\dots,M$ models:
\begin{enumerate}
\item   If $p>p_0$, where $p_0$ is the number of variables to be screened, screen $p_0$ predictors based on the screening coefficients
$\hat{\boldsymbol\omega}$, which yields for model $k$ the screening index set
$I_k=\{j_1^k,\dots,j_{p_0}^k\}\subset \{1,\dots,p\}$.
If probabilistic screening should
be employed, draw the predictors sequentially without replacement using an
initial vector of probabilities $\boldsymbol{\pi}$ with $\pi_j\propto |\hat\omega_j|$; otherwise,
select the $p_0$ variables with the largest $|\hat\omega_j|$. If $p \leq p_0$,
perform no screening and set $I_k=\{1,\dots,p\}$.
\item  Project the screened variables to a randomly drawn dimension
$m_k\sim \text{Unif}(\{m_l,\dots,m_u\})$ using \emph{projection matrix}
$\Phi_k$ to obtain $Z_k=X_{\cdot I_k}\Phi_k^\top \in \mathbb{R}^{n\times m_k}$,
where $X_{\cdot I_k}$ contains the columns in $X$ having a column index in
$I_k$.
\item  Fit a \emph{GLM}-type model of $\boldsymbol y$ on $Z_k$ to obtain estimated
coefficients $\widehat{\boldsymbol\gamma}^k\in\mathbb{R}^{m_k}$ and
$\hat {\boldsymbol\beta}_{I_k}^k=\Phi_k^\top\widehat
{\boldsymbol\gamma}^k$, $\hat {\boldsymbol\beta}_{\bar I_k}^k=0$, where
$\bar I_k=\{1,\dots,p\}\setminus I_k$.
\end{enumerate}
\item  For a given threshold $\nu\geq 0$, set all $\hat\beta_j^k$ with $|\hat\beta_j^k|\leq\nu$ to $0$ for all $j,k$.
\item  Choose $M$ and $\nu$ via a validation set or cross-validation
evaluating a two-dimensional
grid of values. If cross-validation is performed,  repeat Steps~2 to 6
for each training fold using the fixed
index set $I_k$ and projection matrix $\Phi_k$, and assess the performance on the
test fold by a  loss function $K(M, \nu)$:
\begin{align*}
(M_{\text{best}},\nu_{\text{best}}) = \text{argmin}_{M,\nu}K(M,\nu).
\end{align*}
\item Combine models of the ensembles either via the coefficients by using a
\emph{simple average} of the regression coefficients
$\hat{\boldsymbol\beta} = \sum_{k=1}^M\hat {\boldsymbol\beta}^k / M$ (this results in averaging of the
models on the \emph{link level}) or via the fitted values by taking a
\emph{simple average} on the \emph{response level} i.e.,~$\hat y = \sum_{k=1}^Mg^{-1}(\boldsymbol x_i^\top\hat{\boldsymbol\beta}^k) / M$ (here we omit the intercept for
the sake of notational simplicity).
\item Output the averaged estimated coefficients and predictions for
$(M_{\text{best}},\nu_{\text{best}})$.
\end{enumerate}
%
The proposed algorithm is flexible, allowing for various choices of
screening coefficients $\hat{\boldsymbol w}$, random projection matrices
$\Phi$ and marginal models. Several variants have been introduced in the literature.
For example, \cite{Dunson2020TargRandProj} propose a version
that uses probabilistic marginal correlation screening along with the sparse random projection method from \cite{ACHLIOPTAS2003JL}, without applying thresholding.
More recent work by \cite{parzer2024sparse,parzer2024glms} extends this approach by employing data-informed random projection that leverages information from the ridge coefficients and complement the projection step with
probabilistic screening based on the same ridge regression coefficients
(see details in Section~\ref{sec-rps}).

%
In order to choose default values for
$p_0$, $m_l$, $m_u$ as well as the maximal number of considered models $M$,
we rely on previous literature.
%%
Regarding the number of screened variables, \cite{parzer2024sparse}
propose choosing $p_0=c\cdot n$ as a multiple of $n$ (thus independent of $p$) and find that the best results are achieved for $2\leq c\leq 4$. We thus set in the algorithm $p_0=2n$ as a default value.

The dimension $m_k$ is random for the purpose of introducing more variability
in the ensemble and reducing the reliance on a fixed (possibly arbitrarily chosen)
goal dimension. \cite{Dunson2020TargRandProj} propose using $m_l=2\log(p)$ and $m_u=3n/4$
while \cite{parzer2024glms} use slightly smaller goal dimensions ($m_l=\log(p)$, $m_u=n/2$), to reduce the dimension of the marginal models to be estimated.
On the other hand, for random projection matrices satisfying the JL lemma, $m_l$ can be derived from the theoretical bounds for preserving the distances between all pairs of points approximately. We employ as default values the ones proposed in \cite{parzer2024sparse}.

A further modification of the algorithm above is to
employ part of the data for estimating the screening coefficients (Step~3) and the remaining data for estimating the ensemble models (Step~4). Such a strategy could avoid issues related to overfitting. Even though \cite{parzer2024sparse} find that in the linear regression case such a splitting approach does not improve performance, the implementation in \pkg{spareg}
allows for data splitting.

\cite{parzer2024sparse} find that, when using their proposed data-driven
random projection and HOLP screening coefficients in a linear regression setting,
the gain in predictive performance when using more than $M=20$ models is marginal.
We use therefore $M=20$ as a default value in the algorithm. However, for ultra
high-dimensional settings, larger values of $M$ can be beneficial to explore
more of the predictor space, particularly when using data-agnostic projections.
In such cases, stronger thresholding (higher $\nu$) is also recommended to
mitigate overfitting, especially when no screening is applied or when the
expected sparsity is medium to high. When using data-agnostic
projections, we encourage users to explore a grid of $M$
up to 50 or 100 to encourage diversity in the ensemble.

The default grid for $\nu$, based on equally spaced quantiles of the absolute
values of the estimated coefficients, generally yields reasonable performance.
However, in high-dimensional settings with high expected sparsity, thresholding
becomes particularly useful for larger $n$ values and especially when employing the CW data-dependent projection with ridge coefficients, as it helps identify zero coefficients in the recovered signal.

Finally, when choosing the optimal $M$ and $\nu$ in Step~6 given a grid of values
for each parameter, one can pick the combination delivering the lowest loss
function or the combination delivering
the sparsest model with a loss within one standard deviation of the minimum
(this corresponds to the commonly employed ``one-standard-error (1se)'' rule
in penalized regression). We note that the cross-validation on the grid of
number of models is not computationally intensive as we only need to train the
maximum number of models considered in the grid on each training fold and
evaluate performance by aggregating the first $M_1$ models, then the first $M_2$, and so on.

\section{Software} \label{sec-software}
%
Package \pkg{spareg} is available on CRAN and can be installed and loaded as folows:
<<eval=FALSE>>=
install.packages("spareg")
library("spareg")
@
<<echo=FALSE>>=
library("spareg")
@
%
An overview of the main functions is provided in Table~\ref{tab:overview},
followed by detailed descriptions in the remainder of this section. For
demonstration, we use a simulated dataset with a continuous response
and $p=2000$ predictors where $n=200$ observations are
used as training and
$n=100$ observations are used as test data. The data is generated using
\code{simulate_spareg_data()}, which simulates by default a linear regression
model with $\sigma^2=83$, intercept $\mu=1$ and
coefficients $\boldsymbol\beta$
with 100 non-zero entries uniformly
sampled from $\{-3,-2,-1,1,2,3\}$:
<<eval=FALSE>>=
set.seed(1234)
example_data <- simulate_spareg_data(n = 200, p = 2000, ntest = 200)
str(example_data)
@
<<echo=FALSE>>=
set.seed(1234)
example_data <- simulate_spareg_data(n = 200, p = 2000, ntest = 200)
str(example_data, list.len = 4)
@
\begin{table}[t!]
\centering
\begin{footnotesize}
\begin{tabularx}{\textwidth}{p{3cm} X}
\toprule
\textbf{Function} & \textbf{Description} \\
\midrule
\multicolumn{2}{l}{\emph{High-level functions--return \proglang{S}3  objects of class \class{spar} and \class{spar.cv}}} \\
% \midrule
\code{spar()} &
Fits the ensemble. Supports parallelization and validation-set tuning for $M$ and $\nu$. \\
\code{spar.cv()} &
Cross-validated wrapper for optimal $M$ and $\nu$ selection. \\
\midrule
\multicolumn{2}{l}{\emph{Screening functions \code{screen\_*(..., control = list())}--return \proglang{S}3  objects of class \class{screencoef}}} \\
%%
\addlinespace[1ex]
%%
\multicolumn{2}{l}{\emph{Common attributes to be passed by \code{...}}}\\
\multicolumn{2}{l}{\hspace{5ex}\code{nscreen}: Variables to retain.  \emph{Default:} $\min(2n, p)$}\\
\multicolumn{2}{l}{\hspace{5ex}\code{split\_data\_prop}: Proportion of data for screening.  \emph{Default:} 1.0 (full)} \\
\multicolumn{2}{l}{\hspace{5ex}\code{type}: Screening logic (probabilistic \code{"prob"} or \code{"fixed"}).  \emph{Default:} \code{"prob"}} \\
%%
\addlinespace[1ex]
%%
\code{screen\_marglik()} &
Univariate marginal likelihood screening using \code{stats::glm()}. \\
\code{screen\_cor()} &
Correlation-based screening using \code{stats::cor}; accepts \code{method} for correlation type in \code{control} argument.  \emph{Default:} \code{method = "pearson"}. \\
\textbf{\code{screen\_glmnet()}} &
Multivariate screening using \code{glmnet::glmnet} and its controls.  \emph{Default:} elastic net mixing parameter \code{alpha = 0} (ridge regression). \\
\midrule
\multicolumn{2}{l}{\emph{Random projection functions \code{rp\_*(..., control = list())}--return \proglang{S}3 objects of class \class{randomprojection}}} \\
%%
\addlinespace[1ex]
%%
\multicolumn{2}{l}{\emph{Common attributes to be passed by \code{...}}}\\
\multicolumn{2}{l}{\hspace{5ex}\code{mslow}: Minimum projection dimension. \emph{Default:} $\log(p)$} \\
\multicolumn{2}{l}{\hspace{5ex}\code{msup}: Maximum projection dimension. \emph{Default:} $\min(n/2, p)$ } \\
%%
\addlinespace[1ex]
%%
\code{rp\_gaussian()} &
Projection with iid normal entries using \code{stats::rnorm()}. \emph{Default:} \code{mean = 0, sd = 1}. \\
\code{rp\_sparse()} &
Sparse projection with entries from the Rademacher distribution (parameter \code{psi} passed through \code{control}). \emph{Default:}  \code{psi = 1}. \\
\code{rp\_cw()}$^\dagger$ &
Sparse embedding matrix. \emph{Default:} \code{rp\_cw(data = TRUE)} for data-driven diagonals, otherwise randomly generated. \\
\midrule
\multicolumn{2}{l}{\emph{Marginal model functions \code{spar\_*(..., control = list())}--return \proglang{S}3  objects of class \class{sparmodel}}} \\
%%
\addlinespace[1ex]
%%
\code{spar\_glm()} &
Fits unregularized GLM for each projected model using \code{stats::glm};
accepts \code{glm} controls in \code{control} argument).\\
\code{spar\_glmnet()}$^\dagger$ &
Fits marginal elastic net models using \code{glmnet::glmnet};
accepts \code{glmnet} controls in \code{control} argument).  \emph{Default:}
small ridge penalty. \\
\midrule
\multicolumn{2}{l}{\emph{\proglang{S}3 methods for \class{spar} and \class{spar.cv} objects}} \\
% \midrule
\code{coef()} &
Extracts raw or aggregated coefficients as an \proglang{S}3  \class{coefspar} object for specific
$(M, \nu)$.  \emph{Default:} mean aggregation for best $(M, \nu)$. \\
\code{predict()} &
Generates predictions aggregated on \code{"link"} or \code{"response"} levels for specific $(M, \nu)$.  \emph{Default:} link-level mean aggregation for best $(M, \nu)$. \\
\code{plot()} &
Returns a \class{ggplot} of (cross-)validation measure, sparsity, residuals  or raw coefficients.
\emph{Default:} (cross-)validation measure as a function of $\nu$. \\
\midrule
\multicolumn{2}{l}{\emph{Extractor and constructor functions}} \\
% \midrule
\code{get\_intercept()} &
Extracts intercept(s) from a \class{coefspar} object. \\
\code{get\_coef()} &
Extracts the (aggregated) coefficients from \class{coefspar} object. \\
\code{get\_model()} &
Retrieves the $(M,\nu)$ model minimizing the (cross-)validation measure; supports one-standard-error choice for \class{spar.cv}. \\
\code{get\_measure()} &
Extracts performance measure and sparsity statistics across the tuning grid from \class{spar} and \class{spar.cv} objects. \\
\code{constructor\_*()} &
Helpers to create \class{screencoef}, \class{randomprojection}, or \class{sparmodel} objects. \\
\bottomrule
\end{tabularx}
\end{footnotesize}
\caption{Summary of key functions in the \pkg{spareg} package. $^\dagger$ marks
the modular components used as defaults in the \fct{spar} and \fct{spar.cv}
functions. By default, no screening is performed.}
\label{tab:overview}
\end{table}

\subsection{Main functions and their arguments}
% %
The two main functions for fitting the SPAR algorithm are:
%
\begin{itemize}
\item \fct{spar} fits the SPAR algorithm without
cross-validation and returns an \class{spar} object.
\item \fct{spar.cv} implements the cross-validated procedure and returns a `\code{spar.cv}' object.
\end{itemize}
%
Additionally, the aliases \fct{spareg} and \fct{spareg.cv} are available for convenience.
Both functions accept the following key arguments: \\
\begin{itemize}
    \item \code{x}--an \( n \times p \) numeric matrix of predictors and
    \code{y}--a numeric vector of length \( n \).
%
    \item \textbf{\code{family}}--an object from \fct{stats::family}; defaults to \fct{gaussian}.

    \item \textbf{\code{model}}--an object of class \class{sparmodel} specifying the model for each element of the ensemble. Defaults to \code{spar\_glm()} for Gaussian-identity and to \code{spar\_glmnet()} for others.

    \item \textbf{\code{rp}}--an object of class \class{randomprojection}. Defaults to \code{rp\_cw(data = TRUE)}.

    \item \textbf{\code{screencoef}}--an object of class \class{screencoef}. Defaults to \code{NULL} (no screening).

    \item \textbf{\code{nnu}}--number of threshold values \( \nu \) to consider. Defaults to 20.

    \item \textbf{\code{nus}}--an optional vector of threshold values \( \nu \). If not provided, defaults to a grid of \code{nnu} values: zero and \code{nnu - 1} equally spaced quantiles of the absolute values of the non-zero marginal coefficients.

    \item \textbf{\code{nummods}}--number of models in the ensemble. Defaults to 20. Also allows for a vector of ensemble sizes.

    \item \textbf{\code{measure}}--performance measure for selecting $M$ and $\nu$. Options: \code{"deviance"} (default), \code{"mse"} (mean squared error), \code{"mae"} (mean absolute error), \code{"class"} (misclassification error), or \code{"1-auc"} (one minus area under the ROC curve).

    \item \textbf{\code{avg\_type}}--stage of aggregation for computing the performance measure---either on \code{"link"} or \code{"response"} level.

    \item \textbf{\code{parallel}}--use parallel backend if available. Defaults to \code{FALSE}.
\end{itemize}
%
Furthermore, \code{spar()} has the arguments \code{xval} and  \code{yval}--validation data for choosing  $M_\text{best}$ and $\nu_\text{best}$ (if not provided, \code{x} and
\code{y} will be employed). Optional argument \code{inds} can be used to pass a list specifying columns of variables to retain in each marginal model
while  projection matrices to be used in each marginal model can be passed as a list through optional argument \code{RPMs}.
%

Function~\code{spar.cv()} has specific arguments
\code{nfolds}, the
number of folds to be used (defaults to 10) and
\code{precompute_mode = c("precompute_all", "precompute_rpm", "full_cv")}. Argument
\code{precompute_mode} indicates whether to use the same random projection matrix and/or
the same indices for screening across folds
(defaults to \code{"precompute_all"}).
%
If one of the first two options is selected,  \code{spar.cv()} uses a
lightweight version of \code{spar()} in each fold. Here, random projection
matrices and/or screening indicators are computed once on
the full dataset before cross-validation.
This is
achieved by an initial call to \code{spar()} on the entire dataset.
If \code{"precompute_all"}, screening indicators remain fixed throughout the
analysis. The random elements of the projection matrices are also held
fixed in all folds for \code{"precompute_all"} and \code{"precompute_rpm"}
However, if data-driven projections are used, the training data for each
fold will be used to update the data-dependent entries of the projection matrix
in all cases.
If \code{precompute_mode = "full_cv"} is selected, the random projection and screening
indicators are recomputed in each fold.
Fixing the screening indicators and/or (random elements of)
projection matrices prior to
cross-validation reduces computational overhead in some settings. On the other
hand, fixing these components removes a source of variability from the
cross-validation process which can lead to a more stable and interpretable
selection of $M$ and $\nu$, as the cross-validation results are not confounded
by variations in the screening or projection steps. However, this approach can
be affected by target leakage, as estimates obtained on the whole data are used
in the cross-validation exercise. This is particularly relevant when performing
screening.
To avoid target leakage \code{precompute_mode = "full_cv"} or
\code{precompute_mode = "precompute_rpm"}
can be selected.
%
\subsection[Extensible S3 objects]{Extensible \proglang{S}3 objects}\label{sec:soft:S3}
The \code{spar()} framework is organized around a common extensible structure comprising three method-specification objects:
\class{screencoef} for screening coefficients, \class{randomprojection} for
random projection matrices and \code{sparmodel} for the marginal model to be
used in the ensemble.

We will explain in the following how each object is created to store the
selected method and its initial \code{control} parameters. These objects are
created before data is passed to \code{spar()}. Once the relevant training data
and response family are available, \code{spar()} invokes the object's
\code{update_fun()} to resolve context- or data-dependent settings, adjust
controls, and store any required state;
it then invokes the corresponding generation function \code{generate_fun()} to
perform the primary computation.
As a convention, \code{control} is a list of parameters to be passed to the main
generating function, while parameters related to the procedure itself (such as
the number of variables to be screened or the minimum dimension to which the
variables should be projected) are to be saved as attributes to the objects.
The shared interfaces allow users to extend the framework with custom screening,
projection, and modeling methods while preserving consistent behavior across the
three object types.

\subsubsection{Screening coefficients}\label{sec:soft:screen}
%
Screening coefficients in \pkg{spareg} are implemented as \proglang{S}3 classes \class{screencoef}.
These objects are created using functions like \code{screen_*()} which accept as arguments:
\code{...} (i.e.,~ellipsis) and
\code{control} (a list of parameters for computing the screening coefficients).
Additional arguments passed through the ellipsis will be saved as
attributes of the object.

%
The following screening coefficients are implemented in \pkg{spareg}:
%
\begin{itemize}
\item \code{screen_marglik()}--computes the screening coefficients by the
coefficient of $\boldsymbol X_{j}$ for $j =1,\dots,p$  in a univariate GLM using the
\fct{stats::glm} function.
It allows to pass a list of controls through the \code{control} argument to
\fct{stats::glm}
e.g.,~\code{screen_marglik(control = list(family = binomial(probit)))}.
%
\item  \code{screen_cor()}--computes the screening coefficients by the
correlation between $\boldsymbol y$ and $\boldsymbol X_{j}$ using the function \fct{stats::cor}. It
allows to pass a list of controls through the \code{control} argument to
\fct{stats::cor}--e.g.,~\code{screen_cor(control = list(method = "spearman"))}.
%
\item  \code{screen_glmnet()}--computes by default a ridge coefficient with a small penalty:
$$
\hat\omega:= \text{argmin}_{\beta\in \mathbb{R}^p}\text{min}_{{\beta_0}
\in\mathbb{R}}\sum_{i=1}^n -\ell(\beta;y_i,\boldsymbol x_i) + \frac{\varepsilon}{2}
\sum_{j=1}^p{\beta}_j^2, \, \varepsilon > 0.
$$
This function uses \code{glmnet()} of package
\pkg{glmnet} \citep{pkg:glmnet, pkg:glmnet_glm} with
elastic net mixing parameter $\alpha =0$ (ridge penalty),
a small \code{lambda.min.ratio} (smallest penalty as a fraction of the maximum value in the grid) and
penalty selection based on \cite{parzer2024glms}: the smallest penalty where the deviance ratio is
$< 0.99$ (Gaussian) or $< 0.8$ (other families).
Custom controls can be passed  via the \code{control} argument
e.g.,~\code{screen_glmnet(control = list(alpha = 0.5))}.
\end{itemize}
%
The following attributes can be passed to the \code{screen_*()} function through \code{...}:
%
\begin{itemize}
\item \code{nscreen}--number of variables to be retained
after screening; defaults to $2n$. If
$2n > p$, no screening is performed.
%
\item \code{split_data_prop}--proportion of the data to be used for
computing the screening coefficient. The remaining data will be used for
estimating the marginal models; default is to use the whole data for
screening coefficients and the marginal model estimation.
%
\item \code{type}--either \code{"prob"} (for
probabilistic screening--default)  or \code{"fixed"} (for using a fixed set of
\code{nscreen} variables).
%
\item \code{reuse_in_rp}--indicates whether the screening
coefficients should be reused at a later stage in the construction of a data-dependent random
projection. Defaults to \code{FALSE}. This attribute is ignored if a data-agnostic random projection is used.
\end{itemize}
%
%
All implemented \code{screen_*()} functions return an object of class
\class{screencoef}
which in turn is a list containing an \code{update_fun}
function which in the default case updates the object with the \code{family} information,
a \code{generate_fun}
function, the \code{control} list passed by the user in
\code{screen_*()}, and an optional \code{name} string. A user-friendly \fct{print} method of the \class{screencoef} object is also provided.
% %
% <<>>=
% screen_marglik()
% @
% \begin{itemize}
% \item an optional character \code{name} used for printing;
% %
% \item \code{generate_fun()}--function to generate a vector of screening coefficients of length $p$. Takes \code{x} (standardized predictors), \code{y}  (standardized responses), and \code{object} (the \class{screencoef} object itself).
% %
% \item  \code{control}, which is the control list passed by the user in
% \code{screen_*()}. These controls are arguments to be used in
% \code{generate_fun()}.
% \end{itemize}
%
% The structure of the object is the following:
% <<>>=
% str(obj)
% @
% Function \code{generate_fun()} defines the generation of the screening
% coefficient. Note that it considers the controls in \code{object$control} when
% calling the \code{stats::glm()} function (unless it is provided, the
% \code{family} argument in \code{stats::glm()} will be set to the ``global''
% family of the SPAR algorithm which is assigned
% inside the \code{spar()} function as an attribute for the \class{screencoef} object).
For convenience, \code{constructor\_screencoef()}, a constructor function
is provided, which can be used to create new \code{screen_*()} functions.
An example is presented in Section~\ref{sec-extensscrcoef}.
%
%
\subsubsection{Random projections}\label{sec:soft:rps}
%
Random projection objects are
implemented as \proglang{S}3
classes \class{randomprojection} and are created by functions like \code{rp_*()}.
These functions take the ellipsis \code{...} and a list of controls \code{control} as arguments.
%
% \begin{Code}
% rp_*(..., control = list())
% \end{Code}
%
The package includes the following random projections:
\begin{itemize}
\item \code{rp_gaussian()}--entries are iid from the normal distribution (defaults to standard normal).
%
\item \code{rp_sparse()}--sparse matrix as in \cite{ACHLIOPTAS2003JL};
$\psi$ can be specified
via \code{control} e.g.,~\code{rp_sparse(control = list(psi = 1/3))}.  Defaults to  \code{psi = 1}.
%
\item \code{rp_cw()}--sparse embedding random projection in
\cite{Clarkson2013LowRankApprox}.
If the random projection is specified as \code{rp_cw(data = TRUE)} and the
screening coefficient object has the attribute \code{reuse_in_rp = TRUE}, the
diagonal entries of $D$ in Equation~\ref{eq:rp_cw}
are set to the screening coefficients.
Otherwise, if \code{rp_cw(data = TRUE)}, the data-driven modification to the CW matrix proposed in \cite{parzer2024glms} is used.
\end{itemize}
%
Common attributes which can be passed via \code{...} to the \code{rp_*()} functions include:
%
\begin{itemize}
\item \code{mslow}--minimum dimension for projecting the predictors; defaults to
$\log(p)$. Note that for JL-type random projections, \code{mslow} can be chosen
from JL bounds.
%
\item \code{msup}--maximum dimension for projecting the predictors; defaults to
$n/2$ or. if $p < n/2$, to $p$.
%
\end{itemize}
%
% Note that for JL-type random projections, \code{mslow} can be chosen from JL bounds
% which preserve the distances between
% all pairs of points within a factor $(1 \pm \epsilon)$ (e.g.,~\cite{ACHLIOPTAS2003JL} show that for a given $\epsilon\in(0,1)$, the JL lemma is satisfied for
% $m_0 = \log n(4 + 2\tau)/(\epsilon^2/2 - \epsilon^3/3)$ with probability
% $1 - n^{-\tau}$).
%
The \code{rp_*()} functions return an object of class \class{randomprojection}
which is a list of five elements. A constructor function
\code{constructor\_randomprojection()} for creating these objects is provided
(as will be shown in Section~\ref{sec-ext-rp}).
The most important elements contain
\code{generate_fun()}, the function for generating the random projection
matrix, and \code{control}, the list of controls passed by the user in
\code{rp_*()}.
As previously mentioned, objects of class \class{randomprojection}
are constructed before the data are available. Consequently,
settings that depend on the data cannot always be resolved by the
constructor. Before calling \code{generate_fun()}, \fct{spar} calls
\code{update_fun} to resolve such settings and to store any (data-dependent)
information in the object (such as \class{family} or precomputed screening coefficients). If argument \code{update_fun} is not specified when using a constructor, the default will only add the family information to the  \class{randomprojection} object.
Further optional elements include \code{name}, or function
\code{update_rpm_w_data()}, which specifies how to update the pre-computed
matrices (when \code{precompute_mode} is \code{precompute_all} or \code{precompute_rpm})  with data information in each fold of the cross-validation procedure.
Note that \code{update_fun()} updates the RP specification before
generating matrices while \code{update_rpm_w_data()} modifies an already-generated matrix, usually during cross-validation.
% \begin{itemize}
% \item a optional character \code{name}
% for printing;
% %
% \item \code{generate_fun()} function for generating the random projection matrix . This
% function should have arguments
% \begin{itemize}
% \item \code{rp}, which is itself a \class{randomprojection}
% object;
% \item \code{m}, the target dimension;
% \item a vector of indices \code{included_vector}
% which indicates the column index of the original variables in the \code{x}
% matrix to be projected using the random projection. This is needed due to the
% fact that screening can be employed pre-projection.
% \item \code{x} and \code{y}, the matrix of standardized
% predictors and vector of  (standardized)
% responses--can be \code{NULL} if the random projection to be generated is data-agnostic.
% \end{itemize}
% It returns a matrix or
% a sparse matrix of class \class{dgCMatrix} of the \pkg{Matrix} package
% \citep{pkg:Matrix} with \code{m} rows and \code{length(included_vector)} columns.
% %
% \item \code{control}, which is the control list in \code{rp_*()}. These controls
% are arguments needed in \code{generate_fun()} in order to generate the desired
% random projection.
% \end{itemize}
% For the case where the random projection should incorporate some information
% related to the data, two more elements can be potentially relevant.
% \begin{itemize}
% \item Function \code{update_fun()} updates the \class{randomprojection}
% object with relevant information from the arguments of the \code{spar()} or \code{spar.cv()}
% function call. If it is not provided, it defaults to updating the
% \class{randomprojection} object with a further attribute
% \code{family_string} which is character version of the \code{family} argument
% (e.g.,~in the default case it will be \code{"gaussian(identity)"}).
% In constructions like those of \cite{parzer2024glms},
% certain data-dependent quantities--such as ridge coefficients--are estimated once at the start and reused across all projections. To support this, the \code{update_fun()} can be used to attach data information as attributes to the \class{randomprojection} object, which are then accessed by \code{generate_fun()}. If provided, \code{update_fun()} should accept \code{...} as its argument; the \code{spar()} function will then internally pass all its own arguments to it. The function must return a \class{randomprojection} object.
% %
% \item In the case where a list of predefined \code{RPMs} is provided in \code{spar()},
% for the data-driven random projections we allow the user to
% specify whether some parts of the given \code{RPMs} should be updated with
% the provided data. This is possible through another optional function
% \code{update_rpm_w_data()}.
% This is particularly relevant for the cross-validation procedure,
% which employs the random projection matrices
% generated by calling the \code{spar()} function on the whole data set
% before starting the cross-validation exercise. For example, in our
% implementation of the data-driven \code{rp_cw(data = TRUE)},
% in each fold, we only update the list of \code{RPMs} by adjusting the
% diagonal elements to the vector of screening coefficients computed on
% the training data for the current fold, but do not modify the
% random elements of the projection matrices in each fold, to reduce the
% computational burden.
% Defaults  to \code{NULL}. If not provided, the values of the provided
% \code{RPM}s do not change.
% \end{itemize}
%
% For illustration, consider the built-in \code{rp_gaussian()} function.

% We turn to looking at the structure of the object:
% <<>>=
% unclass(obj)
% @
% The \code{generate_fun()} function returns a matrix with
% \code{m} rows and \code{length(included_vector)} columns.
% Note that \code{included_vector} gives the indices of the
% variables which have been selected by the screening
% procedure. In this case, the random projection does
% not use any data information and we are only interested in the
% length of this vector.

% The function \code{update_fun()} only converts the \class{family} object to a
% string and adds it as an attribute. Function \code{update_rpm_w_data()} is
% \code{NULL} as this random projection is data-agnostic.

%
\subsubsection{Marginal models}\label{sec:soft:sparmodel}
%
The \pkg{spareg} package uses the \class{sparmodel} class to define marginal
models for each ensemble element. These models assume a linear predictor
(a linear combination of projected variables).
Functions \code{spar_*()} create these objects and
have arguments \code{...} (i.e.,~the ellipsis) and
\code{control} (list used
in the main function for building the model). All arguments passed through the
ellipsis will be saved as attributes of the object.

The two functions implemented are
\begin{itemize}
\item \code{spar_glm()}--fits unregularized GLMs using \code{stats::glm()}
(default for Gaussian family with identity link).
\item \code{spar_glmnet()}--fits regularized GLMs using function
\code{glmnet::glmnet()}
(default for all other families using $\alpha=0$
with the minimal penalty considered in \pkg{glmnet}\footnote{
  In package \pkg{glmnet} $\lambda_\text{min}=\lambda_\text{max}/100$ where
  $\lambda_\text{max}$ is the smallest value that shrinks all coefficients
  to zero and is calculated as the maximum absolute inner product of the
  standardized predictors and the scores evaluated at the null model, divided
  by the number of observations times $\alpha$. For ridge regression
  $\alpha=0.001$ is used for computational purposes.}).
\end{itemize}
By default, \code{spar_glm()} is used for Gaussian family with identity link,
while \code{spar_glmnet()} is used for all other families.

Similar to the already presented \proglang{S}3 objects, an object of class `\code{sparmodel}' is a list containing
 an optional character \code{name} for printing, \code{generate_fun()}--a function used to generate the coefficients
and intercept value for each model in the ensemble, \code{update_fun()}--an (optional) function to modify the \class{sparmodel}
object at the start of the SPAR algorithm (e.g.,~adjusting the \class{family}
object for \pkg{glmnet} compatibility) and
\code{control}, the list of controls passed by the user in \code{spar_*()}.

%
An example for implementing a new marginal model class using the constructor
function \code{constructor\_sparmodel()} is presented in
Section~\ref{sec-ext-model}.
%
\subsection{Methods}
%
Methods \fct{print}, \fct{plot}, \fct{coef}, \fct{predict} are available for both `\code{spar}'
and `\code{spar.cv}' classes.
%
\subsubsection[print]{\fct{print}}
The \fct{print} method returns information on $M_\text{best}$ and $\nu_\text{best}$, the number of active predictors
(i.e.,~predictors which have at least a non-zero coefficient across the marginal models)
and a summary of the non-zero non-standardized coefficients. For `\code{spar.cv}'
it shows  information for the $(M, \nu)$
combination chosen by the best and the 1se rule.
<<first_example>>=
set.seed(12)
spar_res <- spar(example_data$x, example_data$y,
  xval = example_data$xtest, yval = example_data$ytest,
  nummods = c(5, 10, 15, 20, 25, 30))
spar_res
@
%
<<first_example_cv, echo = FALSE>>=
set.seed(12)
spar_cv <- spar.cv(example_data$x, example_data$y, nfolds = 3,
  nummods = c(5, 10, 15, 20, 25, 30))
@
%
\subsubsection[coef]{\fct{coef}}
Method \fct{coef} extracts coefficients  with arguments:
\begin{itemize}
\item \code{nu},\code{nummod}--threshold values and number of marginal models
to use in the calculation of the coefficients. Defaults to the values with
minimal validation error.
%
\item \code{aggregate}--aggregation method for coefficients  across the
\code{nummod} ensemble models. Options: \code{"mean"}
(default),  \code{"median"} or \code{"none"}(skips the aggregation and returns
the matrix of raw coefficients for each  marginal model--useful
for understanding the variability of the coefficients across the ensemble models).
\item \code{opt_par}--only for \class{spar.cv} with
\code{precompute_mode = "precompute_all"}. In this case, it chooses
between \code{"best"} or \code{"1se"}.   Otherwise, coefficients cannot be
generated without refitting on whole data with \code{"best"} or \code{"1se"}
parameters by calling \fct{spar} with the desired $(M, \nu)$ combination.
\end{itemize}
%
Method \fct{coef} returns an object of class \class{coefspar}, for which
user-friendly \fct{print} and \fct{summary} methods as well as extractor functions
\fct{get\_intercept} and \fct{get\_coef}
are provided.
%
<<eval=FALSE>>=
coef(spar_res)
@
<<echo=FALSE>>=
print(coef(spar_res), show = 5L)
@
<<>>=
get_intercept(coef(spar_res, nummod = 5, aggregate = "none"))
@
%
\subsubsection[predict]{\fct{predict}}
Functionality for computing predictions is provided through the method
\fct{predict} with arguments:
\begin{itemize}
\item \code{xnew}--new predictor matrix.
%
\item \code{type}--the type of predictions, either on
\code{"response"} level (default) or on \code{"link"} level.
%
\item  \code{avg_type}--stage of aggregation:  at the level of regression
coefficients (i.e.,~on the \code{link} level) or at the fitted values level
(i.e.,~\code{response} level).
% \begin{itemize}
%     \item \code{"link"} (default): combines models via an aggregation of the regression coefficients (i.e.,~on the \emph{link level}).
%     \item \code{"response"}: combines models via the fitted values by taking computing the predictions of each model in the ensemble and then aggregating on the \emph{response level}.
% \end{itemize}
%
\item \code{nu}, \code{nummod},  --	 threshold values and number of models
used to compute the averaged coefficients; defaults to optimal values.
%
\item \code{aggregate}--aggregation method (\code{"mean"} or \code{"median"}).
Defaults to mean aggregation.
%
\item \code{opt_par}--only for \class{spar.cv} objects fitted with
\code{precompute_mode = "precompute_all"}. In this case,it chooses
between \code{"best"} or \code{"1se"}.   Otherwise, predictions cannot be
generated without refitting on whole data with \code{"best"} or \code{"1se"}
parameters. To generate predictions, \fct{spar} can be called to refit the
ensembles with the desired $(M, \nu)$ combination.
\end{itemize}
%
% <<eval=FALSE>>=
% pred <- predict(spar_res, xnew = example_data$xtest)
% pred_cv <- predict(spar_cv, xnew = example_data$xtest, opt_par = "1se")
% @
%
\subsubsection[plot]{\fct{plot}}
%
Plotting functionality is provided through the \fct{plot} method with arguments:
\begin{itemize}
\item \code{plot\_type}--one of:
\begin{itemize}
\item \code{"val_measure"} plots the (cross-)validation \code{measure} for
either a grid of \code{nu} values for a fixed number of models \code{nummod} or
vice versa.
%
\item \code{"val_numactive"} plots the number of active variables for either a grid
of \code{nu} values for a fixed number of models \code{nummod} or vice versa.
%
\item \code{"res_vs_fitted"} produces a residuals-vs-fitted plot. The residuals
are computed as $r_i=y_i- \widehat y_i$, where $\widehat y_i$ is the prediction computed
on response level (available for \class{spar} objects and
\class{spar.cv} object with \code{precompute_mode = "precompute_all"}).
%
\item \code{"coefs"}  produces a plot of the standardized
coefficient values for each predictor in each marginal model (before thresholding).
For each predictor, the coefficients are sorted from largest to
smallest in absolute value across the marginal models and their values are
represented in the plot using a color scale  (available for \class{spar}
objects and \class{spar.cv} object with \code{precompute_mode = "precompute_all"}).
\end{itemize}
\item \code{plot_along}--for
\code{plot_type = "val_measure"} or
\code{"val_numactive"} only, it indicates whether to vary \code{nu} or \code{nummod} on the $x$~axis.
%
\item \code{nummod}--fixed value for number of models when
\code{plot_along = "nu"} for \code{plot_type = "val_measure"} or
\code{"val_numactive"};  if \code{plot_type = "res_vs_fitted"}, it is used in \code{predict()}.
%
\item \code{nu}--fixed value for  $\nu$ when \code{plot_along = "nummod"}
for \code{plot_type = "val_measure"} or \code{"val_numactive"}; if
\code{plot_type = "res_vs_fitted"}, it is used in  \code{predict()}.
\item \code{xfit} and \code{yfit}--predictor matrix and response vector for \code{plot_type = "res_vs_fitted"}.
%%
\item \code{prange}, \code{coef_order}--only for \code{plot_type = "coefs"}, optional customization parameters. Use \code{prange} to select a specific
range of predictors to display, and \code{coef_order} to rearrange their order.
%%
\item \code{opt_par}--only for \code{plot_type = "res_vs_fitted"} and
\class{spar.cv} objects fitted with
\code{precompute_mode = "precompute_all"}. It indicates whether
predictions should be computed using the \code{"best"} or \code{"1se"}
selection rule.
\end{itemize}
%
The \fct{plot} methods return objects of class \class{ggplot} \citep{ggplotR}.
The four types of plots described above for a \class{spar} object are
shown in Figure~\ref{fig:plot_example}.
%
% <<eval=FALSE>>=
% plot(spar_res)
% plot(spar_res, plot_type = "val_numactive")
% plot(spar_res,  plot_type = "coefs")
% plot(spar_res, plot_type = "res_vs_fitted", xfit = example_data$xtest,
%   yfit = example_data$ytest)
% @
%
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=\textwidth}
<<plotmethod, echo=FALSE, fig=TRUE, width=14, height=7,prefix.string=spareg>>=
library(ggplot2)
p1 <- plot(spar_res)+
  theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
p2 <- plot(spar_res, plot_type = "val_numactive") +
  theme(axis.text.x = element_text(angle = 45,vjust = 0.5))
p4 <- plot(spar_res, plot_type = "res_vs_fitted", xfit = example_data$xtest,
  yfit = example_data$ytest)
p3 <- plot(spar_res,  plot_type = "coefs")
p <- ggpubr::ggarrange(p1, p2, p3, p4,
  ncol = 2, nrow = 2)
p
@
\caption{\fct{plot} methods for \class{spar} object. They refer to the \code{plot\_type} of \code{"val\_measure"} (top left), \code{"val\_numactive"} (top right),  \code{"coefs"} (bottom left), and   \code{"res\_vs\_fitted"} (bottom right), respectively.
The red dots in the top
figures show the best $\nu$ chosen on the validation set for the optimal number
of models.
\label{fig:plot_example}}
\end{figure}
%
\subsection{Parallelization}
%
The package supports parallelization in the estimation of the
ensemble models via package \pkg{foreach} \citep{pkg:foreach}.
This is possible by setting the parameter \code{parallel = TRUE} in \code{spar()} or
\code{spar.cv()},
assuming that a parallel backend for \pkg{foreach} is registered
using the \code{registerDoParallel()}
function of package \pkg{doParallel} \citep{pkg:doParallel}.
To ensure reproducibility the \pkg{doRNG} \citep{pkg:doRNG}
can be used.
%
Our experiments suggest that parallelization
pays off especially for data sets with a larger number of observations~$n$.
A minimal example with a correlation-based
probabilistic screening and a Gaussian random projection matrix
is provided below:
<<parallel_example,message=FALSE,eval=FALSE>>=
ex4 <- simulate_spareg_data(n = 1000, p = 2000, ntest = 1000, seed = 123)
library(doParallel)
library(doRNG)
cl <- makeCluster(2, type = "PSOCK")
registerDoParallel(cl)
registerDoRNG(seed = 123)
spar_res_par <- spar(ex4$x, ex4$y, screencoef = screen_cor(),
  rp = rp_gaussian(), nummods = 50, parallel = TRUE)
stopCluster(cl)
@
%
\section{Extensibility} \label{sec-extensibility}

\subsection{Screening coefficients} \label{sec-extensscrcoef}
%
We demonstrate how screening coefficients implemented in package
\pkg{VariableScreening} can easily be incorporated in the framework of
\pkg{spareg}.

We start by defining the function for generating the screening coefficients
using function \code{screenIID()} in \pkg{VariableScreening}.
<<>>=
generate_scr_sirs <- function(object, x, y, ...) {
  ctrl <- object$control[names(object$control) %in%
                           names(formals(VariableScreening::screenIID))]
  res_screen <- do.call(function(...)
    VariableScreening::screenIID(x, y, ...), ctrl)
  res_screen$measurement
}
@
%
Note that \code{screenIID()} also takes \code{method} as an argument.
To allow for flexibility, we do not
fix the screening method in \code{generate_scr_sirs()} but rather allow the
user to pass a method through the \code{control} argument in the
\code{screen_*()} function. This function is created using the helper
\code{constructor_screencoef()}:
<<sirs_chunk>>=
screen_sirs <- constructor_screencoef(generate_fun = generate_scr_sirs,
  name = "screen_sirs")
@
%
We now call \code{spar()} with the newly created screening procedure.
We consider the method SIRS of \cite{zhu2011model}, which ranks the predictors
by their correlation with the rank-ordered response; we do not perform
probabilistic variable screening but employ the top $2n$ variables in
each marginal model.
The employed random projection is \code{rp_sparse()} where we
set $\psi=1/\sqrt{p}$, as proposed by \cite{LiHastie2006VerySparseRP}.
%
<<screen_sirs>>=
set.seed(123)
spar_example <- spar(example_data$x, example_data$y,
  screencoef = screen_sirs(type = "fixed", control = list(method = "SIRS")),
  rp = rp_sparse(psi = 1/sqrt(ncol(example_data$x))), measure = "mse")
@
%
\subsection{Random projections}\label{sec-ext-rp}
%
To extend the package, we incorporate the random projection method proposed by  \cite{cannings2017random},
which uses the Haar measure to generate orthogonal matrices.
This method involves drawing entries from a standard normal distribution and applying QR decomposition.
%
<<>>=
simulate_haar <- function(m, p) {
  R0 <- matrix(1/sqrt(p) * rnorm(p * m), nrow = p, ncol = m)
  t(qr.Q(qr(R0), complete = FALSE))
}
@
%
\cite{cannings2017random} propose selecting ``good'' random projections
to optimize out-of-sample prediction. Specifically, for each of the models in
the ensemble, $B$ Haar random projections are generated. The projection with
the lowest error on a test set (a random subset of proportion $\xi$ of the
data) is selected for projecting the variables. While the $B$ projections are
data-agnostic, the entire dataset is required to identify the best projection.

In the implementation, we
fit a GLM as a marginal model on the training data.
The performance is evaluated using deviance on the test set. Default values
can be specified in \code{update_fun()}, which must have arguments
\code{object}, \code{x} and \code{y} (predictor matrix and response vector
already standardized and supplied by `spar()`), \code{family} and  \code{...},
whereas all other relevant arguments from \code{spar()} will be passed
through the ellipsis.
% related to the family. Given that we will rely again on
% \code{glmnet::glmnet()} in each iteration~$b$,
% the family object is modified to a
% string containing the name of the family for families Gaussian, binomial and
% Poisson with canonical link (this will lead to using faster specialized
% algorithms). This modified family is
% added to the list of controls as \code{fit\_family}.
<<update_cannings>>=
update_rp_cannings <- function(object, x, y, family, ...) {
  if (is.null(object$control$alpha)) object$control$alpha <- 0
  if (is.null(object$control$B))  object$control$B <- 50
  if (is.null(object$control$xi)) object$control$xi <- 0.25
  if (is.null(object$control$family)) object$control$family <- family
  object
}
@
%
Next we specify the function for generating the random projection.
<<generate_cannings>>=
generate_cannings <- function(object, x, y, m, included_vector, ...) {
  xs <- x[, included_vector]; n <- nrow(x);  p <- ncol(xs)
  id_test <- sample(n, size = n * object$control$xi)
  xtrain <- xs[-id_test, ];  xtest <- xs[id_test,]
  ytrain <- y[-id_test];  ytest <- y[id_test]
  control_glm <-
    object$control[names(object$control) %in% names(formals(stats::glm.fit))]
  best_val <- Inf
  for (b in seq_len(object$control$B)) {
    RM <- simulate_haar(m, p)
    xrp <- tcrossprod(xtrain, RM)
    mod <- do.call(function(...)
      glm.fit(x = cbind(1, xrp), y = ytrain, ...), control_glm)
    eta_test <- cbind(1, tcrossprod(xtest, RM)) %*% coef(mod)
    pred <- object$control$family$linkinv(eta_test)
    out_dev <- sum(object$control$family$dev.resid(ytest, pred, wt = 1))
    if (out_dev < best_val) {
      best_val <- out_dev; best_RM <- RM
    }
  }
  return(best_RM)
}
@
Putting it all together, we get:
<<define_rp_cannings>>=
rp_cannings <- constructor_randomprojection(
  name = "RP Cannings and Samworth 2017",
  generate_fun = generate_cannings, update_fun = update_rp_cannings)
@
%
We can now estimate SPAR for a binomial model, after simulating
data from a linear
model and then transforming the response to a binary scale:
<<data_cannings>>=
set.seed(1234)
ex2 <- simulate_spareg_data(n = 200, p = 500, ntest = 100)
ys <- (ex2$y > 0) + 0; ysval <- (ex2$ytest > 0) + 0
@
We use only one value for the number of models $M=50$ \citep[which is in line
to recommendations in][]{cannings2017random}, no screening and no thresholding.
We set the maximum dimension for the random projection to $m_u = 20$
to avoid generating large matrices and for stability in the marginal binomial
GLMs.
<<rp_cannings_1,eval=TRUE>>=
set.seed(1234)
spar_cs <- spar(x = ex2$x, y = ys, family = binomial(),
  rp = rp_cannings(msup = 20), model = spar_glm(),
  nus = 0, nummods = 50, xval = ex2$xtest,  yval = ysval)
@
We can compare with results obtained by using the data-agnostic
sparse embedding matrix \code{rp_cw(data = FALSE)}:
<<rp_cannings_2, eval=TRUE>>=
set.seed(1234)
spar_cw <- spar(x = ex2$x, y = ys, family = binomial(),
  rp = rp_cw(data = FALSE, msup = 20), model = spar_glm(),
  nus = 0, nummods = 50, xval = ex2$xtest, yval = ysval)
@
The two approaches are compared by directly looking at the validation measure
achieved by the two ensembles on the provided validation set.
The extractor function \code{get_measure()} can be used to extract the
following information from the \class{spar} object:
<<compare>>=
rbind("cannings" = get_measure(spar_cs), "cw" = get_measure(spar_cw))
@
%
\subsection{Marginal models}\label{sec-ext-model}
%
Finally, we extend the package by implementing a robust GLM as a marginal model
using the package \pkg{robustbase} \citep{pkg:robustbase}.

We define a \code{generate_fun()} that uses \code{robustbase::lmrob()} for Gaussian
models with identity link and \code{robustbase::glmrob()}
for other family-link combinations.
The function must take by design the arguments \code{object} (of class \class{sparmodel}),
\code{x} and \code{y} (standardized predictors and response vector passed by \code{spar()}),
\code{z} (matrix of reduced predictors) as well as \code{...}.
It outputs a list with two elements: \code{gammas} (vector of coefficients)
and the \code{intercept} value.
<<>>=
model_glmrob <- function(object, x, y, z, ...) {
  requireNamespace("robustbase")
  fam <- object$control$family
  f <- ifelse(fam$family == "gaussian" & fam$link == "identity",
    robustbase::lmrob, robustbase::glmrob)
  glmrob_res <- do.call(function(...)
    f(y ~ as.matrix(z), ...), object$control)
  cfs <- coef(glmrob_res)
  list(gammas = cfs[-1], intercept = cfs[1])
}
@
%
Function \code{spar_glmrob()} that builds the `\code{sparmodel}'
object can be created by:
%
<<spar_glmrob>>=
spar_glmrob <- constructor_sparmodel(generate_fun = model_glmrob)
@
%
We demonstrate this on simulated Poisson data with 25\% contaminated predictors:
%
<<generate_example_data3>>=
set.seed(123)
example_data3 <- simulate_spareg_data(n = 100, p = 1000,
  ntest = 100, snr = 10, beta_vals = c(-1, 1)/10)
ypois <- round(exp(example_data3$y)) + 1; x <- example_data3$x;
np <- ncol(x) * nrow(x)
id_outliers_x <- sample(seq_len(np),  0.25 * np)
x[id_outliers_x] <- x[id_outliers_x] + 50
@
%
We then compare standard GLMs to SPAR with robust GLM marginal models,
using no screening,  \code{rp_gaussian(msup = 25)}---instead of the default
\code{msup} of $n/2=50$ to improve stability---and MAE as the selection
criterion.
<<run_glmrob_glm>>=
set.seed(1234)
spar_rob_res <- spar(x, ypois, family = poisson(),
  model = spar_glmrob(), rp = rp_gaussian(msup = 25), measure = "mae")
set.seed(1234)
spar_res <- spar(x, ypois, family = poisson(),
  model = spar_glm(), rp = rp_gaussian(msup = 25), measure = "mae")
@
%
The \code{get_model()} function extracts the optimal ensemble
from a \class{spar}  or \class{spar.cv} object.
It selects either the best model (\code{opt_par = "best"}, default)
or, for \class{spar.cv} objects,
the model based on the 1se rule (\code{opt_par = "1se"}).
<<compare_glmrob_glm>>=
best_rob <- get_model(spar_rob_res, opt_par = "best")
best_glm <- get_model(spar_res, opt_par = "best")
@
The validation measure for each of these models is given below.
The robust version achieves lower MAE, with all predictors retained in the
ensemble (each has at least one non-zero coefficient across models).
<<compare_glmrob_glm>>=
rbind("robglm" = get_measure(best_rob), "glm" = get_measure(best_glm))
@
%
\ifPositive

\section{Illustrations} \label{sec-illustrations}
%
\subsection{Face image data}
%
We illustrate \pkg{spareg} on the Isomap faces data set
containing $n = 698$
grayscale face images of size $p = 64 \times 64 = 4096$
of a single human face under varying poses and lighting conditions.
The face orientation measured by the horizontal viewing angles
is used as response variable.
The data set is publicly available for download at
\url{https://web.archive.org/web/20150922051706/http://isomap.stanford.edu/face_data.mat.Z}.
<<urls, echo=FALSE>>=
if (is_paper) {
  url1 <- "https://web.archive.org/web/20150922051706/"
  url2 <- "http://isomap.stanford.edu/face_data.mat.Z"
  ## On Ubuntu or MacOS use code below
  ## On Windows, download locally and unzip
  if (!file.exists(file.path(tmpdir, "face_data.mat"))) {
    # system('uncompress face_data.mat.Z')
    library("R.matlab")
    download.file(paste0(url1, url2),
                  file.path(tmpdir, "face_data.mat.Z"))
    system(sprintf('uncompress %s', file.path(tmpdir, "face_data.mat.Z")))
  }
}
@
<<unzip_file,echo=FALSE,message=FALSE>>=
if (is_paper) {
  library("R.matlab")
  facedata <- readMat(file.path(tmpdir,"face_data.mat"))
  x <- t(facedata$images)
  y <- facedata$poses[1,]
  x[, apply(x, 2, sd) < 0.01] <- 0
}
@
%
After downloading and decompressing the \code{.mat} file, we import it into \proglang{R}
using package \pkg{R.matlab} \citep{pkg:rmatlab}. Moreover, the data set
contains (nearly) constant columns. Functions \code{spar()} and \code{spar.cv()}
ignore constant columns from the estimation,
but we also set columns with a low standard deviation to zero. One image in
the data set is presented in the left panel of Figure~\ref{fig:faces_predictions}.
<<read_matlab,eval=FALSE>>=
library("R.matlab")
facedata <- readMat("face_data.mat")
x <- t(facedata$images); y <- facedata$poses[1, ]
x[, apply(x, 2, sd) < 0.01] <- 0
@
%
%
We split the data into 75\% training vs 25\%
test sample:
<<std_x,echo=FALSE>>=
if (is_paper) {
  set.seed(123)
  ntot <- length(y); ntest <- ntot * 0.25
  testind <- sample(ntot, ntest, replace = FALSE)
  xtrain <- as.matrix(x[-testind, ]); ytrain <- y[-testind]
  xtest <- as.matrix(x[testind, ]); ytest <- y[testind]
}
@
<<eval=FALSE>>=
set.seed(123)
testind <- sample(length(y), length(y) * 0.25, replace = FALSE)
xtrain <- as.matrix(x[-testind, ]); ytrain <- y[-testind]
xtest <- as.matrix(x[testind, ]); ytest <- y[testind]
@
We estimate the SPAR algorithm with cross-validation as in
\cite{parzer2024sparse} (with \code{lambda.min.ratio = 0.001}
to ensure convergence of the algorithm when computing the ridge screening
coefficients), and use a grid of $M=5,10,20,50$ models.
<<echo=FALSE>>=
if (is_paper) {
  file <- file.path(tmpdir, "faces_res.rda")
  if (!file.exists(file)) {
    set.seed(123)
    control_glmnet <- list(lambda.min.ratio = 0.001)
    library("spareg")
    Sys.time()
    spar_faces <- spar.cv(
      xtrain, ytrain,
      model = spar_glm(),
      screencoef = screen_glmnet(control = control_glmnet,
                                 reuse_in_rp = TRUE),
      rp = rp_cw(data = TRUE),
      nummods = c(5, 10, 20, 50),
      measure = "mse")
    Sys.time()
    save(spar_faces, file = file)
  } else {
    load(file)
  }
}
@
<<eval=FALSE>>=
set.seed(123)
control_glmnet <- list(lambda.min.ratio = 0.001)
spar_faces <- spar.cv(xtrain, ytrain,
  model = spar_glm(), rp = rp_cw(data = TRUE),
  screencoef = screen_glmnet(control = control_glmnet, reuse_in_rp = TRUE),
  nummods = c(5, 10, 20, 50), measure = "mse")
@
%\subsubsection{Results}
The \fct{print} method returns the following information about the fitted
\class{spar.cv} object. The combination $(M=20, \nu=0)$
 delivers the lowest MSE while the 1se rule chooses $(M=5, \nu=0.00424)$, with
 small differences in performance.
%
<<eval=FALSE>>=
spar_faces
@
<<echo=FALSE>>=
if (is_paper) spar_faces
@
%
The \code{plot} method for `\code{spar.cv}' objects displays by default the
measure employed in the cross-validation (in this case MSE)
for a grid of thresholding values $\nu$, where the number of models is fixed to
the value found to perform best in cross-validation
exercise
(Figure~\ref{fig:facesplot_valmeasure} top panel). Results for $M=5$, which
is the value chosen by the 1se rule, can be visualized by setting
\code{nummod = 5} (Figure~\ref{fig:facesplot_valmeasure} bottom panel).
%
<<plotfacesmeasure, eval=FALSE>>=
plot(spar_faces)
plot(spar_faces, nummod = 5)
@
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=0.95\textwidth}
<<plotfacesmeasurereal_best, echo=FALSE, fig=TRUE, width=8, height=3,prefix.string=spareg>>=
if (is_paper) {
  ## best
  p <- plot(spar_faces, digits = 3L) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                     hjust = 1))
  p
}
@
<<plotfacesmeasurereal_1se, echo=FALSE, fig=TRUE, width=8, height=3,prefix.string=spareg>>=
if (is_paper) {
  ## 1se
  p2 <- plot(spar_faces, nummod=5, digits = 3L) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                     hjust = 1))
  p2
}
@
\caption{Mean squared error vs. a grid of threshold values $\nu$ for fixed
number of models, once for the number of models chosen by the best CV measure
($M=20$, top panel) and one for the number of models chosen by the
1se rule ($M=5$, bottom panel) for the \emph{Isomap faces} data set.
Red points correspond to the best (minimal MSE) and sparsest ensemble model
within 1se of the best, respectively. Confidence bands represent one standard
deviation in the MSE across folds. \label{fig:facesplot_valmeasure}}
\end{figure}
%
The coefficients of the different variables (in this example pixels) obtained by
averaging over the coefficients of the marginal models can be extracted by \fct{coef}. For the ensemble chosen by the 1se rule we have:
%
<<facecoef1se, eval=FALSE>>=
face_coef_1se <- coef(spar_faces, opt_par = "1se")
summary(face_coef_1se)
@
<<facecoef1se,echo=FALSE>>=
if (is_paper) {
  face_coef_1se <- coef(spar_faces, opt_par = "1se")
  summary(face_coef_1se)
}
@
%
We visualize, in the right panel of Figure~\ref{fig:faces_predictions},
the effect of each pixel $\hat\beta^\text{1se}_j x^\text{new}_{i,j}$ in
predicting the face orientation in the  image given in the left panel.
We observe that the face contours and hairlines are most
predictive of face orientation.
%

Finally, we compare the performance of the estimated version of the SPAR algorithm
with another version proposed in \cite{Dunson2020TargRandProj}, where the sparse
random projection in \cite{ACHLIOPTAS2003JL} is combined with
probabilistic marginal correlation screening and no thresholding is performed.
We also use the values of $m_l=2\log p$ and $m_u=3/4 n$ as recommended in \cite{Dunson2020TargRandProj} to define the lower and upper bounds for the random projection dimension.
%
<<tarp_comparison_true,eval=FALSE>>=
set.seed(123)
m_l <- 2 * log(ncol(xtrain)); m_u <- 3/4 * nrow(xtrain)
tarp_faces <-  spar.cv(xtrain, ytrain, screencoef = screen_cor(),
  rp = rp_sparse(psi = 1/3, mslow = m_l, msup = m_u), model = spar_glm(),
  nus = 0, nummods = c(5, 10, 20, 50), measure = "mse")
@
<<tarp_comparison_true,echo=FALSE>>=
if(is_paper) {
  set.seed(123)
  m_l <- 2 * log(ncol(xtrain)); m_u <- 3/4 * nrow(xtrain)
  tarp_faces <-  spar.cv(xtrain, ytrain,
                         model = spar_glm(),
                         screencoef = screen_cor(),
                          nus = 0,
                         rp = rp_sparse(psi = 1/3, mslow = m_l, msup = m_u),
                         nummods = c(5, 10, 20, 50), measure = "mse")
  ynew_tarp <- predict(tarp_faces, xnew = xtest, opt_par = "1se")
}
@
%
The predictive performance on the \emph{test data} can be assessed by employing
the \code{predict()} method on the `\code{spar.cv}' object.
We will use the sparser solution chosen by the \code{opt_par = "1se"} rule and
the default mean aggregation for the predictions:
%
<<test_tarp_comparison_false,eval=FALSE>>=
ynew_spar <- predict(spar_faces, opt_par = "1se", xnew = xtest)
ynew_tarp <- predict(tarp_faces, opt_par = "1se", xnew = xtest)
c("SPAR" = mean((ytest - ynew_spar)^2),
  "TARP" = mean((ytest - ynew_tarp)^2))
@
<<test_tarp_comparison_true,echo=FALSE>>=
if(is_paper) {
ynew_spar <- predict(spar_faces, opt_par = "1se", xnew = xtest)
ynew_tarp <- predict(tarp_faces, opt_par = "1se", xnew = xtest)
  c("SPAR" = mean((ytest - ynew_spar)^2),
    "TARP" = mean((ytest - ynew_tarp)^2))}
@
%
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=0.45\textwidth}
<<plotfacesobs179, fig=TRUE,echo=FALSE,width=4.5,height=4.5,prefix.string=spareg>>=
if (is_paper) {
  i <- 179
  p <- ggplot(data.frame(X = rep(1:64,each=64),
                         Y = rep(64:1,64),
                         Z = facedata$images[,i]),
              aes(X, Y, fill = Z)) +
    geom_tile() +
    theme_void() +
    ggtitle(paste0("y = ", round(facedata$poses[1, i],1))) +
    theme(legend.position = "none",
          plot.title = element_text(hjust = 0.5))
  p
}
@
<<plotobsvspred179,fig=TRUE,echo=FALSE,width=4.5,height=4.5,prefix.string=spareg>>=
if (is_paper) {
  id <- 3
  p4 <- ggplot(data.frame(X = rep(1:64, each = 64),
                          Y = rep(64:1, 64),
                          effect = xtest[id,] * face_coef_1se$beta),
               aes(X, Y, fill = effect)) +
    geom_tile() +
    theme_void() +
    scale_fill_gradient2() +
    ggtitle(bquote(hat(y) == .(round(ynew_spar[id]))), ) +
    theme(plot.title = element_text(hjust = 0.5))
  p4
}
@
\caption{Left: Image corresponding to one observation in the
\emph{Isomap faces} data set. Right: The effect of each pixel
$\hat\beta^\text{1se}_j x^\text{new}_{i,j}$ in predicting the face
orientation in left panel. \label{fig:faces_predictions}}
\end{figure}
%
\subsection{Darwin data set}
%
The Darwin data set \citep{CILIA2022darwin} contains a binary response for
Alzheimer's disease (AD) together with extracted features from 25 handwriting
tests (18 features per task) for 89 AD patients and 85 healthy people ($n=174$).
The data set can be downloaded from
\url{https://archive.ics.uci.edu/dataset/732/darwin} as a \code{.zip} file.
The corresponding \code{.csv} file can then be imported into \proglang{R}.
<<readdarwinreal, echo=FALSE,message=FALSE>>=
if (is_paper) {
  # if (file.exists(file.path(tmpdir, "data.zip"))) {
  #  unzip(file.path(tmpdir, "data.zip")exdir = )
  #} else {
  if (!file.exists(file.path(tmpdir, "darwin.zip"))) {
    download.file("https://archive.ics.uci.edu/static/public/732/darwin.zip",
                  file.path(tmpdir, "darwin.zip"))
    unzip(zipfile = file.path(tmpdir, "darwin.zip"),
          files = "data.csv",
          exdir = tmpdir)
  }
  darwin_tmp <- read.csv(file.path(tmpdir, "data.csv"),
                         stringsAsFactors = TRUE)
}
@
<<readdarwinzip,eval=FALSE>>=
darwin_tmp <- read.csv("data.csv", stringsAsFactors = TRUE)
@
Before proceeding with the analysis, the data is screened for
multivariate outliers using the detect deviating cells (DDC) algorithm in package
\pkg{cellWise} \citep{rcellwise}.
<<darwinimpute,eval=FALSE>>=
darwin_orig <- list(y = as.numeric(darwin_tmp$class) - 1,
  x = darwin_tmp[, !(colnames(darwin_tmp) %in% c("ID", "class"))])
tmp <- cellWise::DDC(as.matrix(darwin_orig$x),
  list(returnBigXimp = TRUE, tolProb = 0.999, silent = TRUE))
darwin <- list(x = tmp$Ximp, y = darwin_orig$y)
@
<<echo=FALSE,message=FALSE>>=
if (is_paper) {
  darwin_orig <- list(
    x = darwin_tmp[, !(colnames(darwin_tmp) %in% c("ID", "class"))],
    y = as.numeric(darwin_tmp$class) - 1)
  tmp <- suppressMessages(cellWise::DDC(
    as.matrix(darwin_orig$x),
    list(returnBigXimp = TRUE, tolProb = 0.999, silent = TRUE)))
  darwin <- list(x = tmp$Ximp, y = darwin_orig$y)
}
@
% The structure of the data is:
% <<darwinstr,eval=FALSE>>=
% str(darwin)
% @
% <<echo=FALSE,message=FALSE>>=
% if (is_paper) str(darwin)
% @
%
We estimate the SPAR algorithm with the screening and random projection
introduced in \cite{parzer2024glms} for binomial family and logit link,
using one minus area under the ROC curve as the cross-validation measure:
<<spardarwin, eval=FALSE>>=
set.seed(1234)
spar_darwin <- spar.cv(darwin$x, darwin$y, family = binomial(logit),
  screencoef = screen_glmnet(reuse_in_rp = TRUE), nummods = c(10, 20, 30, 50),
  measure = "1-auc")
spar_darwin
@
<<spardarwinreal, echo=FALSE>>=
if (is_paper) {
  file <- file.path(tmpdir, "darwin_res.rda")
  if (!file.exists(file)) {
    set.seed(1234)
    spar_darwin <- spareg::spar.cv(darwin$x, darwin$y,
                                   family = binomial(logit),
                                   screencoef = screen_glmnet(reuse_in_rp = TRUE),
                                   nummods = c(10, 20, 30, 50),
                                   measure = "1-auc")
    save(spar_darwin, file = file)
  } else {
    load(file = file)
  }
  spar_darwin
}
@
No thresholding achieves the best measure and translates to almost all variables
being active (some variables can be inactive at $\nu=0$
as they may never be screened).
The 1se rule  achieves a much sparser model with a cross-validation measure
of around 0.075. The average number of active variables for the grid of thresholds
$\nu$ and $M = 20$ models is obtained by the following call and
shown in Figure~\ref{fig:darwin_activevars}.
<<plotdarwinact, eval=FALSE>>=
plot(spar_darwin, plot_type = "val_numactive")
@
%
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=0.8\textwidth}
<<plotdarwinactnumactive,fig=TRUE,echo=FALSE,width=8,height=3,prefix.string=spareg>>=
if (is_paper) {
  p <- plot(spar_darwin, plot_type = "val_numactive", digits=3L)
  p
}
@
\caption{Average number of active variables vs. threshold values
$\nu$ for $M=20$ models for the \emph{Darwin} data set. The red points
correspond to the average number of active variables for the ensemble with the
lowest cross-validation measure and to the one chosen by the 1se
rule. \label{fig:darwin_activevars}}
\end{figure}
%

Finally, standardized coefficients of the predictors across the maximum number of
considered marginal models (in this case $M=50$) can be visualized with
option \code{plot_type = "coefs"} in method \code{plot()}.
In this data set, the predictors are ordered by task, where the first 18
covariates represent different features measured for the first task, the
following 18 covariates represent the same features for the second task and
so on.
The predictors naturally group by feature, allowing for a more intuitive
assessment of how each feature (e.g., time spent on a task) relates to the
likelihood of Alzheimer's disease. Reordering the coefficients by feature--rather than by task--can be done by using the reordering argument
\code{coef\_order}. The resulting plot is shown in
Figure~\ref{fig:darwin_coefs}.
We add dashed vertical lines
to separate the features,  and color indicates coefficient sign (blue for positive, red for negative).
Consistent sign patterns across tasks suggest stable feature influences.
The presence of distinct blocks of blue or red lines reveals that the features
measured for different tasks generally exert a consistent--either positive or negative--influence on the probability of Alzheimer’s disease.
Line height shows selection frequency, with important features appearing in most models.
%%
\begin{figure}[t!]
\centering
\setkeys{Gin}{width=0.8\textwidth}
<<plotdarwincoef,fig=TRUE,echo=FALSE,width=7,height=4.5,prefix.string=spareg>>=
if (is_paper) {
  ntasks <- 25
  nfeat <- 18
  reorder_ind <- c(outer(
    (seq_len(ntasks) - 1) * nfeat,
    seq_len(nfeat), "+"))
  feat_names <- sapply(colnames(darwin$x)[seq_len(nfeat)],
                       function(name) substr(name, 1, nchar(name) - 1))
  p <- plot(spar_darwin,"coefs",coef_order = reorder_ind) +
    geom_vline(xintercept = 0.5 + seq_len(ntasks - 1) * ntasks,
               alpha = 0.2, linetype = 2) +
    annotate("text",x = (seq_len(nfeat) - 1) * ntasks + 12,
             y = 40,
             label=feat_names, angle = 90,
             size = 3)
  p
}
@
\caption{Coefficient plot for all variables and all considered $M=50$ models
in the SPAR ensemble for the \emph{Darwin} data set. The coefficients are
standardized, before thresholding. \label{fig:darwin_coefs}}
\end{figure}
%%
\fi

\section{Conclusion} \label{sec-conclusion}

Package \pkg{spareg} can be employed for modeling high-dimensional data in a GLM
framework, especially in settings where the number of predictors is much larger
than the number of observations. The package provides an implementation of a
computationally-efficient algorithm for sparse
projected and average regression (SPAR),
which combines variable screening and random projection in an ensemble of
GLMs to reduce dimensionality of the predictors.
Package \pkg{spareg} provides flexible classes for:
i)~specifying the
coefficient based on which screening should be performed (either in a classical
fashion, where the predictors with the largest screening coefficient
are selected for subsequent analysis, or in a probabilistic fashion, where
variables are sampled for inclusion with probabilities proportional
to their screening coefficients),
ii)~generating the random projection
to be employed in each marginal model, and
iii)~specifying the marginal models to be
used in the ensemble. Screening coefficients based on
marginal correlation between the predictors and the response, marginal
coefficients from a GLM, and ridge coefficients are provided in the package.
Moreover, several random projections are implemented:
the Gaussian and sparse matrices, which are data-agnostic and satisfy the JL
lemma, and the data-driven projection proposed in \cite{parzer2024sparse}
for linear
regression and extended to GLMs in \cite{parzer2024glms}.  This data-driven
approach, where information about the relationship among the responses and the
predictors is incorporated in the random projection by replacing the non-zero
elements of the random projection of \cite{Clarkson2013LowRankApprox} with the
ridge coefficients obtained with a small penalty, has the
advantage of ensuring that the true
regression coefficients can be recovered approximately after the projection.
Methodologically, the SPAR algorithm with ridge screening coefficients and
data-driven random projection  of \cite{parzer2024glms}
has been demonstrated to perform effectively in terms of predictive power and
variable ranking in settings with correlated predictors and across different
degrees of sparsity of the coefficient
vector, making it a suitable method for applications with high-dimensional correlated
predictors where the sparsity of the problem is unknown. This is the motivation
behind setting the random projection in \cite{parzer2024glms} as a default
in \pkg{spareg}.

Moreover, the flexibility and adaptability of the \pkg{spareg} package in
dealing with different types of screening, random projection and GLMs
make it an attractive choice for practitioners and researchers in a
wide variety of data settings.
It encourages exploration
of new methods for variable screening and random projections or the
combination of existing approaches to tailor solutions to specific data
requirements.

\section*{Computational details}

The results in this paper were obtained using \proglang{R}
\Sexpr{paste(R.Version()[6:7], collapse = ".")}.
\proglang{R} itself and all packages used are available from CRAN
at \url{https://CRAN.R-project.org/}.

\ifPositive
\else
<<session_info>>=
sessionInfo()
@
\fi

\ifPositive
\section*{Acknowledgments}

The authors are grateful to the two anonymous reviewers for their constructive
comments and suggestions, which significantly improved the quality of the
manuscript. We also extend our thanks to Jenny Fan for her careful review of
the manuscript and her insightful feedback and suggestions.
Roman Parzer and Laura Vana-G\"ur acknowledge funding from the Austrian
Science Fund (FWF) for the project ``High-dimensional statistical
learning: New methods to advance economic and sustainability policies''
(ZK 35), jointly carried out by WU Vienna University of Economics and
Business, Paris Lodron University Salzburg, TU Wien, and the Austrian
Institute of Economic Research (WIFO).

\fi


%% -- Bibliography -------------------------------------------------------------
%% - References need to be provided in a .bib BibTeX database.
%% - All references should be made with \cite, \citet, \citep, \citealp etc.
%%   (and never hard-coded). See the FAQ for details.
%% - JSS-specific markup (\proglang, \pkg, \code) should be used in the .bib.
%% - Titles in the .bib should be in title case.
%% - DOIs should be included where available.

\bibliography{spareg}
\ifPositive
\appendix

\section[Simulation Studies]{Simulation Study} \label{app:simulation-study}
\subsection[Computational efficiency and memory usage of spareg]{Computational efficiency and memory usage of \pkg{spareg}}
%%
To evaluate the computational efficiency and memory usage of the algorithm
implemented in \pkg{spareg}, we conducted a large-scale simulation study for the
linear regression setting comparing its performance against Elastic Net (ElNet)
under varying sample sizes ($n$) and numbers of predictors ($p$). The study
included 100 repetitions for each of 9 simulation settings. For investigating
the effect of increasing $p$ we cover combinations
of $n = 200$ and $p = 500, 1000, 2000, 5000, 10000, 20000$.
For investigating  the effect of increasing $n$ we cover combinations
$p = 2000$ and $n = 100, 200, 400, 800$. For each setting,
synthetic data from a linear regression model were generated using the
\code{simulate_spareg_data()} function with
default parameters: $100$ non-zero coefficients, a signal-to-noise ratio of 10,
and a compound symmetry correlation structure ($\rho = 0.5$).
We evaluated 7 methods:
%%
\begin{itemize}
\item SPAR (default): no screening, data-dependent random projection
  \code{rp\_cw(data = TRUE)},
\item SPAR CV (default): no screening,
  \code{rp_cw(data = TRUE)}, 10-fold cross-validation,
\item SPAR screen: screening with ridge
  coefficients \code{screen_glmnet(reuse_in_rp = TRUE)},
  \code{rp\_cw(data = TRUE)}, where the screening coefficients are used in
  constructing the CW random projection matrix),
\item SPAR screen CV (\code{screen\_glmnet(reuse_in_rp=TRUE)} and
  \code{rp\_cw(data = TRUE)}, where the ridge coefficients are used in
  constructing the CW random projection matrix,  and 10-fold cross-validation),
\item SPAR screen cor: \code{screen\_cor(reuse_in_rp=TRUE)} and \code{rp_cw(data = TRUE)}, where the correlation coefficients are used in the CW random projection matrix,
\item SPAR CW: no screening, data-agnostic \code{rp_cw(data = FALSE)},
\item SPAR Gaussian: no screening, \code{rp_gaussian()},
\item ElNet: Elastic Net with 10-fold CV for $\lambda$ and $\alpha=1/2$,
\item ElNet CValpha: Elastic Net with 10-fold CV for both  $\lambda$ and $\alpha$.
\end{itemize}
%%
Performance was assessed using the metrics:
%%
\begin{itemize}
\item computation time (elapsed time in seconds),
\item memory usage (in megabytes),
\item relative mean squared prediction error (rMSPE, relative to a model containing only
an intercept) on a test set
of $n_\text{test}=200$ observations,
\item to evaluate how well the methods are able to rank the variables, we employ the absolute value of the estimated coefficient vector to compute the partial area under the receiver operating characteristic curve (pAUC), where the number of false positives is limited to $n/2$.
\end{itemize}
%%
The results for increasing $n$ are presented in
Figure~\ref{fig:bench_sim_n} while Figure~\ref{fig:bench_sim_p} shows the results
for increasing $p$.
Increasing $n$ leads to the largest increase in computation time and memory usage for the SPAR CV methods.
 In terms of predictive performance, SPAR and its variants using the data-driven
 projections outperform on average ElNet for lower $n$, while
 elastic net is superior when $n$ increases.

When increasing $p$, we observe that SPAR and its variants maintain their computational efficiency and memory usage, while ElNetCV's scales worse.
Performance deteriorates in terms of both rMSPE and pAUC when increasing $p$
but SPAR and its data-driven projection variants do remain on average superior
to ElNet in terms of predictive performance and variable ranking.
%
\begin{figure}[t!]
\centering
\includegraphics[width=0.45\textwidth]{sim_bench/CompTime_n.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/Memory_n.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/rMSPE_n.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/pAUC_n.pdf}
\caption{Impact of increasing $n$ with fixed $p=2000$ for different methods in terms of
computation time in seconds (log10 $x$ and $y$~axis, top left panel), memory usage in megabytes (log10 $x$ and $y$~axis, top right panel),
rMSPE (bottom left panel), pAUC (bottom right panel).
\label{fig:bench_sim_n}}
\end{figure}
%
\begin{figure}[t!]
\centering
\includegraphics[width=0.45\textwidth]{sim_bench/CompTime_p.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/Memory_p.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/rMSPE_p.pdf}
\includegraphics[width=0.45\textwidth]{sim_bench/pAUC_p.pdf}
\caption{Impact of increasing $p$ with fixed $n=200$ for different methods in terms of
computation time in seconds (log10 $x$ and $y$~axis, top left panel), memory usage in megabytes (log10 $x$ and $y$~axis, top right panel),
rMSPE (bottom left panel), pAUC (bottom right panel).
\label{fig:bench_sim_p}}
\end{figure}

\fi
\end{document}