Towards Mechanised Proofs in
Double-Pushout Graph Transformation
Robert Söldner      Detlef Plump
Department of Computer Science, University of York, York, UK

R. Söldner and D. Plump
Towards Mechanised Proofs in Double-Pushout Graph Transformation

We formalise the basics of the double-pushout approach to graph transformation in the proof
assistant Isabelle/HOL and provide associated machine-checked proofs.
Specifically, we formalise graphs, graph morphisms and rules, and a definition of direct derivations
based on deletion and gluing. We then formalise graph pushouts and prove with Isabelle's help that
both deletions and gluings are pushouts. We also prove that pushouts are unique up to isomorphism.
The formalisation comprises around 2000 lines of source text. Our motivation is to pave the way for
rigorous, machine-checked proofs in the theory of the double-pushout approach, and to lay the
foundations for verifying graph transformation systems and rule-based graph programs by
interactive theorem proving.

[1]
authorJiří Adámek,
authorHorst Herrlich &
authorGeorge Strecker
(year1990): titleAbstract and Concrete

[2]
authorLuciano Baresi,
authorVahid Rafe,
authorAdel T. Rahmani &
authorPaola Spoletini
(year2008): titleAn Efficient Solution for
Model Checking Graph Transformation Systems.
journalElectronic Notes in Theoretical Computer
Science volume213(number1), pp.
pages3–21, 10.1016/j.entcs.2008.04.071.

[3]
authorYves Bertot &
authorPierre Castéran
(year2004): titleInteractive Theorem Proving and
Program Development — Coq'Art: The Calculus of Inductive Constructions.
seriesTexts in Theoretical Computer Science,
publisherSpringer, 10.1007/978-3-662-07964-5.

[4]
authorJon Haël Brenas,
authorRachid Echahed &
authorMartin Strecker
(year2018): titleVerifying graph transformations
with guarded logics.
In: booktitle2018 International Symposium on
Theoretical Aspects of Software Engineering (TASE),
organizationIEEE, pp. pages124–131,

[5]
authorGraham Campbell,
authorBrian Courtehoute &
authorDetlef Plump
(year2022): titleFast Rule-Based Graph
journalScience of Computer Programming
volume214, 10.1016/j.scico.2021.102727.

[6]
authorSimone André da Costa Cavalheiro, authorLuciana Foss & authorLeila Ribeiro
(year2017): titleTheorem proving graph grammars
with attributes and negative application conditions.
journalTheoretical computer science
volume686, pp. pages25–77,

[7]
authorMark Dowson (year1997):
titleThe Ariane 5 software failure.
journalACM SIGSOFT Software Engineering
Notes volume22(number2), p. pages84,

[8]
authorHartmut Ehrig
(year1979): titleIntroduction to the Algebraic
Theory of Graph Grammars.
In: booktitleProc. Graph-Grammars and Their
Application to Computer Science and Biology, seriesLecture Notes in Computer Science volume73,
publisherSpringer-Verlag, pp. pages1–69,

[9]
authorHartmut Ehrig,
authorKarsten Ehrig,
authorUlrike Prange &
authorGabriele Taentzer
(year2006): titleFundamentals of Algebraic Graph
seriesMonographs in Theoretical Computer Science,
publisherSpringer, 10.1007/3-540-31188-2.

[10]
authorAmir Hossein Ghamarian,
authorMaarten de Mol,
authorArend Rensink,
authorEduardo Zambon &
authorMaria Zimakova
(year2012): titleModelling and analysis using
journalInternational Journal on Software Tools
for Technology Transfer volume14(number1), pp.
pages15–40, 10.1007/s10009-011-0186-x.

[11]
authorAnnegret Habel,
authorJürgen Müller &
authorDetlef Plump
(year2001): titleDouble-Pushout Graph
Transformation Revisited.
journalMathematical Structures in Computer
Science volume11(number5), pp.
pages637–688, 10.17/S0960129501003425.

[12]
authorAnnegret Habel &
authorKarl-Heinz Pennemann
(year2009): titleCorrectness of high-level
transformation systems relative to nested conditions.
journalMathematical Structures in Computer
Science volume19(number2), pp.
pages245–296, 10.1017/S0960129508007202.

[13]
authorJ. Harrison (year2003):
titleFormal verification at Intel.
In: booktitleProc. 18th Annual IEEE Symposium
of Logic in Computer Science, publisherIEEE Computer Society,

[14]
authorKazuhiro Inaba,
authorSoichiro Hidaka,
authorZhenjiang Hu,
authorHiroyuki Kato &
authorKeisuke Nakano
(year2011): titleGraph-Transformation
Verification Using Monadic Second-Order Logic.
In: booktitleProc. of the 13th International
ACM SIGPLAN Symposium on Principles and Practices of Declarative Programming
(PPDP '11), publisherACM Press, p. pages17–28,

[15]
authorN.G. Leveson &
authorC.S. Turner
(year1993): titleAn investigation of the
Therac-25 accidents.
volume26(number7), pp. pages18–41,

[16]
authorTobias Nipkow &
authorGerwin Klein
(year2014): titleConcrete Semantics — With
publisherSpringer, 10.1007/978-3-319-10542-0.

[17]
authorLars Noschinski
(year2015): titleA graph library for
journalMathematics in Computer Science
volume9(number1), pp. pages23–39.

[18]
authorLawrence C. Paulson,
authorTobias Nipkow &
authorMakarius Wenzel
(year2019): titleFrom LCF to
journalFormal Aspects of Computing
volume31(number6), pp. pages675–698,

[19]
authorChristopher M. Poskitt &
authorDetlef Plump
(year2014): titleVerifying Monadic Second-Order
Properties of Graph Programs.
In: booktitleProc. International Conference on
Graph Transformation (ICGT 2014), seriesLecture Notes
in Computer Science volume8571, publisherSpringer,
pp. pages33–48, 10.1007/978-3-319-09108-2_3.

[20]
authorArend Rensink
(year2008): titleExplicit State Model Checking
for Graph Grammars.
In: booktitleConcurrency, Graphs and Models,
Essays Dedicated to Ugo Montanari on the Occasion of His 65th Birthday,
seriesLecture Notes in Computer Science
volume5065, publisherSpringer, pp.
pages114–132, 10.1007/978-3-540-68679-8_8.

[21]
authorArend Rensink,
authorÁkos Schmidt &
authorDániel Varró
(year2004): titleModel checking graph
transformations: A comparison of two approaches.
In: booktitleProc. International Conference on
Graph Transformation (ICGT 2004), organizationSpringer, pp.
pages226–241, 10.1007/978-3-540-30203-2_17.

[22]
authorRobert Söldner &
authorDetlef Plump
(year2022): titleTowards Mechanised Proofs in
Double-Pushout Graph Transformation.
In: booktitleProc. International Workshop on
Graph Computation Models (GCM 2022).

[23]
authorMartin Strecker
(year2018): titleInteractive and automated
proofs for graph transformations.
journalMathematical Structures in Computer
Science volume28(number8), pp.
pages1333–1362, 10.1017/S096012951800021X.

[24]
authorJan Stückrath
(year2016): titleVerification of
Well-Structured Graph Transformation Systems.
Ph.D. thesis, schoolUniversität Duisburg-Essen.

[25]
authorDániel Varró
(year2004): titleAutomated formal verification
of visual modeling languages by model checking.
journalSoftware & Systems Modeling
volume3(number2), pp. pages85–113,

[26]
authorMakarius Wenzel:
titleThe Isabelle/Isar Reference Manual.

[27]
authorMarkus Wenzel
(year1999): titleIsar — A Generic
Interpretative Approach to Readable Formal Proof Documents.
In: booktitleTheorem Proving in Higher
Order Logics, seriesLecture Notes in Computer
Science, publisherSpringer, pp. pages167–183.

[28]
authorGia S. Wulandari &
authorDetlef Plump
(year2021): titleVerifying Graph Programs with
Monadic Second-Order Logic.
In: booktitleProc. International Conference on
Graph Transformation (ICGT 2021), seriesLecture Notes
in Computer Science volume12741, publisherSpringer,
pp. pages240–261, 10.1007/978-3-030-78946-6_13.


#  Subset selection for linear mixed models

Daniel R. Kowal  Dobelman Family Assistant Professor, Department of
Statistics, Rice University, Houston, TX (daniel.kowal@rice.edu).

###### Abstract

Linear mixed models (LMMs) are instrumental for regression analysis with
structured dependence, such as grouped, clustered, or multilevel data.
However, selection among the covariates—while accounting for this structured
dependence—remains a challenge. We introduce a Bayesian decision analysis for
subset selection with LMMs. Using a Mahalanobis loss function that
incorporates the structured dependence, we derive optimal linear coefficients
for (i) any given subset of variables and (ii) all subsets of variables that
satisfy a cardinality constraint. Crucially, these estimates inherit shrinkage
or regularization and uncertainty quantification from the underlying Bayesian
model, and apply for any well-specified Bayesian LMM. More broadly, our
decision analysis strategy deemphasizes the role of a single “best” subset,
which is often unstable and limited in its information content, and instead
favors a collection of near-optimal subsets. This collection is summarized by
key member subsets and variable-specific importance metrics. Customized subset
search and out-of-sample approximation algorithms are provided for more
scalable computing. These tools are applied to simulated data and a
longitudinal physical activity dataset, and demonstrate excellent prediction,
estimation, and selection ability.

Keywords: Bayesian analysis; hierarchical models; prediction; regression;
variable selection

##  1 Introduction

Linear mixed models (LMMs) enable regression analysis in the presence of
structured dependence, such as longitudinal data, grouped or clustered
observations, or spatio-temporal effects. LMMs are widespread in both Bayesian
and classical statistical analysis and include many hierarchical models and
linear regression as special cases. We consider LMMs of the general form

$\bm{y}=\bm{X}\bm{\beta}+\bm{Z}\bm{u}+\bm{\epsilon},$ (1)

where $\bm{y}$ is the $N$-dimensional response, $\bm{X}$ is the $N\times p$
matrix of covariates, $\bm{\beta}$ is the $p$-dimensional vector of fixed
effects regression coefficients, $\bm{Z}$ is the $N\times q$ random effects
design matrix, $\bm{u}$ is the $q$-dimensional vector of random effects
regression coefficients, and $\bm{\epsilon}$ is the $N$-dimensional
observation error. Model (1) is paired with the assumptions that $\bm{u}$ and
$\bm{\epsilon}$ are uncorrelated and mean zero with
$\mbox{Cov}(\bm{u})=\bm{\Sigma_{u}}$ and
$\mbox{Cov}(\bm{\epsilon})=\bm{\Sigma_{\epsilon}}$. Most commonly, the random
effects $\bm{u}$ and the errors $\bm{\epsilon}$ are endowed with Gaussian
distributions, but our approach does not require any specific distributional
assumptions beyond these moments.

The benefit of the LMM (1) is that it marries the classical linear regression
term $\bm{X}\bm{\beta}$ with a random effects term $\bm{Z}\bm{u}$ to capture
structural dependence unexplained by $\bm{X}\bm{\beta}$. More formally, (1)
can be expressed in the marginal form $\bm{y}=\bm{X}\bm{\beta}+\bm{\nu},$
where $\bm{\nu}\coloneqq\bm{Z}\bm{u}+\bm{\epsilon}$ has mean zero and
covariance $\bm{Z}\bm{\Sigma_{u}}\bm{Z}^{\prime}+\bm{\Sigma_{\epsilon}}$. The
covariance of $\bm{\nu}$ incorporates elements of the random effects design
$\bm{Z}$, the random effects covariance $\bm{\Sigma_{u}}$, and the observation
error covariance $\bm{\Sigma_{\epsilon}}$. LMMs are capable of modeling a
broad variety of dependence structures; specific examples are given in Section
2.1.

Regardless of the structured dependence in the LMM, a core goal of regression
analysis is _selection_ among the $p$ (fixed effects) covariates $\bm{x}$.
Selection provides interpretable summaries of the data, reduced storage
requirement, and often better prediction and lower estimation variability. We
emphasize four main priorities that motivate our approach:

  1. (P1)

The selection criteria and accompanying performance metrics should account for
the _structured dependence_ modeled by the LMM;

  2. (P2)

Selection should be applied _jointly_ across covariates rather than
_marginally_ for each covariate;

  3. (P3)

Selection of a _single_ “best” subset of covariates should be accompanied by
an analysis of “near-optimal” subsets of covariates; and

  4. (P4)

The inference and selection procedure should be computationally scalable in
$N$ and $p$.

P1 simply states that any structured dependence worth modeling in the LMM must
also be included in the selection and evaluation process—which renders many
existing tools ineligible. P2 notes that variables selected using marginal
criteria, such as hypothesis tests of the form $H_{0j}:\beta_{j}=0$ or
posterior inclusion probabilities from sparse Bayesian models, do not
necessarily satisfy any joint optimality criteria. Hence, reporting the
marginally-selected variables as a joint subset of variables often lacks
justification. More directly, P2 is satisfied only for subset selection. Yet
subset selection is accompanied by other challenges, including selection
instability and computational scalability. P3 addresses the instability of
subset selection: the “best” subset often changes dramatically under minor
perturbations or resampling of the data. This effect is most pronounced in the
presence of correlated covariates, weak signals, or small sample sizes, and
undermines the elevated status of a “best” subset. By instead collecting
“near-optimal” subsets, we acquire more information about the competing
(predictive) explanations. Lastly, P4 recognizes the computational burdens of
subset search and demands tools that are feasible for moderate to large $N$
and $p$.

Variable selection for LMMs has most commonly relied on penalized maximum
likelihood estimation. Foster et al., (2007) and Wang et al., (2011)
incorporated random effects within an adaptive lasso estimation procedure to
account for genetic and experimental effects in quantitative trait loci
analysis and plant population studies, respectively. Bondell et al., (2010)
and Ibrahim et al., (2011) selected fixed and random effects jointly using a
modified Cholesky decomposition with adaptive lasso or SCAD penalties. These
Cholesky parametrizations are order-dependent, so permutations of the columns
of $\bm{Z}$ can produce different estimates and selections. Müller et al.,
(2013) also noted that the accompanying algorithms can be slow and fail to
converge, and reviewed alternative strategies such as information criteria.
Fan and Li, (2012) selected fixed effects by marginalizing over the random
effects and maximizing a penalized (marginal) log-likelihood. The primary
limitation is the need for a “proxy matrix” for the inverse marginal
covariance (of $\bm{\nu}$); Fan and Li, (2012) simply used a multiple of the
identity matrix, but this ignores the random effects covariance structure. In
general, such penalized estimators can address priorities P1, P2, and P4, but
not P3: they focus on selecting a single “best” subset, and the accompanying
(forward) search paths are too restrictive to enumerate a sufficiently rich
collection of competitive subsets.

From a Bayesian perspective, Chen and Dunson, (2003) and Kinney and Dunson,
(2007) proposed sparsity-inducing spike-and-slab priors for both the fixed and
random effects. These priors are compatible with our approach. The primary
distinction is the mechanism for selection: Chen and Dunson, (2003) and Kinney
and Dunson, (2007) compute posterior probabilities for all possible submodels.
However, this strategy is computationally prohibitive and unreliable for small
to moderate $p+q$, since only a small fraction of possible subsets can be
visited regularly within the stochastic search Gibbs sampler. Hence, P4 is not
satisfied. Marginal criteria such as posterior inclusion probabilities or
hard-thresholding resolve these challenges, but fail to satisfy P2.

More broadly, Lindley, (1968) and Hahn and Carvalho, (2015) have argued that
selection is a _decision problem_ distinct from model specification. Sparsity
or shrinkage priors cannot alone select variables: the prior is a component of
the Bayesian model while the selection process requires its own criteria,
typically a loss function that balances accuracy with sparsity. This decision
analysis approach to selection has proven useful for functional regression
(Kowal and Bourgeois,, 2020), seemingly unrelated regressions (Puelz et al.,,
2017), and graphical models (Bashir et al.,, 2019), among others. However,
these methods were not designed for LMMs and therefore fail to satisfy P1. In
addition, with the exception of Kowal, 2021a , these decision analysis
approaches use (variations of) $\ell_{1}$-penalties and suffer from the same
restrictive search paths as in classical penalized regression, which fails to
satisfy P3.

We propose a Bayesian approach for subset search and selection in LMMs that
satisfies P1–P4. Using decision analysis with a predictive loss function that
directly incorporates the structured dependence in (1), we derive and compute
the optimal linear coefficients for (i) any _given_ subset
$\mathcal{S}\subseteq\\{1,\ldots,p\\}$ of variables and (ii) all subsets of
variables that satisfy a _cardinality constraint_ $|\mathcal{S}|\leq k$ (P1,
P2). These optimal coefficients are computable for any Bayesian LMM and
inherit model-based regularization and posterior predictive uncertainty
quantification. Linear coefficients are compared across subsets using out-of-
sample predictive performance metrics that leverage both the structural
dependencies and the predictive uncertainty from the Bayesian LMM. From these
metrics, we construct the _acceptable family_ of near-optimal subsets, which
collects those subsets that perform nearly as well as the “best” subset with
nonnegligble probability under the Bayesian LMM (P3). The acceptable family is
more informative and robust than the “best” subset—which itself is a
member—and is summarized using other key member subsets and variable
importance metrics. Customized subset search and out-of-sample approximation
algorithms are provided to enable scalable computing (P4).

We focus on subset selection of _fixed effects_ covariates, but note that the
distinction between fixed and random effects is less pertinent for Bayesian
modeling. Unlike frequentist LMMs that place a prior only on the random
effects, Bayesian models require a prior on all parameters. Here, we consider
“fixed effects” as those covariates designated for selection, while “random
effects” capture the structured dependencies unmodeled by the fixed effects.

The methodology is applied to moderate-to-vigorous physical activity (MVPA)
data from the 2005-2006 National Health and Nutrition Examination Survey
(NHANES). Repeated measurements of daily MVPA were recorded for each subject
for one to seven days, along with several subject-specific demographic,
health, and behavioral variables. The goal is to analyze which of these
variables predict MVPA while adhering to priorities P1–P4 and accounting for
the structured dependence implied by the longitudinal observations.

The paper is outlined as follows: Section 2 develops the methodology and
algorithms; Section 3 provides results for simulated data; Section 4 presents
an application to physical activity data; Section 5 concludes. Supporting
information includes a document with additional simulation results, and
additional results from the NHANES application, proofs of all results, and
computational details; and R code to reproduce the simulation study and data
analysis. An R package is available at
https://github.com/drkowal/BayesSubsets.

##  2 Methods

###  2.1 Predictive decision analysis for linear mixed models

Bayesian analysis of LMMs pairs the model (1) with suitable priors on
$\bm{\beta}$ and $\bm{u}$ and a distributional choice for $\bm{\epsilon}$ to
determine the likelihood, which is typically Gaussian. Specific choices will
depend on the formulation of (1) and are discussed subsequently; for now, we
denote a generic Bayesian LMM by $\mathcal{M}$. The Bayesian model
$\mathcal{M}$ induces a data-generating process via the posterior predictive
distribution,

$p_{\mathcal{M}}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\mid\bm{y}\\}=\int
p_{\mathcal{M}}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\mid\bm{\theta}\\}\
p_{\mathcal{M}}(\bm{\theta}\mid\bm{y})\ d\bm{\theta},$ (2)

where $\bm{\theta}$ denotes the model $\mathcal{M}$ parameters including
$\bm{\beta}$, $\bm{u}$, and any covariance parameters. The terms in the
integrand are defined by the likelihood in (1) evaluated at the covariate
values $\bm{\tilde{X}}$ and $\bm{\tilde{Z}}$ and the joint posterior
distribution under $\mathcal{M}$. Informally, (2) describes the distribution
of future or unobserved data $\bm{\tilde{y}}$ at the design matrices
$\bm{\tilde{X}}$ and $\bm{\tilde{Z}}$ conditional on the observed data
$\bm{y}$ and according to model $\mathcal{M}$. The choice of $\bm{\tilde{X}}$
and $\bm{\tilde{Z}}$ can target covariate values or subpopulations of interest
and determines the type of predictive observations, such as predictions for a
new group or new measurements on an existing group. Absent other
considerations, our default is the observed matrices, $\bm{\tilde{X}}=\bm{X}$
and $\bm{\tilde{Z}}=\bm{Z}$

While the posterior predictive distribution formalizes the model-based
uncertainty about unobserved data
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$, _predictive decision
analysis_ determines the actions—point or interval predictions or estimators,
selection among hypotheses, etc.—that provide optimal data-driven decision-
making under $\mathcal{M}$. Here, the goals are to (i) compute optimal linear
coefficients for any given subset $\mathcal{S}\subseteq\\{1,\ldots,p\\}$ of
variables, (ii) conduct an efficient search over candidates subsets, and (iii)
evaluate and compare predictive performance among subsets—all while adhering
to the priorities P1–P4. Predictive decision analysis requires a loss function
of the form
$\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}\\}$,
which enumerates the cost of an action $\bm{\delta}$ when
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$ is realized. In accordance
with P1 and P2, we deploy a Mahalanobis loss function

$\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{\mathcal{S}};\bm{\psi}\\}=\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}\|_{\bm{\Omega}_{\bm{\psi}}}^{2}$
(3)

where $\bm{\delta}_{\mathcal{S}}$ is the $p$-dimensional linear coefficients
with zeros for any index $j\not\in\mathcal{S}$ and the norm
$\|\bm{v}\|_{\bm{\Omega}_{\bm{\psi}}}^{2}=\bm{v}^{\prime}\bm{\Omega}_{\bm{\psi}}\bm{v}$
depends on a positive definite weighting matrix $\bm{\Omega}_{\bm{\psi}}$ that
can depend on model parameters $\bm{\psi}$.

For LMMs, a natural choice of $\bm{\Omega_{\psi}}$ is the inverse marginal
covariance of $\bm{\nu}$,

$\bm{\Omega}_{\bm{\psi}}=(\bm{\tilde{Z}}\bm{\Sigma_{u}}\bm{\tilde{Z}}^{\prime}+\bm{\Sigma_{\epsilon}})^{-1}$
(4)

with $\bm{\psi}=(\bm{\Sigma_{u}},\bm{\Sigma_{\epsilon}})$. While the central
quantity
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}$
in (3) explicitly measures the linear predictive ability of a subset of
variables $\mathcal{S}$, the choice of (4) incorporates weighting to account
for the structured dependencies that are unknown yet modeled by the random
effects under the LMM. With (4), the Mahalanobis loss (3) resembles a
multivariate Gaussian (negative) log-likelihood. However, this mathematical
similarity should not be confused with a distributional assumption: the
Mahalanobis predictive loss (3) inherits a joint posterior predictive
distribution
$p_{\mathcal{M}}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\psi}\mid\bm{y}\\}$
under $\mathcal{M}$.

For any _given_ subset $\mathcal{S}$, the optimal coefficients are obtained by
minimizing the posterior expected loss under $\mathcal{M}$:

$\bm{\hat{\delta}}_{\mathcal{S}}\coloneqq\arg\min_{\bm{\delta}_{\mathcal{S}}}\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{\mathcal{S}};\bm{\psi}\\},$
(5)

which averages over the joint uncertainty in
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$ and $\bm{\psi}$ conditional on
the data $\bm{y}$ and according to the model $\mathcal{M}$. The solution to
(5) is derived explicitly:

######  Lemma 1.

When
$\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\|_{\bm{\Omega}_{\bm{\psi}}}^{2}<\infty$,
the optimal coefficients in (5) for a _given_ subset
$\mathcal{S}\subseteq\\{1,\ldots,p\\}$ is given by the nonzero entries

$\bm{\hat{\delta}}_{\mathcal{S}}=(\bm{\tilde{X}}_{\mathcal{S}}^{\prime}\bm{\hat{\Omega}}\bm{\tilde{X}}_{\mathcal{S}})^{-1}\bm{\tilde{X}}_{\mathcal{S}}^{\prime}\bm{\hat{y}^{\Omega}}$
(6)

with zeros for indices $j\not\in\mathcal{S}$, where
$\bm{\tilde{X}}_{\mathcal{S}}$ subsets the columns of $\bm{\tilde{X}}$ based
on $\mathcal{S}$ and
$\bm{\hat{\Omega}}\coloneqq\mathbb{E}_{[\bm{\psi}\mid\bm{y}]}\bm{\Omega}_{\bm{\psi}}$
and
$\bm{\hat{y}^{\Omega}}\coloneqq\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\\{\bm{\Omega}_{\bm{\psi}}\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\\}$
are posterior expectations under $\mathcal{M}$.

A generalized inverse may be substituted when the solution (6) is nonunique.

Lemma 1 explicitly derives the optimal Bayesian estimator under Mahalanobis
loss for any _given_ subset $\mathcal{S}$. The optimal
$\bm{\hat{\delta}}_{\mathcal{S}}$ is a “fit to the fit” from $\mathcal{M}$,
and therefore inherits shrinkage or regularization from the Bayesian LMM. For
illustration, consider a fixed and known weighting matrix $\bm{\Omega}$: the
pseudo-response variable is $\bm{\hat{y}^{\Omega}}=\bm{\Omega}\bm{\hat{y}}$
where
$\bm{\hat{y}}\coloneqq\mathbb{E}_{[\bm{\tilde{y}}\mid\bm{y}]}\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})=\bm{\tilde{X}}\bm{\hat{\bm{\beta}}}+\bm{\tilde{Z}}\bm{\hat{u}}$
for $\bm{\hat{\beta}}\coloneqq\mathbb{E}_{[\bm{\beta}\mid\bm{y}]}\bm{\beta}$
and $\bm{\hat{u}}\coloneqq\mathbb{E}_{[\bm{u}\mid\bm{y}]}\bm{u}$. The
regularization from $\mathcal{M}$—usually applied via the priors for
$\bm{\beta}$ and $\bm{u}$—is valuable for point prediction and estimation, and
its absence in classical subset selection is detrimental (Hastie et al.,,
2020).

The optimal coefficients in (6) resemble generalized least squares (GLS)
estimators for linear regression, including LMMs. The primary challenge in GLS
estimation is that the inverse covariance or weight matrix
$\bm{\Omega_{\psi}}$ is unknown. Feasible GLS iteratively estimates the
covariance and the linear coefficients via plug-in estimation, which is
suboptimal. For LMMs, Fan and Li, (2012) substituted a multiple of the
identity matrix for $\bm{\Sigma_{u}}$ in (4) in order to avoid estimation of
this covariance. These concessions are avoided in our approach: we solve a GLS
optimization problem, but compute model-based expectations jointly over the
unknown parameters—including the necessary inverse covariance matrix. The
estimate of $\bm{\Omega_{\psi}}$ derives from the Bayesian LMM (1), which can
benefit from the model-based regularization induced by the choice of shrinkage
or sparsity priors under $\mathcal{M}$.

###  2.2 The Mahalanobis weight matrix

To illustrate the use of the weighting matrix $\bm{\Omega_{\psi}}$, we
consider several examples. Since
$\bm{\Omega}_{\bm{\psi}}=\bm{\Sigma_{\epsilon}}^{-1}-\bm{\Sigma_{\epsilon}}^{-1}\bm{\tilde{Z}}(\bm{\Sigma_{u}}^{-1}+\bm{\tilde{Z}}^{\prime}\bm{\Sigma_{\epsilon}}^{-1}\bm{\tilde{Z}})^{-1}\bm{\tilde{Z}}^{\prime}\bm{\Sigma_{\epsilon}}^{-1}$
by the Woodbury identity, the common assumption of
$\bm{\Sigma_{\epsilon}}=\sigma_{\epsilon}^{2}\bm{I}_{N}$ results in the
simplification

$\bm{\Omega}_{\bm{\psi}}=\sigma_{\epsilon}^{-2}(\bm{I}_{N}-\bm{\tilde{Z}}\bm{\Sigma_{u^{*}}}^{-1}\bm{\tilde{Z}}^{\prime}),$
(7)

where
$\bm{\Sigma_{u^{*}}}\coloneqq\sigma_{\epsilon}^{2}\bm{\Sigma_{u}}^{-1}+\bm{\tilde{Z}}^{\prime}\bm{\tilde{Z}}$.
The Mahalanobis predictive loss (3) then decomposes as

$\sigma_{\epsilon}^{2}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{\mathcal{S}};\bm{\psi}\\}=\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}\|_{2}^{2}-\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}\|_{\bm{\tilde{Z}}\bm{\Sigma_{u^{*}}}^{-1}\bm{\tilde{Z}}^{\prime}}^{2}$

which isolates the contribution from the squared error loss and the
Mahalanobis loss based only on $\bm{\tilde{Z}}$ and
$\bm{\Sigma_{u^{*}}}$—i.e., the critical terms in the random effects
component.

The optimal coefficients in (6) require computation of $\bm{\hat{\Omega}}$ and
$\bm{\hat{y}^{\Omega}}$ under $\mathcal{M}$. We further consider two important
examples: the random intercept model (Section 2.2.1) and the random slope
model (Section 2.2.2).

####  2.2.1 Random intercept model

Consider longitudinal observations $\\{y_{ij}\\}_{j=1}^{m_{i}}$ on each
subject $i=1,\ldots,n$, so $N=\sum_{i=1}^{n}m_{i}$. The within-subject
correlations are often modeled using the random intercept model

$y_{ij}=\bm{x}_{i}^{\prime}\bm{\beta}+u_{i}+\epsilon_{ij},$ (8)

usually with $u_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{u}^{2})$
and $\epsilon_{ij}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma_{\epsilon}^{2})$. The crucial role of $u_{i}$ cannot
be ignored: since
$\mbox{Corr}(y_{ij},y_{ij^{\prime}}\mid\bm{x}_{i},\bm{\beta})=\sigma_{u}^{2}/(\sigma_{u}^{2}+\sigma_{\epsilon}^{2})$,
$\sigma_{u}$ accounts for the within-subject correlation that remains
unexplained by the covariates $\bm{x}_{i}$. Model (8) is a special case of (1)
with $\bm{\Sigma_{u}}=\sigma_{u}^{2}\bm{I}_{N}$ and
$\bm{Z}=\mbox{bdiag}\\{\bm{1}_{m_{i}}\\}_{i=1}^{n}$ is a block diagonal matrix
with $n$ $m_{i}$-dimensional vectors of ones.

For predictive decision analysis, let $\bm{\tilde{x}}_{i}$ denote the target
covariate values and $\tilde{m}_{i}$ the number of observations for each
subject $i=1,\ldots,\tilde{n}$, which determines $\bm{\tilde{Z}}$. The
subject-specific predictive variables are
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})=(\bm{\tilde{y}}_{1}^{\prime},\ldots,\bm{\tilde{y}}_{\tilde{n}}^{\prime})^{\prime}$
with
$\bm{\tilde{y}}_{i}=(\tilde{y}_{i1},\ldots,\tilde{y}_{i\tilde{m}_{i}})^{\prime}$
and the fixed effects covariate matrix is
$\bm{\tilde{X}}=(\bm{1}_{m_{1}}^{\prime}\otimes\bm{\tilde{x}}_{1},\ldots,\bm{1}_{m_{\tilde{n}}}^{\prime}\otimes\bm{\tilde{x}}_{\tilde{n}})^{\prime}$.
To compute $\bm{\Omega_{\psi}}$, observe that
$\bm{\Sigma_{u^{*}}}=\mbox{bdiag}\\{\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\bm{1}_{\tilde{m}_{i}}^{\prime}\bm{1}_{\tilde{m}_{i}}\\}_{i=1}^{n}=\mbox{diag}\\{\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i}\\}_{i=1}^{n}$
and
$\bm{\Sigma_{u^{*}}}^{-1}=\mbox{diag}\\{(\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i})^{-1}\\}_{i=1}^{n}$,
so the Mahalanobis weight matrix (up to $\sigma_{\epsilon}^{2}$) is

$\sigma_{\epsilon}^{2}\bm{\Omega_{\psi}}=\mbox{bdiag}\Big{\\{}\bm{I}_{\tilde{m}_{i}}-\frac{1}{\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i}}\bm{1}_{\tilde{m}_{i}}\bm{1}_{\tilde{m}_{i}}^{\prime}\Big{\\}}_{i=1}^{n}$
(9)

and does not require any numerical matrix inversions. Given (9), the
Mahalanobis predictive loss simplifies to

$\sigma_{\epsilon}^{2}\|\bm{\tilde{y}}-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}\|_{\bm{\Omega}_{\bm{\psi}}}^{2}=\sum_{i=1}^{\tilde{n}}\Big{[}\sum_{j=1}^{\tilde{m}_{i}}({\tilde{y}}_{ij}-\bm{\tilde{x}}_{i}^{\prime}\bm{\delta}_{\mathcal{S}})^{2}-\frac{1}{\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i}}\Big{\\{}\sum_{j=1}^{\tilde{m}_{i}}(\tilde{y}_{ij}-\bm{\tilde{x}}_{i}^{\prime}\bm{\delta}_{\mathcal{S}})\Big{\\}}^{2}\Big{]}$
(10)

which clearly isolates the difference between the Mahalanobis loss and squared
error loss. In particular, (10) incorporates the sign of the errors
$e_{ij}\coloneqq{\tilde{y}}_{ij}-\bm{\tilde{x}}_{i}^{\prime}\bm{\delta}_{\mathcal{S}}$.
For example, suppose $\tilde{n}=1$ and $\tilde{m}_{1}=2$, so the Mahalanobis
loss (up to $\sigma_{\epsilon}^{2}$) is
$e_{1}^{2}+e_{2}^{2}-(\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+2)^{-1}(e_{1}+e_{2})^{2}$.
The squared error loss $e_{1}^{2}+e_{2}^{2}$ is invariant to the signs of the
errors. However, the second term in (10) includes a reduction in the loss by a
factor of $(e_{1}+e_{2})^{2}$, which is larger when the errors have the same
sign. Compared to the squared error loss, this Mahalanobis loss is more
forgiving for errors in the same direction—and this is accentuated when
$\sigma_{u}$ is large—which reflects the within-subject correlation induced by
the underlying model (8).

The posterior expectation $\bm{\hat{\Omega}}$ of (9) is straightforward to
compute, for example given posterior samples of
$\\{\sigma_{\epsilon}^{2},\sigma_{u}^{2}\\}$. To compute the posterior
expectation of
$\bm{\Omega_{\psi}}\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$, the block
diagonality simplifies this term to $n$ blocks of the form
$\sigma_{\epsilon}^{-2}\\{\bm{I}_{\tilde{m}_{i}}-(\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i})^{-1}\bm{1}_{\tilde{m}_{i}}\bm{1}_{\tilde{m}_{i}}^{\prime}\\}\bm{\tilde{y}}_{i}=\sigma_{\epsilon}^{-2}\bm{\tilde{y}}_{i}-\\{\sigma_{\epsilon}^{-2}(\sigma_{\epsilon}^{2}/\sigma_{u}^{2}+\tilde{m}_{i})^{-1}\sum_{j=1}^{\tilde{m}_{i}}\tilde{y}_{ij}\\}\bm{1}_{\tilde{m}_{i}}.$
The posterior expectation of each $\tilde{m}_{i}$-dimensional vector is easily
computable given posterior samples of
$\\{\sigma_{\epsilon}^{2},\sigma_{u}^{2},\bm{\tilde{y}}_{i}\\}_{i=1}^{\tilde{n}}$.

######  Remark.

These simplifications also provide a scalable Gibbs sampling algorithm for a
Gaussian random intercept model with large $N,p$. We apply a _joint_ sampling
step for all fixed and random effects that (nearly) maintains the
computational scalability of Bayesian linear regression _without_ the random
intercepts. For simplicity, fix $m_{i}=m$ and let $\bm{Y}=[y_{ij}]$ denote the
$n\times m$ matrix of observations. The strategy is to decompose
$[\bm{\beta},\\{u_{i}\\}_{i=1}^{n}\mid\bm{y},-]=[\bm{\beta}\mid\bm{y},-][\\{u_{i}\\}_{i=1}^{n}\mid\bm{y},\bm{\beta},-]$
and draw from the constituents of the product. Under the prior $\bm{\beta}\sim
N(\bm{0},\bm{\Sigma}_{\beta})$ and $u_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma_{u}^{2})$, the regression coefficients satisfy
$[\bm{\beta}\mid\bm{y},-]\sim
N(\bm{Q}_{\beta}^{-1}\bm{\ell}_{\beta},\bm{Q}_{\beta}^{-1})$ with
$\bm{Q}_{\beta}=\omega^{tot}\bm{X}^{\prime}\bm{X}+\bm{\Sigma}_{\beta}^{-1}$
and $\bm{\ell}_{\beta}=\bm{X}^{\prime}\bm{Y}{\bm{\omega}_{\cdot}}^{tot}$,
where $\omega^{tot}=\sum_{k,j}\bm{\Omega}_{kj}$ is the grand sum and
${\bm{\omega}_{\cdot}}^{tot}=(\omega_{\cdot 1}^{tot},\ldots,\omega_{\cdot
m}^{tot})^{\prime}$ for $\omega_{\cdot
j}^{tot}=\sum_{k=1}^{m}\bm{\Omega}_{kj}$ is the column sums of $\bm{\Omega}$
defined in (9). Notably, this distributional form matches the canonical
posterior distribution of the regression coefficients for (non-LMM) Gaussian
linear regression, which admits efficient sampling methods for large $n,p$
(Bhattacharya et al.,, 2016; Nishimura and Suchard,, 2018). The random
intercepts are sampled independently via $[u_{i}\mid\bm{y},\bm{\beta},-]\sim
N(Q_{u_{i}}^{-1}\ell_{u_{i}},Q_{u_{i}}^{-1})$ with
$Q_{u_{i}}=m\sigma_{\epsilon}^{-2}+\sigma_{u}^{-2}$ and
$\ell_{u_{i}}=\sigma_{\epsilon}^{-2}\sum_{j=1}^{m}(y_{ij}-\bm{x}_{i}^{\prime}\bm{\beta})$
for $i=1,\ldots,n$. Most important, these sampling steps for high-dimensional
Bayesian random intercept regression are comparable to those for high-
dimensional Bayesian linear regression, and only add minimal additional
computations related to summations of $\bm{\Omega}$ and (parallelizable) draws
of the scalar random intercepts $u_{i}$. The remaining sampling steps for the
variance components are standard but depend on the choice of priors. These
results also apply to Gibbs samplers for Gaussian mixture models (e.g.,
Dirichlet process mixtures of Gaussians) for $\epsilon_{ij}$ and/or $u_{i}$.

####  2.2.2 Random slope model

Subject-specific slopes are common in hierarchical or multilevel models. By
applying (1) with $\bm{Z}=\mbox{bdiag}\\{\bm{x}_{i}^{\prime}\\}_{i=1}^{n}$,
the random slope model allows for subject-specific deviations from the
population-level coefficients $\bm{\beta}$ (including a subject-specific
intercept):

$y_{i}=\bm{x}_{i}^{\prime}\bm{\beta}_{i}+\epsilon_{i},\quad\bm{\beta}_{i}\coloneqq\bm{\beta}+\bm{u}_{i}.$
(11)

Model (11) is often accompanied by shrinkage priors on $\bm{\beta}$ and
$\bm{u}_{i}$ to regularize against unnecessary predictors and unnecessary
heterogeneity, respectively. Predictive decision analysis with Mahalanobis
loss enables coefficient estimation and subset selection for $\bm{x}$ (see
Section 2.4) while adjusting for the heterogeneities induced by the random
effects $\bm{u}_{i}$.

When $\bm{\Sigma_{\epsilon}}=\sigma_{\epsilon}^{2}\bm{I}_{N}$, the key term
$\bm{\Sigma_{u^{*}}}^{-1}$ in the inverse covariance (7) is directly available
from the Sherman-Morrison formula,
$\bm{\Sigma_{u^{*}}}^{-1}=\sigma_{\epsilon}^{-2}\mbox{bdiag}\big{\\{}\bm{\Sigma_{u_{i}}}-\bm{\Sigma_{u_{i}}}\bm{\tilde{x}}_{i}\bm{\tilde{x}}_{i}^{\prime}\bm{\Sigma_{u_{i}}}/(\sigma_{\epsilon}^{2}+\bm{\tilde{x}}_{i}^{\prime}\bm{\Sigma_{u_{i}}\bm{\tilde{x}}_{i}})\big{\\}}_{i=1}^{n}$.
The accompanying Mahalanobis weight matrix then simplifies to the diagonal
matrix $\bm{\Omega_{\psi}}=\mbox{diag}\\{\omega_{i}\\}_{i=1}^{n}$ with
$\omega_{i}\coloneqq
1/(\sigma_{\epsilon}^{2}+\bm{\tilde{x}}_{i}^{\prime}\bm{\Sigma_{u_{i}}}\bm{\tilde{x}}_{i})$,
which is computable without numerical matrix inversions. The implied
Mahalanobis predictive loss is the weighted least squares
$\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})-\bm{\tilde{X}}\bm{\delta}_{\mathcal{S}}\|_{\bm{\Omega}_{\bm{\psi}}}^{2}=\sum_{i=1}^{\tilde{n}}\omega_{i}\\{\tilde{y}_{i}(\bm{\tilde{x}}_{i})-\bm{\tilde{x}}_{i}^{\prime}\bm{\delta}_{\mathcal{S}}\\}^{2}$.
The subject-specific weights $\omega_{i}$ are primarily driven by
$\bm{\tilde{x}}_{i}^{\prime}\bm{\Sigma_{u_{i}}}\bm{\tilde{x}}_{i}$, where
$\bm{\Sigma_{u_{i}}}$ is the covariance of the subject-specific deviations
$\bm{u}_{i}=\bm{\beta}_{i}-\bm{\beta}$. The posterior expectations required by
Lemma 1 are straightforward:
$\bm{\Omega_{\psi}}=\mbox{diag}\\{\omega_{i}\\}_{i=1}^{n}$ and
$\bm{\Omega_{\psi}}\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$ is an
$\tilde{n}$-dimensional vector with elements
$\\{\omega_{i}\tilde{y}_{i}(\bm{\tilde{x}}_{i})\\}_{i=1}^{\tilde{n}}$, both of
which are easily computable given posterior samples of
$\\{\sigma_{\epsilon}^{2},\bm{\Sigma_{u_{i}}},\tilde{y}_{i}(\bm{\tilde{x}}_{i})\\}_{i=1}^{\tilde{n}}$.

###  2.3 Subset search for linear mixed models

Although Lemma 1 produces the optimal linear coefficients for a _given_ subset
$\mathcal{S}$, it does not guide the subset _search_ or _selection_ process.
To remedy this, we append the Mahalanobis loss function (3) with a
_cardinality constraint_ and define the optimal action

$\bm{\hat{\delta}}_{k}\coloneqq\arg\min_{\bm{\delta}}\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta};\bm{\psi}\\}\quad\mbox{subject
to}\quad\|\bm{\delta}\|_{0}\leq k$ (12)

so $\bm{\hat{\delta}}_{k}$ provides the optimal coefficients among all subsets
with _at most_ $k$ variables. The solution in (12) resemble the “best subset
selection” problem in classical regression (e.g., Miller,, 1984), suitably
modified for Bayesian decision analysis.

The cardinality constraint diverges from the ubiquitous strategy among
decision analysis methods for variable selection, which is to append the loss
function (e.g., (3)) with an $\ell_{1}$-penalty to encourage sparsity among
the coefficients (Hahn and Carvalho,, 2015). Such a strategy may be viewed as
a convex relaxation of (12). However, the $\ell_{1}$-penalty introduces
additional regularization—beyond the regularization from $\mathcal{M}$—and can
overshrink true signals. Adaptive lasso-type adjustments are available (Kowal
et al.,, 2021) but cannot circumvent this issue entirely. Further, the
(adaptive) lasso-based search paths are highly constrained within the space of
all possible subsets, and therefore cannot enumerate a sufficiently broad
collection of competitive subsets to satisfy P3.

We instead target (12) directly, and provide a substantial simplification of
the solution:

######  Theorem 1.

The optimal coefficients (12) using the loss (3) and the cardinality
constraint $\|\bm{\delta}\|_{0}\leq k$ (with $k\leq p$) are

$\bm{\hat{\delta}}_{k}=\arg\min_{\bm{\delta}}\|\bm{y}^{*}-\bm{X}^{*}\bm{\delta}\|_{2}^{2}\quad\mbox{subject
to}\quad\|\bm{\delta}\|_{0}=k$ (13)

where $\bm{y}^{*}\coloneqq\bm{\hat{\Omega}}^{-1/2}\bm{\hat{y}^{\Omega}}$,
$\bm{X}^{*}\coloneqq\bm{\hat{\Omega}}^{1/2}\bm{\tilde{X}}$, and
$(\bm{\hat{\Omega}}^{1/2})^{\prime}\bm{\hat{\Omega}}^{1/2}=\bm{\hat{\Omega}}$.

The expected predictive Mahalanobis loss in (12) is reduced to a squared error
loss involving pseudo-data $\bm{y}^{*}$ and $\bm{X}^{*}$. Most important, the
squared error representation in (13) enables application of state-of-the-art
subset search algorithms for _classical_ linear regression (Furnival and
Wilson,, 2000; Bertsimas et al.,, 2016) to the setting of (12). The pseudo-
data $\bm{y}^{*}$ and $\bm{X}^{*}$ are a one-time computing cost, while the
matrix square root $\bm{\hat{\Omega}}^{1/2}$ often admits fast Cholesky
decompositions (e.g., block diagonality in Section 2.2.1) or direct
computations (e.g., diagonality in Section 2.2.2) depending on the form of the
LMM (1). In addition, Theorem 1 reduces the search space from $2^{k}$ subsets
to $p\choose{k}$ subsets. For any subset of size $k$, we simply apply Lemma 1
to compute the optimal linear coefficients as in (6).

Despite these advantageous results, Theorem 1 also highlights the limitations
of the representation in (12). First, this solution does not consider P3:
there may be many near-optimal subsets of smaller sizes, yet all subsets with
$|\mathcal{S}|<k$ are immediately discarded. Second, this solution does not
favor parsimony: the optimal coefficients $\bm{\hat{\delta}}_{k}$ are the
_largest_ allowable subset under the cardinality constraint. Hence, optimizing
over all possible subsets is achieved by setting $k=p$, which yields a trivial
solution:

######  Corollary 1.

The optimal coefficients under the loss (3) and computed _across all possible
subsets_ $\mathcal{S}\subseteq\\{1,\ldots,p\\}$ are
$\bm{\hat{\delta}}_{\mathcal{\widehat{S}}}\coloneqq\arg\min_{\mathcal{S},\bm{\delta}}\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta};\bm{\psi}\\}=\bm{\hat{\delta}}_{\\{1,\ldots,p\\}}=(\bm{\tilde{X}}^{\prime}\bm{\hat{\Omega}}\bm{\tilde{X}})^{-1}\bm{\tilde{X}}^{\prime}\bm{\hat{y}^{\Omega}}$
with $\mathcal{\widehat{S}}=\\{1,\ldots,p\\}$.

Clearly, selection via direct optimization is inadvisable: the selected subset
includes all variables and therefore is invariant to the data or the model.

In conjunction, Theorem 1 and Corollary 1 imply the need to (i) search over
multiple cardinalities $k=1,\ldots,p$ and (ii) develop alternative metrics to
compare subsets of distinct sizes. Even with the simplifications of Theorem 1
and the accompanying subset search algorithms (Furnival and Wilson,, 2000;
Bertsimas et al.,, 2016), it is often necessary to restrict the search space
when $p$ is moderate or large. We do so by bounding the maximum subset size
$s_{max}\leq p$ and the number of subsets $s_{k}\leq{p\choose k}$ of each size
$k$.

First, we pre-screen to select the $s_{max}=\min\\{p,35\\}$ covariates that
have the largest effect sizes under the LMM. Although this pre-screening
applies a _marginal_ criterion, it is based on a _joint_ model under
$\mathcal{M}$. In that sense, this procedure is similar to the most popular
Bayesian variable selection strategies based on posterior inclusion
probabilities or hard-thresholding. In our case, this is a coarse pre-
screening technique, not a terminal selection procedure.

Second, we apply the _branch-and-bound algorithm_ (BBA; Furnival and Wilson,,
2000) to filter to the “best” $s_{k}$ subsets of each size
$k=1,\ldots,s_{max}$. BBA searches through a tree-based enumeration of all
possible subsets (up to size $s_{max}$), yet avoids an exhaustive subset
search by carefully eliminating non-competitive subsets (or branches)
according to least squares. Hence, application of BBA requires a least squares
representation of the expected predictive Mahalanobis loss in (12), which we
provide below:

######  Lemma 2.

Let $\bm{\delta}_{1}$ and $\bm{\delta}_{2}$ denote linear coefficients. When
$\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\|\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\|_{\bm{\Omega}_{\bm{\psi}}}^{2}<\infty$,
we have the ordering equivalence
$\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{1};\bm{\psi}\\}\leq\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{2};\bm{\psi}\\}$
if and only if
$\|\bm{y}^{*}-\bm{X}^{*}\bm{\delta}_{1}\|_{2}^{2}\leq\|\bm{y}^{*}-\bm{X}^{*}\bm{\delta}_{2}\|_{2}^{2}.$

The key implication of Lemma 2 is that we may directly apply BBA using the
pseudo-data $\bm{y}^{*}$ and $\bm{X}^{*}$ (defined in Theorem 1) to obtain the
“best” $s_{k}$ subsets of each size $k=1,\ldots,s_{max}$. Ideally, $s_{k}$
should be set to the largest size possible for a given computing budget. We
use the default values $s_{k}=15$ or $s_{k}=100$ and apply the efficient BBA
implementation in the leaps package in R. However, Lemma 2 also enables any
other subset search strategy based on least squares (e.g., Bertsimas et al.,,
2016).

###  2.4 Acceptable families for near-optimal subsets

Subset selection via the decision analysis in (12) is incomplete: the solution
returns only the “best” model of each size $k$ and trivially prefers the
largest possible subset. Additional tools are needed to (i) compare subsets of
distinct sizes and (ii) collect the _near-optimal_ subsets in accordance with
P3. For these tasks, we use _out-of-sample_ predictive performance and adapt
the _acceptable family_ of Kowal, 2021b  for the LMM setting. Informally, the
acceptable family is the collection of all subsets that (nearly) match the
predictive performance of the “best” subset with nonnegligble probability
under $\mathcal{M}$. By studying this collection of near-optimal subsets, we
deemphasize the role of a single “best” subset in favor of many distinct yet
predictively-competitive alternatives. The acceptable family has been applied
for Bayesian subset selection (Kowal, 2021a, ), $\ell_{1}$-penalized selection
(Kowal et al.,, 2021), and targeted variable selection (Kowal, 2021b, ), but
none have considered LMMs.

The acceptable family is built by evaluating out-of-sample predictive
performance, which requires careful consideration for LMMs. For repeated or
longitudinal observations, it must be determined whether to evaluate
predictions for new subjects or for new measurements on existing subjects. For
concreteness, we proceed under the longitudinal setting of Section 2.2.1 and
evaluate predictions on new subjects. Modifications for other cases are
available.

Consider $n$ subjects with $m_{i}$ observations per subject, $i=1,\ldots,n$.
We implement a Bayesian $K$-fold cross-validation procedure, where the $K$
folds are taken across subjects $i=1,\ldots,n$. Let
$\mathcal{I}_{k}\subset\\{1,\ldots,n\\}$ denote the $k$th validation set,
where each subject point appears in one validation set,
$\cup_{k=1}^{K}\mathcal{I}_{k}=\\{1,\ldots,n\\}$. By default, we use $K=10$
validation sets that are equally-sized, mutually exclusive, and selected
randomly from $\\{1,\ldots,n\\}$. For each subset $\mathcal{S}$, we define the
out-of-sample _empirical loss_

$\mathcal{L}_{\mathcal{S}}\coloneqq\frac{1}{K}\sum_{k=1}^{K}\mathcal{L}_{\mathcal{S}}(k),\quad\mathcal{L}_{\mathcal{S}}(k)\coloneqq\frac{1}{|\mathcal{I}_{k}|}\mathcal{L}(\bm{y}_{\mathcal{I}_{k}},\bm{\hat{\delta}}_{\mathcal{S}}^{-\mathcal{I}_{k}};\bm{\hat{\psi}}^{-\mathcal{I}_{k}}),$
(14)

where
$\bm{y}_{\mathcal{I}_{k}}\coloneqq\\{\bm{y}_{i}\\}_{i\in\mathcal{I}_{k}}$
denotes the response variables on the validation data with
$\bm{y}_{i}=(y_{i1},\ldots,y_{im_{i}})^{\prime}$,
$\bm{\hat{\delta}}_{\mathcal{S}}^{-\mathcal{I}_{k}}\coloneqq\arg\min_{\bm{\delta}_{\mathcal{S}}}\mathbb{E}_{[\bm{\tilde{y}},\bm{\psi}\mid\bm{y}_{-\mathcal{I}_{k}}]}\mathcal{L}(\bm{\tilde{y}}_{\mathcal{I}_{k}},\bm{\delta}_{\mathcal{S}};\bm{\psi}\\}$
are the optimal coefficients (5) but estimated using only the training data
$\bm{y}_{-\mathcal{I}_{k}}\coloneqq\\{\bm{y}_{i}\\}_{i\not\in\mathcal{I}_{k}}$,
and, with abuse of notation, $\bm{\hat{\psi}}^{-\mathcal{I}_{k}}$ in (14)
indicates the Mahalanobis loss (3) with weighting matrix
$\bm{\hat{\Omega}}^{-\mathcal{I}_{k}}\coloneqq\mathbb{E}_{[\bm{\psi}\mid\bm{y}_{-\mathcal{I}_{k}}]}\bm{\Omega}_{\bm{\psi}}$.
The empirical loss (14) resembles classical $K$-fold cross-validation with a
point estimate for each Mahalanobis loss weighting matrix. From this quantity,
we define the “best” subset for out-of-sample point prediction,

$\mathcal{S}_{min}\coloneqq\arg\min_{\mathcal{S}}\mathcal{L}_{\mathcal{S}},$
(15)

so that $\bm{\hat{\delta}}_{\mathcal{S}_{min}}$ are the optimal linear
coefficients for the subset $\mathcal{S}_{min}$ that minimizes (14).

To define the acceptable family, we first introduce the out-of-sample
_predictive loss_ analogous to (14):

$\widetilde{\mathcal{L}}_{\mathcal{S}}\coloneqq\frac{1}{K}\sum_{k=1}^{K}\widetilde{\mathcal{L}}_{\mathcal{S}}(k),\quad\widetilde{\mathcal{L}}_{\mathcal{S}}(k)\coloneqq\frac{1}{|\mathcal{I}_{k}|}\mathcal{L}(\bm{\tilde{y}}_{\mathcal{I}_{k}}^{-\mathcal{I}_{k}},\bm{\hat{\delta}}_{\mathcal{S}}^{-\mathcal{I}_{k}};\bm{\psi}^{-\mathcal{I}_{k}})$
(16)

where $\bm{\tilde{y}}_{\mathcal{I}_{k}}^{-\mathcal{I}_{k}}\sim
p_{\mathcal{M}}[\\{\bm{\tilde{y}}(\bm{x}_{i},\bm{z}_{i})\\}_{i\in\mathcal{I}_{k}}\mid\bm{y}_{-\mathcal{I}_{k}}]$
denotes the predictive variables in the validation set conditional on the
training data and $\bm{\psi}^{-\mathcal{I}_{k}}$ similarly conditions only on
the training data. Unlike the empirical loss ${\mathcal{L}}_{\mathcal{S}}$,
the predictive loss $\widetilde{\mathcal{L}}_{\mathcal{S}}$ incorporates out-
of-sample predictive _uncertainty_ under $\mathcal{M}$, as well as the
uncertainty regarding relevant model parameters $\bm{\psi}$. The uncertainty
reflects the fact that the validation data $\bm{y}_{\mathcal{I}_{k}}$
represent only one possible realization of observables at the covariate values
$\\{\bm{x}_{i},\bm{z}_{i}\\}_{i\in\mathcal{I}_{k}}$. The (out-of-sample)
predictive distribution provides alternative model-based realizations, and
hence is informative for quantifying the uncertainty of out-of-sample
predictive performance.

Using the predictive loss, the acceptable family is defined as those subsets
that are “near-optimal” relative to $\mathcal{S}_{min}$:

$\mathbb{A}_{\eta,\varepsilon}\coloneqq\big{\\{}\mathcal{S}:\mathbb{P}_{\mathcal{M}}\big{(}\widetilde{\mathcal{D}}_{\mathcal{S}_{min},\mathcal{S}}<\eta\big{)}\geq\varepsilon\big{\\}},\quad\eta\geq
0,\varepsilon\in[0,1]$ (17)

where $\widetilde{\mathcal{D}}_{\mathcal{S}_{min},\mathcal{S}}\coloneqq
100\times(\widetilde{\mathcal{L}}_{\mathcal{S}}-\widetilde{\mathcal{L}}_{\mathcal{S}_{min}})/\widetilde{\mathcal{L}}_{\mathcal{S}_{min}}$
is the percent increase in predictive loss for subset $\mathcal{S}$ relative
to $\mathcal{S}_{min}$, $\eta\geq 0\%$ is the margin, and
$\varepsilon\in[0,1]$ is the probability level. Equivalently, a subset
$\mathcal{S}$ is acceptable if and only if there exists a lower
$(1-\varepsilon)$ posterior prediction interval for
$\widetilde{\mathcal{D}}_{\mathcal{S}_{min},\mathcal{S}}$ that includes $\eta$
(Kowal, 2021b, ). Subsets are _not_ acceptable if there is insufficient
predictive probability under $\mathcal{M}$ that the out-of-sample accuracy is
within a predetermined margin of the “best” subset. $\mathcal{S}_{min}$ is
necessarily a member of $\mathbb{A}_{\eta,\epsilon}$ for any
$(\eta,\varepsilon)$, so the acceptable family is always nonempty. Larger
values of $\eta$ and smaller values of $\varepsilon$ expand the acceptable
family; we select $\eta=0$ and $\varepsilon=0.10$ by default and conduct
sensitivity analyses (see also Kowal, 2021b, ; Kowal et al.,, 2021; Kowal,
2021a,  for further sensitivity evaluations).

The acceptable family is related to fence methods for model selection (Jiang
et al.,, 2008), which seek to eliminate “incorrect” models using likelihood
criteria. These methods rely on asymptotic arguments or bootstrap
computations, while our approach emphasizes out-of-sample predictive
performance with (finite-sample) predictive uncertainty quantification under
the LMM. Fence methods primarily focus on selection of a _single_ model,
rather than analysis of the _collection_ of near-optimal models or subsets.
Most critically, in our empirical examples the fence methods failed to
converge for $p\geq 10$ (using the R package fence), while the proposed
approach is highly scalable in both $n$ and $p$.

We summarize the acceptable family using two strategies. First, we report two
key members: the “best” subset $\mathcal{S}_{min}$ and the _smallest_
acceptable subset,

$\mathcal{S}_{small}\coloneqq\arg\min_{\mathcal{S}\in\mathbb{A}_{\eta,\varepsilon}}|\mathcal{S}|,$
(18)

which is the smallest subset of covariates that satisfies the near-optimality
condition in (17). Typically, we find
$|\mathcal{S}_{small}|\ll|\mathcal{S}_{min}|$ which is expected: selection
based on minimizing cross-validation error is known to produce models that are
unnecessarily complex. Although we caution against overreliance on any single
subset, $\mathcal{S}_{small}$ is a uniquely valuable summary of the acceptable
family: smaller subsets are _not_ acceptable, and therefore
$\mathcal{S}_{small}$ offers a notion of the “necessary” variables for near-
optimal prediction. When $\mathcal{S}_{small}$ is nonunique, the acceptable
subsets of size $|\mathcal{S}_{small}|$ should be evaluated in concert; a
unique choice of $\mathcal{S}_{small}$ is given by the acceptable subset of
this size that achieves the smallest empirical loss (14).

Second, we summarize $\mathbb{A}_{\eta,\varepsilon}$ using the variable
importance metric for each covariate $j$:

$\mbox{VI}_{\rm
incl}(j)\coloneqq|\mathbb{A}_{\eta,\varepsilon}|^{-1}\sum_{\mathcal{S}\in\mathbb{A}_{\eta,\varepsilon}}\mathbb{I}\\{j\in\mathcal{S}\\},$
(19)

which can also be generalized for two or more covariates (Kowal, 2021a, ).
This quantity is most informative at each endpoint: $\mbox{VI}_{\rm
incl}(j)\approx 1$ implies that covariate $j$ belongs to (nearly) all
acceptable subsets and is therefore an essential or _keystone_ covariate,
while $\mbox{VI}_{\rm incl}(j)\approx 0$ suggests that covariate $j$ is
irrelevant for (nearly) all acceptable subsets. By design, $\mbox{VI}_{\rm
incl}(j)$ provides a _variable-specific_ summary of the acceptable family of
subsets. This metric is broadly related to stability selection (Meinshausen
and Bühlmann,, 2010) and frequentist aggregation of variable importance across
many “good” models (Dong and Rudin,, 2019).

To compute the out-of-sample quantities in (14) and (16) under $\mathcal{M}$,
we use an importance sampling algorithm. This algorithm requires only the _in-
sample_ posterior under the LMM and hence avoids the intensive processing of
re-fitting $\mathcal{M}$ for each of the $K$ folds. The algorithm is detailed
in the supporting information and modifies previous approaches (Kowal, 2021b,
; Kowal et al.,, 2021; Kowal, 2021a, ) for LMMs and Mahalanobis loss.

###  2.5 Predictive uncertainty quantification for each action

For any subset $\mathcal{S}$, we provide uncertainty quantification for the
optimal linear coefficients $\bm{\delta}_{\mathcal{S}}$ using the predictive
distribution under $\mathcal{M}$. Specifically, we modify (5) to remove the
expectation under
$p_{\mathcal{M}}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\mid\bm{y}\\}$
and therefore preserve the predictive uncertainty quantification:

$\bm{\tilde{\delta}}_{\mathcal{S}}\coloneqq\arg\min_{\bm{\delta}_{\mathcal{S}}}\mathbb{E}_{[\bm{\psi}\mid\bm{y}]}\mathcal{L}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}),\bm{\delta}_{\mathcal{S}};\bm{\psi}\\}=(\bm{\tilde{X}}_{\mathcal{S}}^{\prime}\bm{\hat{\Omega}}\bm{\tilde{X}}_{\mathcal{S}})^{-1}\bm{\tilde{X}}_{\mathcal{S}}^{\prime}\bm{\hat{\Omega}}\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}}).$
(20)

This mechanism for uncertainty quantification generalizes the predictive
projection approach from Kowal, 2021a  to account for Mahalanobis loss. In
particular, (20) includes marginalization over $\bm{\Omega_{\psi}}$ to ensure
that the resulting quantity is exclusively a _posterior predictive_ variable
with a distribution induced by
$p_{\mathcal{M}}\\{\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})\mid\bm{y}\\}$.
However, (20) can be modified to include the uncertainty of
$\bm{\Omega_{\psi}}$ by replacing $\bm{\hat{\Omega}}$ with
$\bm{\Omega_{\psi}}$. Posterior samples of $\bm{\tilde{\delta}}_{\mathcal{S}}$
only require posterior predictive samples of
$\bm{\tilde{y}}(\bm{\tilde{X}},\bm{\tilde{Z}})$—which can be shared among all
subsets $\mathcal{S}$ of interest—and the solution to a GLS problem (20). In
particular, we use (20) to compute interval estimates for the linear
coefficients associated with $\mathcal{S}_{min}$ and $\mathcal{S}_{small}$.

##  3 Simulation study

We evaluate the proposed LMM subset selection techniques using simulated data
from a Gaussian random intercept model. First, we generate $p$ correlated
fixed effects covariates from marginal standard normal distributions with
$\mbox{Cor}(x_{i,j},x_{i,j^{\prime}})=(0.75)^{|j-j^{\prime}|}$ for
$i=1,\ldots,n$ and $j=1,\ldots,p$. The $p$ columns are randomly permuted and
augmented with an intercept. The true linear coefficients $\bm{\beta}^{*}$ are
constructed by setting $\beta_{0}^{*}=-1$ and fixing $p_{*}=5$ nonzero
coefficients, with $\lceil p_{*}/2\rceil$ equal to $1$ and $\lfloor
p_{*}/2\rfloor$ equal to $-1$, and the rest at zero. Let
$y_{i}^{*}\coloneqq\bm{x}_{i}^{\prime}\bm{\beta}^{*}$ denote the true
expectation. For a given intraclass correlation $\rho_{*}$ and signal-to-noise
ratio SNR, define
$\sigma_{tot}^{2}\coloneqq\mbox{var}(\\{y_{i}^{*}\\}_{i=1}^{n})/\mbox{SNR}$
and let $\sigma_{u}^{2}\coloneqq{\rho_{*}\sigma_{tot}^{2}}$ and
$\sigma_{\epsilon}^{2}\coloneqq{\sigma_{tot}^{2}-\sigma_{u}^{2}}$. The data
are generated as $y_{ij}=y_{i}^{*}+u_{i}+\epsilon_{ij}$ where
$u_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{u}^{2})$ and
$\epsilon_{ij}\sim N(0,\sigma_{\epsilon}^{2})$ for $j=1,\ldots,m$ and
$i=1,\ldots,n$. We consider $\rho_{*}=0.25$, $m=4$, $\mbox{SNR}=1$ (see the
supporting information for SNR = 5), $p\in\\{15,200\\}$, and
$n\in\\{75,150,300\\}$. We repeat the data-generating process 100 times for
each design.

We implement a Bayesian LMM using the sampler from Section 2.2.1 and horseshoe
priors on the fixed effects (see Section 4). Using $\mathcal{M}$, we extract
the acceptable family $\mathbb{A}_{0,0.10}$ with $s_{k}=15$. We compute point
predictions for all acceptable subsets and evaluate $\mathcal{S}_{min}$ and
$\mathcal{S}_{small}$ for variable selection and uncertainty quantification
via the 90% intervals from $\bm{\tilde{\delta}}_{\mathcal{S}}$. The primary
Bayesian competitor is given by the usual actions under $\mathcal{M}$:
posterior expectations for point predictions, 90% highest posterior density
(HPD) intervals of $\bm{\beta}$ for uncertainty quantification, and selection
based on whether the 95% HPD intervals for each $\beta_{j}$ exclude zero. As a
secondary Bayesian competitor, we compute the point predictions and interval
estimates from $\mathcal{S}_{small}$ under a (non-LMM) Gaussian linear
regression model using squared error loss (Kowal, 2021a, ), which ignores the
longitudinal aspect of the data. Lastly, we compare against classical
selection methods that do not account for the random effects. Specifically, we
apply the adaptive lasso (tuning parameter selected via 10-fold cross-
validation and the one-standard-error rule) and classical subset selection
(using AIC) to the data $\\{(\bm{x}_{i},\bar{y}_{i})\\}_{i=1}^{n}$ for
$\bar{y}_{i}\coloneqq m^{-1}\sum_{j=1}^{m}y_{ij}$. We attempted to include
fence-based variable selection for LMMs (Jiang et al.,, 2008), but the R
package fence failed to converge for any simulation settings with $p\geq 10$.

Point prediction accuracy is evaluated using Mahalanobis loss for $y_{i}^{*}$,
where the weight matrix (9) uses the true parameters for
$\sigma_{\epsilon}^{2}$ and $\sigma_{u}^{2}$. The simulation-averaged results
are in Table 1. $\mathcal{S}_{small}$ consistently provides the best or near-
best point predictions, followed by the posterior mean under $\mathcal{M}$.
Further, $\mathcal{S}_{small}$ usually selects fewer variables than all
competitors (not shown). Hence, $\mathcal{S}_{small}$ offers substantial
reductions in the subset size while maintaining near-optimal prediction
accuracy—which is precisely the goal of the smallest acceptable subset.

$(n,p)$ | lasso | subset | $\mathcal{M}$ | SE | $\mathcal{S}_{min}$ | $\mathcal{S}_{small}$ | $\mathbb{A}^{0}$ | $\mathbb{A}^{0.1}$ | $\mathbb{A}^{0.5}$ | $\mathbb{A}^{0.9}$ | $\mathbb{A}^{1}$
---|---|---|---|---|---|---|---|---|---|---|---
$(75,15)$ | 0.182 | 0.116 | 0.108 | 0.100 | 0.109 | 0.097 | 0.070 | 0.099 | 0.109 | 0.117 | 0.170
$(75,200)$ | 0.507 | 0.644 | 0.426 | 0.351 | 0.410 | 0.315 | 0.279 | 0.375 | 0.399 | 0.410 | 0.444
$(150,15)$ | 0.104 | 0.056 | 0.050 | 0.042 | 0.053 | 0.037 | 0.032 | 0.048 | 0.054 | 0.058 | 0.066
$(150,200)$ | 0.275 | 0.341 | 0.127 | 0.151 | 0.194 | 0.111 | 0.108 | 0.160 | 0.184 | 0.196 | 0.201
$(300,15)$ | 0.057 | 0.025 | 0.023 | 0.019 | 0.025 | 0.015 | 0.014 | 0.022 | 0.026 | 0.027 | 0.028
$(300,200)$ | 0.096 | 0.203 | 0.048 | 0.078 | 0.105 | 0.056 | 0.054 | 0.083 | 0.099 | 0.107 | 0.110
Table 1:  Left: average Mahalanobis loss for each method; the smallest (best)
value is bolded. Right: for each simulation, we compute the $q$th quantile of
the true Mahalanobis loss for predicting $y_{i}^{*}$ for all acceptable
subsets, and then average that quantity across simulations to obtain
$\mathbb{A}^{q}$. $\mathcal{M}$ refers to the posterior mean under the
Bayesian LMM and SE is the smallest acceptable subset using linear regression
and squared error loss (Kowal, 2021a, ). $\mathcal{S}_{small}$ performs
exceptionally well across all settings, while even the worst acceptable
subsets $(\mathbb{A}^{1})$ outperform the frequentist competitors.

More broadly, we evaluate the overall predictive performance of the acceptable
family $\mathbb{A}_{0,0.10}$ using the $q$th quantile of the true Mahalanobis
loss for each acceptable subset at each simulation, and then average that
quantity across simulations to obtain $\mathbb{A}^{q}$. For example,
$\mathbb{A}^{1}$ is the worst possible performance in the acceptable family,
i.e., if an oracle were to select the _worst acceptable_ subset at each
simulation. The main takeaways (see Table 1) are (i) $\mathcal{S}_{small}$
typically outperforms even $\mathbb{A}^{0.1}$ and therefore is consistently in
the top 10% of acceptable subsets, and (ii) even the worst acceptable subsets
outperform the frequentist competitors. These results confirm our notion of
_near-optimality_ of the acceptable family.

The 90% interval estimates for $\bm{\beta}^{*}$ are evaluated in Figure 1,
which reports the mean interval widths and the empirical coverage; narrow
intervals that provide the correct nominal coverage are preferred. The
intervals from $\mathcal{S}_{small}$ are clearly the best for these cases: the
intervals maintain 90% coverage and are much narrower than competing methods.
In particular, the intervals from both $\mathcal{S}_{small}$ and
$\mathcal{S}_{min}$ are far more precise (i.e., narrower) than the 90% HPD
intervals under $\mathcal{M}$.

Figure 1:  Mean 90% interval widths (boxplots) with empirical coverage
(annotations) for $\bm{\beta}^{*}$. Non-overlapping notches indicate
significant differences between medians. The proposed intervals based on
$\mathcal{S}_{small}$ are significantly narrower than the usual HPD intervals
under $\mathcal{M}$ yet maintain the nominal 90% coverage.

Lastly, we evaluate the (marginal) selection capabilities using true positive
rates (TPRs) and true negative rates (TNRs) in Table 2. $\mathcal{S}_{small}$
provides consistently high TPRs and TNRs, while the 95% HPD intervals under
$\mathcal{M}$ are far too conservative for selection (low TPRs). Both
selection mechanisms are based on the same Bayesian LMM $\mathcal{M}$, but
$\mathcal{S}_{small}$ is decisively better. The improvements over classical
subset selection are also substantial.

$(n,p)$ |  | lasso | subset | posterior HPD | $\mathcal{S}_{min}$ | $\mathcal{S}_{small}$
---|---|---|---|---|---|---
$(75,15)$ | TPR | 0.95 | 0.98 | 0.86 | 0.99 | 0.95
TNR | 0.93 | 0.80 | 0.98 | 0.63 | 0.95
$(75,200)$ | TPR | 0.94 | 0.91 | 0.57 | 0.94 | 0.91
TNR | 0.96 | 0.92 | 1.00 | 0.91 | 0.98
$(150,15)$ | TPR | 0.99 | 1.00 | 0.99 | 1.00 | 0.99
TNR | 0.98 | 0.82 | 0.99 | 0.68 | 0.97
$(150,200)$ | TPR | 0.99 | 1.00 | 0.95 | 1.00 | 1.00
TNR | 0.97 | 0.93 | 1.00 | 0.92 | 0.98
$(300,15)$ | TPR | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
TNR | 1.00 | 0.85 | 0.99 | 0.69 | 0.98
$(300,200)$ | TPR | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
TNR | 0.99 | 0.92 | 1.00 | 0.92 | 0.98
Table 2:  True positive rates (TPR) and true negative rates (TNR) for
synthetic data with $p_{*}+1=6$ active covariates including the intercept.
$\mathcal{S}_{small}$ provides consistently high TPRs and TNRs, while the 95%
HPD intervals under $\mathcal{M}$ are far too conservative (low TPRs for
smaller $n$).

##  4 Application

We apply our subset selection analysis to moderate-to-vigorous physical
activity (MVPA) data from NHANES 2005-2006. Intraday activity was measured on
each subject using hip-worn accelerometers for one to seven days.
$\mbox{MVPA}_{ij}$ is defined as the number of minutes with at least 2020
activity counts for subject $i$ on day $j$, and typically corresponds to more
intensive activities that include vigorous walking or running (Fishman et
al.,, 2016). The goal is to determine the subject-specific factors that
predict MVPA. However, these longitudinal data feature repeated measurements
on each participant, and this within-subject dependence must be accounted for
in both modeling and decision analysis.

We specifically analyze older (ages 65-80) and Hispanic (Mexican American or
Other Hispanic) individuals. Fixed effects include body mass index (BMI), age,
gender (male or female), education level (less than high school, completed
high school only, or some college and above), total cholesterol, HDL
cholesterol, systolic blood pressure, smoking status (never, former, or
current), drinking status (never, moderate, or heavy), and presence of
diabetes. After filtering to individuals with at least one day of activity
data, days with at least 10 hours of accelerometer wear time, activity
measurements that were correctly “calibrated” and “reliable” as flagged by
NHANES, and individuals with no mobility problems, the resulting analysis
dataset has $N=243$ measurements on $n=61$ individuals with $p=13$ covariates.

We model $y_{ij}=\log(\mbox{MVPA}_{ij}+1)$ using a Gaussian random intercept
model (see Section 2.2.1) with horseshoe priors for the fixed effects
$\bm{\beta}$, a Jeffreys prior for $\sigma_{\epsilon}^{2}$, and a uniform
prior for $\sigma_{u}\sim\mbox{Unif}(0,100)$. The Gibbs sampler generated
10,000 samples after a burn-in of 5,000; traceplots indicated no lack of
convergence and the effective sample sizes were sufficiently large. A 95% HPD
interval for the within-subject correlation,
$\sigma_{u}^{2}/(\sigma_{u}^{2}+\sigma_{\epsilon}^{2})$, is $(0.18,0.48)$,
which suggests moderate within-subject autocorrelation.

Using the posterior and predictive samples from the Bayesian LMM, we compute
and study the acceptable family. Since $p$ is not large, we filter from the
$2^{p}=4096$ possible subsets to the “best” $s_{k}=100$ models of each size
$k=1,\ldots,s_{max}=p+1$ (the intercept is always included), which produces
973 candidate subsets. Figure 2 summarizes the predictive performance among
these candidates using Mahalanobis predictive loss. The 80% intervals that
include $\eta=0$ (horizontal line) correspond to acceptable subsets with
$\varepsilon=0.1$. Each subset of size two performs 10-35% worse than
$\mathcal{S}_{min}$. One subset of size three outperforms the rest and is
within 3-7% of $\mathcal{S}_{min}$; this subset would be acceptable for
margins $\eta>3\%$ or smaller $\varepsilon\leq 0.01$, which corresponds to
wider intervals in Figure 2. The smallest acceptable subset for $\eta=0\%$ has
four variables, and notably performs as well or better than the larger
subsets; this subset is also unchanged for $\varepsilon\in[0.02,0.16]$.

Figure 2:  Prediction 80% intervals (lines) and expectations (points) for
$\widetilde{\mathcal{D}}_{\mathcal{S}_{min},\mathcal{S}}$ across candidate
subsets $\mathcal{S}$ (including the intercept). The analogous empirical
quantity, ${\mathcal{D}}_{\mathcal{S}_{min},\mathcal{S}}\coloneqq
100\times({\mathcal{L}}_{\mathcal{S}}-{\mathcal{L}}_{\mathcal{S}_{min}})/{\mathcal{L}}_{\mathcal{S}_{min}}$,
is also denoted (x-marks). Subsets of the same size are jittered for clarity
of presentation. Intervals that include $\eta=0$ (horizontal line) correspond
to acceptable subsets. Acceptable subsets range from sizes 4 to 12, including
$|\mathcal{S}_{small}|=4$ (solid line) and $|\mathcal{S}_{min}|=6$ (dashed
line).

The acceptable family features $|\mathbb{A}_{0,0.10}|=333$ members ranging
from sizes 4 to 12. To summarize $\mathbb{A}_{0,0.10}$, we report the variable
importance metric $\mbox{VI}_{\rm incl}(j)$ in Figure 3. Gender, age, and
total cholesterol are keystone covariates that appear in all acceptable
subsets, and are the only members (plus the intercept) of
$\mathcal{S}_{small}$. Notably, the remaining covariates appear in some—but
not most—of the acceptable subsets. These covariates are not entirely
extraneous, but appear to be interchangeable and not strictly necessary for
acceptable linear prediction. The “best” subset $\mathcal{S}_{min}$ adds
smoking status (current) and diabetes to $\mathcal{S}_{small}$. Yet the
variable importance provides important context for $\mathcal{S}_{min}$:
although smoking status belongs to the “best” subset, it only appears in a
moderate fraction (about 40%) of the acceptable subsets. By comparison,
education level (some college and above) appears in vastly more acceptable
subsets, yet does not belong to $\mathcal{S}_{small}$.

Figure 3:  Variable importance $\mbox{VI}_{\rm incl}(j)$ for the acceptable
family $\mathbb{A}_{0,0.10}$. Gender, age, and total cholesterol are keystone
covariates, while the other covariates appear in some—but not most—of the
acceptable subsets. These results are robust for $\varepsilon\in[0.02,0.16]$.

Lastly, Figure 4 compares the point and interval estimates from
$\mathcal{S}_{small}$ against the Bayesian LMM $\mathcal{M}$ and the adaptive
lasso. Both $\mathcal{S}_{small}$ and $\mathcal{M}$ highlight a positive
effect for total cholesterol—perhaps a realization of the common advice that
individuals with high cholesterol should exercise more—while all three methods
agree on negative effects for gender (female) and age. $\mathcal{S}_{small}$
produces narrower intervals among the nonzero coefficients compared to the HPD
intervals under $\mathcal{M}$, and offers a sparsity in point estimation that
is not available for the posterior means under $\mathcal{M}$. Yet the methods
broadly agree: the selected variables in $\mathcal{S}_{small}$ correspond
exactly to the 90% HPD intervals under $\mathcal{M}$ that exclude zero. The
frequentist intervals from Zhao et al., (2021) are difficult to interpret,
since they often fail to include the (adaptive) lasso-based point estimates
from which they were derived.

Figure 4:  Estimated coefficients and 90% intervals for $\mathcal{S}_{small}$
($\bm{\hat{\delta}}_{\mathcal{S}_{small}}$ and quantiles from
$\bm{\tilde{\delta}}_{\mathcal{S}_{small}}$), the Bayesian LMM $\mathcal{M}$
(posterior means and HPD intervals for $\bm{\beta}$), and the adaptive lasso
(intervals from Zhao et al.,, 2021). There is broad agreement, although
$\mathcal{S}_{small}$ produces narrower intervals and sparse estimates
compared to the the usual Bayesian LMM estimates. These results are robust for
$\varepsilon\in[0.02,0.16]$.

The results are robust to $\varepsilon$: $\mathcal{S}_{small}$ is unchanged
for $\varepsilon\in[0.02,0.16]$, and the variable importances are stable.
$\mathcal{S}_{small}$ omits total cholesterol for $\varepsilon=0.01$ and adds
diabetes for $\varepsilon=0.20$. The number of acceptable subsets decreases
from $|\mathbb{A}_{0,\epsilon}|\in\\{658,456,333,60\\}$ for
$\varepsilon\in\\{0.01,0.05,0.10,0.20\\}$, which is expected: larger values of
$\varepsilon$ provide more lenient admission to the acceptable family.

Note that NHANES data are collected from a complex sampling design, and
population-level inference typically requires survey adjustments. The
oversampled groups in NHANES 2005-2006 are specific age groups (12-19 and 60+
years), races (Black and Mexican Americans), and low-income individuals.
Because we subset by age group and race and further include age and many other
covariates in the model, we expect that the effects of the sampling design are
mitigated.

##  5 Discussion

We have developed a decision analysis strategy for subset selection in
Bayesian LMMs. Using a Mahalanobis predictive loss function to bring forward
the structured dependence from the LMM into the decision analysis, we derived
optimal linear coefficients for (i) any _given_ subset of variables and (ii)
all subsets of variables that satisfy a _cardinality constraint_. The
coefficients are accompanied by predictive uncertainty quantification and
regularization inherited from the underlying Bayesian LMM. Comparing across
subsets, we collected and summarized the _acceptable family_ of subsets that
(nearly) matched the predictive performance of the “best” subset. The proposed
tools demonstrated excellent prediction, estimation, and selection properties
on simulated data, and were applied to a longitudinal dataset to study the key
predictors of MVPA.

Given the acceptable family of near-optimal subsets, it is natural to ask:
“Which subset should be used?” Our response is that, based on predictive
accuracy, _any_ of the acceptable subsets provides a reasonable answer. Absent
additional information (such as individual variable costs), we advocate the
_smallest_ acceptable subset $\mathcal{S}_{small}$, which simultaneously (i)
provides excellent prediction, uncertainty quantification, and selection
capabilities across a variety of challenging simulation settings and (ii)
offers a notion of the necessary variables for near-optimal linear prediction
(when $\mathcal{S}_{small}$ is unique), since smaller subsets are _not_
acceptable by definition. However, our prioritization of P4 underlines the
crucial point that no _single_ subset—including $\mathcal{S}_{small}$—should
be used in isolation to report the variables that “matter”. In particular,
variables excluded from the “best” subset are not necessarily irrelevant,
while variables included in the “best” subset are not necessarily essential.
The acceptable family fills in those gaps to provide a more complete picture,
and is accompanied by suitable summaries.

The Mahalanobis loss (3) is designed for the LMM (1), which is most commonly a
Gaussian LMM. Although we focused primarily on random intercept and random
slope models, the results are broadly applicable among LMMs, including many
functional data and spatial regression models. In addition, modifications for
non-Gaussian _generalized_ LMMs (GLMMs) may be attainable. For Bayesian subset
selection with binary data, Kowal, 2021a  used iteratively-reweighted least
squares (IRLS) to approximate the minimizer of a cross-entropy loss with a
weighted least squares solution. IRLS is a widely popular strategy for
estimating generalized linear models, and can be used to produce optimal
coefficients under the corresponding deviance loss functions. For LMMs, a
natural modification would be to insert a weighting matrix akin to
$\bm{\Omega_{\psi}}$ into the IRLS, thereby extending the proposed tools for
compatibility with GLMMs.

### Acknowledgements

Research was sponsored by the Army Research Office and was accomplished under
Grant Number W911NF-20-1-0184. The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein.

## References

  * Bashir et al., (2019) Bashir, A., Carvalho, C. M., Hahn, P. R., and Jones, M. B. (2019).  Post-processing posteriors over precision matrices to produce sparse graph estimates.  Bayesian Analysis, 14(4):1075–1090.
  * Bertsimas et al., (2016) Bertsimas, D., King, A., and Mazumder, R. (2016).  Best subset selection via a modern optimization lens.  Annals of statistics, 44(2):813–852.
  * Bhattacharya et al., (2016) Bhattacharya, A., Chakraborty, A., and Mallick, B. K. (2016).  Fast sampling with Gaussian scale mixture priors in high-dimensional regression.  Biometrika, 103(4):985–991.
  * Bondell et al., (2010) Bondell, H. D., Krishna, A., and Ghosh, S. K. (2010).  Joint variable selection for fixed and random effects in linear mixed‐effects models.  Biometrics, 66(4):1069–1077.
  * Chen and Dunson, (2003) Chen, Z. and Dunson, D. B. (2003).  Random effects selection in linear mixed models.  Biometrics, 59(4):762–769.
  * Dong and Rudin, (2019) Dong, J. and Rudin, C. (2019).  Variable importance clouds: A way to explore variable importance for the set of good models.  arXiv preprint arXiv:1901.03209.
  * Fan and Li, (2012) Fan, Y. and Li, R. (2012).  Variable selection in linear mixed effects models.  Annals of statistics, 40(4):2043.
  * Fishman et al., (2016) Fishman, E. I., Steeves, J. A., Zipunnikov, V., Koster, A., Berrigan, D., Harris, T. A., and Murphy, R. (2016).  Association between Objectively Measured Physical Activity and Mortality in NHANES.  Medicine and Science in Sports and Exercise, 48(7):1303–1311.
  * Foster et al., (2007) Foster, S. D., Verbyla, A. P., and Pitchford, W. S. (2007).  Incorporating LASSO effects into a mixed model for quantitative trait loci detection.  Journal of agricultural, biological, and environmental statistics, 12(2):300–314.
  * Furnival and Wilson, (2000) Furnival, G. M. and Wilson, R. W. (2000).  Regressions by leaps and bounds.  Technometrics, 42(1):69–79.
  * Hahn and Carvalho, (2015) Hahn, P. R. and Carvalho, C. M. (2015).  Decoupling shrinkage and selection in bayesian linear models: A posterior summary perspective.  Journal of the American Statistical Association, 110(509):435–448.
  * Hastie et al., (2020) Hastie, T., Tibshirani, R., and Tibshirani, R. (2020).  Best Subset, Forward Stepwise or Lasso? Analysis and Recommendations Based on Extensive Comparisons.  Statistical Science, 35(4):579–592.
  * Ibrahim et al., (2011) Ibrahim, J. G., Zhu, H., Garcia, R. I., and Guo, R. (2011).  Fixed and random effects selection in mixed effects models.  Biometrics, 67(2):495–503.
  * Jiang et al., (2008) Jiang, J., Rao, J. S., Gu, Z., and Nguyen, T. (2008).  Fence methods for mixed model selection.  The Annals of Statistics, 36(4):1669–1692.
  * Kinney and Dunson, (2007) Kinney, S. K. and Dunson, D. B. (2007).  Fixed and random effects selection in linear and logistic models.  Biometrics, 63(3):690–698.
  * (16) Kowal, D. R. (2021a).  Bayesian subset selection and variable importance for interpretable prediction and classification.  arXiv preprint arXiv:2104.10150.
  * (17) Kowal, D. R. (2021b).  Fast, Optimal, and Targeted Predictions using Parametrized Decision Analysis.  Journal of the American Statistical Association.
  * Kowal and Bourgeois, (2020) Kowal, D. R. and Bourgeois, D. C. (2020).  Bayesian Function-on-Scalars Regression for High-Dimensional Data.  Journal of Computational and Graphical Statistics, 29(3):1–10.
  * Kowal et al., (2021) Kowal, D. R., Bravo, M., Leong, H., Griffin, R. J., Ensor, K. B., and Miranda, M. L. (2021).  Bayesian Variable Selection for Understanding Mixtures in Environmental Exposures.  Statistics in Medicine, 40(22):4850–4871.
  * Lindley, (1968) Lindley, D. V. (1968).  The Choice of Variables in Multiple Regression.  Journal of the Royal Statistical Society: Series B (Methodological), 30(1):31–53.
  * Meinshausen and Bühlmann, (2010) Meinshausen, N. and Bühlmann, P. (2010).  Stability selection.  Journal of the Royal Statistical Society. Series B: Statistical Methodology, 72(4):417–473.
  * Miller, (1984) Miller, A. J. (1984).  Selection of subsets of regression variables.  Journal of the Royal Statistical Society: Series A (General), 147(3):389–410.
  * Müller et al., (2013) Müller, S., Scealy, J. L., and Welsh, A. H. (2013).  Model selection in linear mixed models.  Statistical Science, 28(2):135–167.
  * Nishimura and Suchard, (2018) Nishimura, A. and Suchard, M. A. (2018).  Prior-preconditioned conjugate gradient method for accelerated Gibbs sampling in ”large n and large p” Bayesian sparse regression.  arXiv preprint arXiv:1810.12437.
  * Puelz et al., (2017) Puelz, D., Hahn, P. R., and Carvalho, C. M. (2017).  Variable selection in seemingly unrelated regressions with random predictors.  Bayesian Analysis, 12(4):969–989.
  * Wang et al., (2011) Wang, D., Eskridge, K. M., and Crossa, J. (2011).  Identifying QTLs and epistasis in structured plant populations using adaptive mixed LASSO.  Journal of agricultural, biological, and environmental statistics, 16(2):170–184.
  * Zhao et al., (2021) Zhao, S., Witten, D., and Shojaie, A. (2021).  In defense of the indefensible: A very naive approach to high-dimensional inference.  Statistical Science, 36(4):562–577.