Efficient Computation of High-Dimensional Penalized Generalized Linear Mixed
Models by Latent Factor Modeling of the Random Effects–A.2

# Efficient Computation of High-Dimensional Penalized Generalized Linear Mixed
Models by Latent Factor Modeling of the Random Effects

Hillary M. Heiling1,∗<EMAIL_ADDRESS>Naim U. Rashid1 Quefeng Li1 Xianlu
L. Peng2 Jen Jen Yeh2,3,4
and Joseph G. Ibrahim1
1Department of Biostatistics  University of North Carolina Chapel Hill  Chapel
Hill  NC
2Lineberger Comprehensive Cancer Center  University of North Carolina at
Chapel Hill  Chapel Hill  NC
3Department of Surgery  University of North Carolina Chapel Hill  Chapel Hill
NC
4Department of Pharmacology  University of North Carolina Chapel Hill  Chapel
Hill  NC

(2022)

###### Abstract

Modern biomedical datasets are increasingly high dimensional and exhibit
complex correlation structures. Generalized Linear Mixed Models (GLMMs) have
long been employed to account for such dependencies. However, proper
specification of the fixed and random effects in GLMMs is increasingly
difficult in high dimensions, and computational complexity grows with
increasing dimension of the random effects. We present a novel reformulation
of the GLMM using a factor model decomposition of the random effects, enabling
scalable computation of GLMMs in high dimensions by reducing the latent space
from a large number of random effects to a smaller set of latent factors. We
also extend our prior work to estimate model parameters using a modified Monte
Carlo Expectation Conditional Minimization algorithm, allowing us to perform
variable selection on both the fixed and random effects simultaneously. We
show through simulation that through this factor model decomposition, our
method can fit high dimensional penalized GLMMs faster than comparable methods
and more easily scale to larger dimensions not previously seen in existing
approaches.

###### keywords:

Factor model; Generalized linear mixed models; variable selection

##  1 Introduction

Modern biomedical datasets are increasingly high dimensional and exhibit
complex correlation structures. Generalized Linear Mixed Models (GLMMs) have
long been employed to account for such dependencies. However, proper
specification of the fixed and random effects is a critical step in estimation
of GLMMs. For instance, omitting important random effects can lead to bias in
the estimated variance of the fixed effects, whereas including unnecessary
random effects could increase the computational difficulty of fitting the GLMM
(Thompson et al., 2017; Gurka et al., 2011; Bondell et al., 2010).

Despite the importance of properly specifying the set of fixed and random
effects in such models, it is often unknown a priori which variables should be
specified as fixed or random in the model, particularly in high dimensional
settings in which the feature space is generally assumed to be sparse. As a
solution to this problem, Rashid et al. (2020) used a Monte Carlo Expectation
Conditional Minimization (MCECM) algorithm to perform variable selection in
high dimensional GLMMs. Their method has since been incorporated into the
glmmPen R package (Heiling et al., 2023). In this glmmPen framework,
performing simultaneous variable selection on fixed and random effects in
GLMMs with input dimensions of 50 or so predictors is feasible. Although this
glmmPen framework extends the feasible dimensionality of performing variable
selection within GLMMs relative to existing methods, new methodology is needed
to alleviate the computational burden as the dimension increases even further
and allow scalability to hundreds or thousands of predictors.

We present a novel reformulation of the GLMM, which we call glmmPen_FA, using
a factor model decomposition of the random effects. This factor model is used
as a dimension reduction tool to represent a large number of latent random
effects as a function of a smaller set of latent factors. By reducing the
latent space of the random effects, this new model formulation enables us to
extend the feasible dimensionality of performing variable selection in GLMMs
to hundreds of predictors. We estimate model parameters and perform
simultaneous selection of fixed and random effects using a Monte Carlo
Expectation Conditional Minimization (MCECM) algorithm. We show through
simulations that through this factor model decomposition, our method can fit
high dimensional penalized GLMMs (pGLMMs) faster than comparable methods and
more easily scale to larger dimensions not previously seen in existing
approaches.

The remainder of this paper is organized as follows. Section 2 reviews the
statistical models and algorithm used to estimate pGLMMs in our new factor
model decomposition framework, termed glmmPen_FA. In section 3, simulations
are conducted to assess the performance of the new glmmPen_FA method. Section
4 describes a motivating case study for the prediction of pancreatic cancer
subtypes using gene expression data and provides results from the application
of our new method to the case study. We close the article with some discussion
in Section 5.

##  2 Methods

###  2.1 Model formulation

In this section, we review the notation and model formulation of our approach.
We consider the case where we want to analyze data from $K$ independent groups
of observations. For each group $k=1,...,K$, there are $n_{k}$ observations
for a total sample size of $N=\sum_{k=1}^{K}n_{k}$. For group $k$, let
$\boldsymbol{y}_{k}=(y_{k1},...,y_{kn_{k}})^{T}$ be the vector of $n_{k}$
independent responses, $\boldsymbol{x}_{ki}=(x_{ki,1},...,x_{ki,p})^{T}$ be
the $p$-dimensional vector of predictors, and
$\boldsymbol{X}_{k}=(\boldsymbol{x}_{k1},...,\boldsymbol{x}_{kn_{k}})^{T}$.
For simplification of notation, we will set $n_{1}=...=n_{K}=n$ without loss
of generality. In GLMMs, we assume that the conditional distribution of
$\boldsymbol{y}_{k}$ given $\boldsymbol{X}_{k}$ belongs to the exponential
family and has the following density:

$f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k};\theta)=\prod_{i=1}^{n}c(y_{ki})\exp[\tau^{-1}\\{y_{ki}\eta_{ki}-b(\eta_{ki})\\}],$
(1)

where $c(y_{ki})$ is a constant that only depends on $y_{ki}$, $\tau$ is the
dispersion parameter, $b(\cdot)$ is a known link function, and $\eta_{ki}$ is
the linear predictor.

As outlined in Rashid et al. (2020), the traditional GLMM formulation of the
linear predictor can be represented as

$\eta_{ki}=\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{\gamma}_{k}=\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{\Gamma}\boldsymbol{\delta}_{k},$
(2)

where $\boldsymbol{\beta}=(\beta_{1},...,\beta_{p})^{T}$ is a $p$-dimensional
vector for the fixed effects coefficients (including the intercept),
$\boldsymbol{\Gamma}$ is the Cholesky decomposition of the random effects
covariance matrix $\boldsymbol{\Sigma}$ such that
$\boldsymbol{\Gamma}\boldsymbol{\Gamma}^{T}=\boldsymbol{\Sigma}$,
$\boldsymbol{\gamma}_{k}=\boldsymbol{\Gamma}\boldsymbol{\delta}_{k}$, where
$\boldsymbol{\delta}_{k}\sim N_{q}(0,I)$, is a $q$-dimensional vector of
unobservable random effects (including the random intercept) for group $k$,
and $\boldsymbol{z}_{ki}$ is a $q$-dimensional subvector of
$\boldsymbol{x}_{ki}$.

In its current representation, the model assumes a $q$-dimensional latent
space. When $q$ is large, the estimation of the covariance matrix
$\boldsymbol{\Sigma}$ = Var($\boldsymbol{\gamma}_{k}$), can be computationally
burdensome to compute due to both the number of parameters needed to estimate
this matrix (for an unstructured covariance matrix, $q(q+1)/2$ parameters are
needed) as well as the need to approximate a $q$-dimensional integral (see
Section 2.3 for details). Prior work such as Fan et al. (2013) and Tran et al.
(2020) have assumed a factor model structure in order to estimate high-
dimensional covariance matrices in other settings, such as the estimation of
sample covariance matrices for time series data and the covariance matrix in
variational inference used to approximate the posterior distribution,
respectively. Here we introduce a novel formulation of the GLMM where we
decompose the random effects $\boldsymbol{\gamma}_{k}$ into a factor model
with $r$ latent common factors ($r\ll q$) such that
$\boldsymbol{\gamma}_{k}=\boldsymbol{B}\boldsymbol{\alpha}_{k}$, where
$\boldsymbol{B}$ is the $q\times r$ loading matrix and
$\boldsymbol{\alpha}_{k}$ represents the $r$ latent common factors. We assume
the latent factors $\boldsymbol{\alpha}_{k}$ are uncorrelated and follow a
$N_{r}(\boldsymbol{0},\boldsymbol{I})$ distribution. We re-write the linear
predictor as

$\eta_{ki}=\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{B}\boldsymbol{\alpha}_{k}.$
(3)

In this representation, the random component of the linear predictor has
variance Var($\boldsymbol{B}\boldsymbol{\alpha}_{k}$) = $\boldsymbol{BB}^{T}$
= $\boldsymbol{\Sigma}$. By assuming that $\boldsymbol{\Sigma}$ is low rank,
we also reduce the dimension of the latent space from $q$ to $r$, which
reduces the dimension of the integral in the likelihood and thereby reduces
the computational complexity of the E-step in the EM algorithm. Further
details are given in Section 2.3.

In order to estimate $\boldsymbol{B}$, let
$\boldsymbol{b}_{t}\in\mathbb{R}^{r}$ be the $t$-th row of $\boldsymbol{B}$
and $\boldsymbol{b}=(\boldsymbol{b}_{1}^{T},...,\boldsymbol{b}_{q}^{T})^{T}$.
We can then reparameterize the linear predictor as

$\eta_{ki}=\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{B}\boldsymbol{\alpha}_{k}=\left(\boldsymbol{x}_{ki}^{T}\hskip
10.0pt(\boldsymbol{\alpha}_{k}\otimes\boldsymbol{z}_{ki})^{T}\boldsymbol{J}\right)\left(\begin{matrix}\boldsymbol{\beta}\\\
\boldsymbol{b}\end{matrix}\right)$ (4)

in a manner similar to Chen and Dunson (2003) and Ibrahim et al. (2011), where
$\boldsymbol{J}$ is a matrix that transforms $\boldsymbol{b}$ to
vec($\boldsymbol{B}$) such that
$\text{vec}(\boldsymbol{B})=\boldsymbol{J}\boldsymbol{b}$ and $\boldsymbol{J}$
is of dimension $(qr)\times(qr)$. The vector of parameters
$\boldsymbol{\theta}=(\boldsymbol{\beta}^{T},\boldsymbol{b}^{T},\tau)^{T}$ are
the main parameters of interest. We denote the true value of
$\boldsymbol{\theta}$ as
$\boldsymbol{\theta}^{*}=(\boldsymbol{\beta}^{*T},\boldsymbol{b}^{*T},\tau^{*})^{T}=\text{argmin}_{\boldsymbol{\theta}}\text{E}_{\boldsymbol{\theta}}[-\ell(\boldsymbol{\theta})]$
where $\ell(\boldsymbol{\theta})$ is the observed log-likelihood across all
$K$ groups such that
$\ell(\boldsymbol{\theta})=\sum_{k=1}^{K}\ell_{k}(\boldsymbol{\theta})$, where
$\ell_{k}(\boldsymbol{\theta})=(1/n)\log\int
f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k};\boldsymbol{\theta})\phi(\boldsymbol{\alpha}_{k})d\boldsymbol{\alpha}_{k}$.

Our main interest lies in selecting the true nonzero fixed and random effects.
In other words, we aim to identify the set

$S=S_{1}\cup S_{2}=\\{j:\beta_{j}^{*}\neq
0\\}\cup\\{t:||\boldsymbol{b}_{t}^{*}||_{2}\neq 0\\},$

where $S_{1}$ and $S_{2}$ represent the true fixed and random effects,
respectively. When $\boldsymbol{b}_{t}=\boldsymbol{0}$, this indicates that
effect of covariate $t$ is fixed across the $K$ groups (i.e. the corresponding
$t$-th row and column of $\boldsymbol{\Sigma}$ is set to $\boldsymbol{0}$).

We aim to solve the following penalized likelihood problem:

$\widehat{\boldsymbol{\theta}}=\text{argmin}_{\boldsymbol{\theta}}-\ell(\boldsymbol{\theta})+\lambda_{0}\sum_{j=1}^{p}\rho_{0}\left(\beta_{j}\right)+\lambda_{1}\sum_{t=1}^{q}\rho_{1}\left(||\boldsymbol{b}_{t}||_{2}\right),$
(5)

where $\ell(\boldsymbol{\theta})$ is the observed log-likelihood for all $K$
groups defined earlier, $\rho_{0}(t)$ and $\rho_{1}(t)$ are general folded-
concave penalty functions, and $\lambda_{0}$ and $\lambda_{1}$ are positive
tuning parameters. The $\rho_{0}(t)$ penalty functions could include the
$L_{1}$ penalty, the SCAD penalty, and the MCP penalty (Friedman et al., 2010;
Breheny and Huang, 2011). For the $\rho_{1}(t)$ penalty, we treat the elements
of $\boldsymbol{b}_{t}$ as a group and penalize them in a groupwise manner
using the group LASSO, group MCP, or group SCAD penalties presented by Breheny
and Huang (2015). These groups of $\boldsymbol{b}_{t}$ are then estimated to
be either all zero or all nonzero. In this way, we select covariates to have
varying effects ($\boldsymbol{\widehat{b}}_{t}\neq\boldsymbol{0}$) or fixed
effects ($\boldsymbol{\widehat{b}}_{t}=\boldsymbol{0}$) across the $K$ groups.

We standardize the fixed effects covariates matrix
$\boldsymbol{X}=(\boldsymbol{X}_{1}^{T},...,\boldsymbol{X}_{K}^{T})^{T}$ such
that
$\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}x_{ki,j}=0$ and
$N^{-1}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}x_{ki,j}^{2}=1$ for $j=1,...,p$.

###  2.2 MCECM algorithm

We solve (5) for some specific $(\lambda_{0},\lambda_{1})$ using a Monte Carlo
Expectation Conditional Minimization (MCECM) algorithm (Garcia et al., 2010).

In the $s^{th}$ iteration of the MCECM algorithm, we aim to evaluate the
expectation of (E-step) and minimize (M-step) the following penalized
Q-function:

$\displaystyle\begin{aligned}
Q_{\lambda}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})&=\sum_{k=1}^{K}E\left\\{-\log(f(\boldsymbol{y}_{k},\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k};\boldsymbol{\theta}|\boldsymbol{d}_{o};\boldsymbol{\theta}^{(s)}))\right\\}+\\\
&\hskip
10.0pt\lambda_{0}\sum_{j=1}^{p}\rho_{0}\left(\beta_{j}\right)+\lambda_{1}\sum_{t=1}^{q}\rho_{1}\left(||\boldsymbol{b}_{t}||_{2}\right)\\\
&=Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})+Q_{2}(\boldsymbol{\theta}^{(s)})+\\\
&\hskip
10.0pt\lambda_{0}\sum_{j=1}^{p}\rho_{0}\left(\beta_{j}\right)+\lambda_{1}\sum_{t=1}^{q}\rho_{1}\left(||\boldsymbol{b}_{t}||_{2}\right),\end{aligned}$
(6)

where $(\boldsymbol{y}_{k},\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k})$ gives
the complete data for group $k$,
$d_{k,o}=(\boldsymbol{y}_{k},\boldsymbol{X}_{k})$ gives the observed data for
group $k$, $\boldsymbol{d}_{o}$ represents the entirety of the observed data,
and

$Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})=-\sum_{k=1}^{K}\int\log[f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k};\boldsymbol{\theta})]\phi(\boldsymbol{\alpha}_{k}|\boldsymbol{d}_{k,o};\boldsymbol{\theta}^{(s)})d\boldsymbol{\alpha}_{k},$
(7)
$Q_{2}(\boldsymbol{\theta}^{(s)})=-\sum_{k=1}^{K}\int\log[\phi(\boldsymbol{\alpha}_{k}|\boldsymbol{d}_{k,o};\boldsymbol{\theta}^{(s)})]d\boldsymbol{\alpha}_{k}$
(8)

In the E-step of the algorithm, we aim to approximate the $r$-dimensional
integral expressed in (7).

####  2.2.1 Monte-Carlo E-step

The integrals in the Q-function do not have closed forms when
$f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k}^{(s,m)};\boldsymbol{\theta})$
is assumed to be non-Gaussian. We approximate these integrals using a Markov
Chain Monte Carlo (MCMC) sample of size M from the posterior density
$\phi(\boldsymbol{\alpha}_{k}|\boldsymbol{d}_{k,o};\boldsymbol{\theta}^{(s)})$.
Let $\boldsymbol{\alpha}_{k}^{(s,m)}$ be the $m^{th}$ simulated
$r$-dimensional vector from the posterior of the latent common factors,
$m=1,...,M$, at the $s^{th}$ iteration of the algorithm for group $k$. The
integral in (7) can be approximated as

$Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})\approx-\frac{1}{M}\sum_{m=1}^{M}\sum_{k=1}^{K}\log
f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\alpha}_{k}^{(s,m)};\boldsymbol{\theta}).$

We use the fast and efficient No-U-Turn Hamiltonian Monte Carlo sampling
procedure (NUTS) from the Stan software (Carpenter et al., 2017; Hoffman and
Gelman, 2014) in order to perform the E-step efficiently.

####  2.2.2 M-step

In the M-step of the algorithm, we aim to minimize

$Q_{1,\lambda}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})=Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})+\lambda_{0}\sum_{j=1}^{p}\rho_{0}\left(\beta_{j}\right)+\lambda_{1}\sum_{t=1}^{q}\rho_{1}\left(||\boldsymbol{b}_{t}||_{2}\right)$
(9)

with respect to
$\boldsymbol{\theta}=(\boldsymbol{\beta}^{T},\boldsymbol{b}^{T},\tau)^{T}$. We
do this by using a Majorization-Minimization algorithm with penalties applied
to the fixed effects and the rows of $\boldsymbol{B}$. Let $s$ represent the
iteration of the MCECM algorithm, and let $h$ represent the iteration within a
particular M-step of the MCECM algorithm. The M-step of the $s^{th}$ iteration
of the MCECM algorithm proceeds as described in Algorithm 1 given in Appendix
A.

####  2.2.3 MCECM algorithm

Algorithm 2 in Appendix A describes the full MCECM algorithm for estimating
the parameters with a particular $(\lambda_{0},\lambda_{1})$. The process of
model selection and finding optimal tuning parameters are described further in
the Web Appendix (see Sections 1.1 and 1.2). For further details on
initialization and convergence, also see the Web Appendix 1.3.

###  2.3 Advantages of glmmPen_FA model formulation

There are several advantages to our proposed factor model decomposition of the
random effects. By representing the random effects with a factor model, we
reduce the latent space from a high $q$-dimensional space to a low
$r$-dimensional. In the more traditional GLMM model formulation,
$Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})$ would be represented as

$Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})=-\sum_{k=1}^{K}\int\log[f(\boldsymbol{y}_{k}|\boldsymbol{X}_{k},\boldsymbol{\delta}_{k};\boldsymbol{\theta})]\phi(\boldsymbol{\delta}_{k}|\boldsymbol{d}_{k,o};\boldsymbol{\theta}^{(s)})d\boldsymbol{\delta}_{k},$
(10)

where $\boldsymbol{\theta}$ includes the fixed effects $\boldsymbol{\beta}$,
the non-zero elements of $\boldsymbol{\Gamma}$ given in (2), and $\tau$, and
the $\boldsymbol{\delta}_{k}$ are $q$-dimensional latent variables. However,
by using the novel model formulation given in equation (3), this changes the
integral of interest such that now
$Q_{1}(\boldsymbol{\theta}|\boldsymbol{\theta}^{(s)})$ expressed in (7) is of
dimension $r\ll q$. This significantly reduces the computational complexity of
estimating this integral in the E-step of the algorithm since we only have to
estimate a latent space of dimension $r$. Consequently, this reduces the
computational time. This also enables us to scale our method to hundreds of
predictors since the practical dimension of the latent space will be much
smaller than the total possible random effects predictors.

Furthermore, this proposed formulation allows for more complex correlation
structures in higher dimensions. In Rashid et al. (2020), the authors
approximated the random effect covariance matrix $\boldsymbol{\Sigma}$ as a
diagonal matrix when the dimensions are large as recommended by Fan and Li
(2012). This approximation was employed in order to reduce the computational
complexity of the algorithm and therefore increase the speed of the model fit.
However, in our new formulation, we do not need to assume
$\boldsymbol{\Sigma}$ is a diagonal matrix when the dimension is high.

###  2.4 Estimation of the number of latent factors

Performing our proposed glmmPen_FA method requires specifying the number of
latent factors $r$. Since $r$ is typically unknown _a priori_ , this value
needs to be estimated. There have been several proposed methods of estimating
$r$ for the approximate factor model (Bai and Ng, 2002; Ahn and Horenstein,
2013; Onatski, 2010; Kapetanios, 2010). We tried the Eigenvalue Ratio method
and Growth Ratio method developed by Ahn and Horenstein (2013) as well as the
method proposed in Bai and Ng (2002). We found the Growth Ratio (GR) method
gave the most accurate estimates of $r$. Therefore, in this section, we will
describe how we implement the GR method to estimate $r$. The GR method is used
in all of our numerical works.

To apply the GR method to our problem, we need a $q\times K$ matrix of
observed random effects. Since we can never observe the random effects, we
instead calculate pseudo random effects by first fitting a penalized
generalized linear model with a small penalty to each group individually. We
then take these group-specific estimates and center them so that all features
have a mean of 0. Let these $q$-dimensional group-specific estimates be
denoted as $\boldsymbol{\hat{\gamma}_{k}}$ for each group $k=1,...,K$. We then
define
$\boldsymbol{G}=(\boldsymbol{\hat{\gamma}_{1}},...,\boldsymbol{\hat{\gamma}_{K}})$
as the final $q\times K$ matrix of pseudo random effects.

Let $\psi_{j}(A)$ be the $j$-th largest eigenvalue of the positive
semidefinite matrix $A$, and let
$\tilde{\mu}_{qK,j}\equiv\psi_{j}(\boldsymbol{G}\boldsymbol{G}^{T}/(qK))=\psi_{j}(\boldsymbol{G}^{T}\boldsymbol{G}/(qK))$.

To find the GR estimator, we first order the eigenvalues of
$\boldsymbol{G}\boldsymbol{G}^{T}/(qK)$ from largest to smallest. Then, we
calculate the following ratios:

$GR(j)\equiv\frac{\log[V(j-1)/V(j)]}{\log[V(j)/V(j+1)]}=\frac{\log(1+\tilde{\mu}_{qK,j}^{*})}{\log(1+\tilde{\mu}_{qK,j+1}^{*})},$
(11)

for $j=1,2,...,U$, where $V(j)=\sum_{l=j+1}^{\min(q,K)}\tilde{\mu}_{qK,l}$,
$\tilde{\mu}_{qK,j}^{*}=\tilde{\mu}_{qK,j}/V(j)$, and $U$ is a pre-defined
constant. Then, we estimate $r$ by

$\widehat{r}_{GR}=\text{max}_{1\leq j\leq U}GR(j)$ (12)

##  3 Simulations

In this section, we examine the performance of the glmmPen_FA algorithm in
performing variable selection under several different conditions. In all of
these simulations, we use a pre-screening step to remove some random effects
at the start of the algorithm, the BIC-ICQ (Ibrahim et al., 2011) criterion
for tuning parameter selection, the MCP penalty (MCP penalty for the fixed
effects, group MCP penalty for the rows of the $\boldsymbol{B}$ matrix), and
the abbreviated two-stage grid search as described in the Web Appendix Section
1.1. In order to determine the robustness of our variable selection procedure
based on the assumed value of $r$, we fit models in one of two ways: we
estimated the number of common factors $r$ using the Growth Ratio estimation
procedure discussed in Section 2.4, or we input the true value of $r$ for the
algorithm to use.

###  3.1 Variable Selection in Binomial Data with 100 Predictors

We examine the performance of the glmmPen_FA algorithm when performing
variable selection in high dimensions of $p=100$ total predictors. We
simulated binary responses from a logistic mixed effects model with $p=100$
predictors. Of $p$ total predictors, we assume that the first 10 predictors
have truly non-zero fixed and random effects, and the other $p-10$ predictors
have zero-valued fixed and random effects. We specified a full model for the
algorithm such the random effect predictors equalled the fixed effect
predictors ( e.g. $q=p$), and our aim was to select the set of true predictors
and random effects.

To simulate the data, we set the sample size to $N=2500$ and number of groups
to $K=25$, with an equal number of subjects per group. We set up the random
effects covariance matrix by specifying a $\boldsymbol{B}$ matrix with
dimensions $(p+1)\times r$, where $p+1$ represents the $p$ predictors plus the
random intercept, and the number of latent common factors $r=\\{3,5\\}$.
Eleven of these $p+1$ rows—corresponding to the true 10 predictors plus the
intercept—had non-zero elements, while the remaining $p-10$ rows were set to
zero. For each value of $r$, we considered $\boldsymbol{B}$ matrices that
produced covariance matrices
$\boldsymbol{\Sigma}=\boldsymbol{B}\boldsymbol{B}^{T}$ with moderate variances
and eigenvalues and large variances and eigenvalues (see Appendix A for
further details). In the simulation results presented in this section, these
$\boldsymbol{B}$ matrices are referred to as the ‘moderate’ and ‘large’
$\boldsymbol{B}$ matrices, respectively. We use both moderate predictor
effects and strong predictor effects, where all 10 of the true fixed effects
have coefficient values of 1 or 2, respectively.

For group $k$, we generated the binary response $y_{ki}$, $i=1,...,n_{k}$ such
that $y_{ki}\sim Bernoulli(p_{ki})$ where
$p_{ki}=P(y_{ki}=1|\boldsymbol{x}_{ki},\boldsymbol{z}_{ki},\boldsymbol{\gamma}_{k},\boldsymbol{\theta})=\exp(\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{\gamma}_{k})/\\{1+\exp(\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{\gamma}_{k})\\}$,
and $\boldsymbol{\gamma}_{k}\sim N_{10}(0,\boldsymbol{B}\boldsymbol{B}^{T})$.
Each condition was evaluated using 100 total simulated datasets.

For individual $i$ in group $k$, the vector of predictors for the fixed
effects was $\boldsymbol{x}_{ki}=(1,x_{ki,1},...,x_{ki,p})^{T}$, and we set
the random effects $\boldsymbol{z}_{ki}=\boldsymbol{x}_{ki}$, where
$x_{ki,j}\sim N(0,1)$ for $j=1,...,p$, and each $\boldsymbol{x}_{j}$ was
standardized as described in Section 2.1.

The results for these simulations are presented in Tables 3.1 and 2. Table 3.1
provides the average true and false positives for both the fixed and random
effects variable selection, the median time in hours to complete the variable
selection procedure, and the average of the mean absolute deviation between
the coefficient estimates and the true coefficients across all simulation
replicates. Table 2 gives the Growth Ratio $r$ estimation procedure results,
including the average estimate of $r$ and the proportion of times that the
Growth Ratio estimate of $r$ was underestimated, correct, or overestimated.
All simulations were completed on a Longleaf computing cluster (CPU Intel
processors between 2.3Ghz and 2.5GHz).

True $r$ | $\beta$ | $\boldsymbol{B}$ |  $r$ Est. |  TP % |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
|  FP % |  |  | Fixef |  |  |  |  |
|  TP % |  |  |  |  |  |  |  |
|  FP % |  |  |  |  |  |  |  |
| $T^{med}$ |  Abs. Dev. |  |  |  |  |  |  |
|  | (Mean) |  |  |  |  |  |  |
3 | 1 | Moderate | GR | 98.50 | 2.00 | 97.20 | 0.22 | 2.05 | 0.26
|  |  | True | 99.00 | 2.14 | 98.40 | 0.16 | 2.36 | 0.26
|  | Large | GR | 95.50 | 2.19 | 98.60 | 0.18 | 2.52 | 0.33
|  |  | True | 95.50 | 2.31 | 98.90 | 0.17 | 2.42 | 0.33
| 2 | Moderate | GR | 100.00 | 2.46 | 89.00 | 0.53 | 1.45 | 0.37
|  |  | True | 100.00 | 2.78 | 90.10 | 0.50 | 2.07 | 0.31
|  | Large | GR | 100.00 | 3.39 | 94.60 | 0.80 | 2.39 | 0.43
|  |  | True | 100.00 | 3.40 | 96.20 | 0.49 | 2.41 | 0.41
5 | 1 | Moderate | GR | 96.80 | 2.02 | 96.20 | 0.04 | 3.56 | 0.35
|  |  | True | 96.70 | 1.86 | 96.80 | 0.03 | 3.60 | 0.35
|  | Large | GR | 90.40 | 2.22 | 96.80 | 0.08 | 4.39 | 0.44
|  |  | True | 90.50 | 1.97 | 96.90 | 0.07 | 4.44 | 0.44
| 2 | Moderate | GR | 100.00 | 2.11 | 89.00 | 0.18 | 2.29 | 0.52
|  |  | True | 100.00 | 2.42 | 88.40 | 0.24 | 2.99 | 0.44
|  | Large | GR | 99.90 | 3.28 | 93.10 | 0.50 | 3.03 | 0.57
|  |  | True | 99.90 | 3.36 | 93.40 | 0.47 | 3.98 | 0.55

Table 1: Variable selection results for the $p=100$ logistic regression
simulations, including true positive (TP) percentages for fixed and random
effects, false positive (FP) percentages for fixed and random effects, the
median time in hours for the algorithm to complete ($T^{med}$), and the
average of the mean absolute deviation (Abs. Dev. (Mean)) between the
coefficient estimates and the true $\beta$ values across all simulation
replicates. Column $\boldsymbol{B}$ describes the general size of both the
variances and eigenvalues of the resulting
$\boldsymbol{\Sigma}=\boldsymbol{BB}^{T}$ random effects covariance matrix.
Column ‘$r$ Est.’ refers to the method used to specify $r$ in the algorithm:
the Growth Ratio estimate or the true value of $r$.

True $r$ | $\beta$ | $\boldsymbol{B}$ | Avg. $r$ |  $r$ Underestimated % |  $r$ Correct % |  $r$ Overestimated %
---|---|---|---|---|---|---
3 | 1 | Moderate | 2.79 | 21 | 79 | 0
|  | Large | 2.95 | 5 | 95 | 0
| 2 | Moderate | 2.21 | 80 | 19 | 1
|  | Large | 2.51 | 49 | 51 | 0
5 | 1 | Moderate | 4.60 | 26 | 72 | 2
|  | Large | 4.83 | 15 | 83 | 2
| 2 | Moderate | 3.83 | 70 | 28 | 2
|  | Large | 4.43 | 46 | 50 | 4
Table 2: Results of the Growth Ratio $r$ estimation procedure for $p=100$
logistic mixed effects simulation results, including the average estimate of
$r$ across simulations and percent of times that the estimation procedure
underestimated $r$, gave the true $r$, or overestimated $r$. Column
$\boldsymbol{B}$ describes the general size of both the variances and
eigenvalues of the resulting $\boldsymbol{\Sigma}=\boldsymbol{BB}^{T}$ random
effects covariance matrix.

We see from Table 3.1 that the glmmPen_FA method is able to accurately select
both the fixed and random effects across a variety of conditions, which is
supported by the true positives generally being above 90% for both the fixed
and random effects and the false positives generally being small: across all
conditions, less than 3.5% for fixed effects and less than 1% for random
effects.

We can see from Table 2 that the Growth Ratio estimation procedure applied to
the pseudo random effect estimates described in Section 2.4 has varying levels
of accuracy depending on the structure of the underlying data. Generally, the
Growth Ratio estimation procedure becomes more accurate as the eigenvalues of
the covariance matrix increase and the true predictor effects are moderate.
From simulations not shown, the estimation of $r$ also generally improves when
either the sample size per group increases, or the total number of predictors
used in the GR estimation procedure decreases. When the eigenvalues of the
covariance matrix decrease or the true predictor effects increase, the Growth
Ratio procedure underestimates $r$ on average. However, when we compare the
true and false positives for the fixed and random effects given in Table 3.1
for these less ideal cases, we see that using the estimated $r$ gave very
similar results to when the true value of $r$ was utilized by the algorithm.
From the average of the mean absolute deviation values, we see that the mis-
specification of $r$ also does not significantly impact the estimation of the
fixed effects coefficients.

###  3.2 Comparison of glmmPen vs glmmPen_FA Methods

As far as we are aware, the glmmPen method developed by Rashid et al. (2020)
and implemented in the glmmPen R package available on CRAN is the only other
method that performs simultaneous fixed and random effects variable selection
in high dimensional GLMMs.

We next compare the performance of this glmmPen method and our novel
glmmPen_FA method developed in this paper. We first compared the performance
of these methods in moderate dimensions. We simulated binary responses from a
logistic mixed effects model much like the procedure described in Section 3.1,
except the total number of predictors used in the analyses was $p=25$ and we
restricted our consideration to $r=3$ common factors for all simulation
scenarios. For the glmmPen_FA method, all values of $r$ used in the algorithm
were from the Growth Ratio estimates of $r$. In these moderate dimensions of
$p=25$ with our given sample size of $N=2500$, it is reasonable to use glmmPen
to perform variable selection in logistic mixed effects models assuming an
unstructured random effects covariance matrix, allowing us to use as directly
comparable model assumptions as possible for these method comparisons.

Table 3.2 gives the average true and false positives for both the fixed and
random effects, the median time in hours to complete the variable selection
procedure, and the average of the mean absolute deviation between the
coefficient estimates and the true coefficients across all simulation
replicates. Table 4 gives the Growth Ratio $r$ estimation procedure results
for the glmmPen_FA method.

$\beta$ | $\boldsymbol{B}$ | Method |  TP % |  |  |  |  |
---|---|---|---|---|---|---|---|---
|  FP % |  | Fixef |  |  |  |  |
|  TP % |  |  |  |  |  |  |
|  FP % |  |  |  |  |  |  |
| $T^{med}$ |  Abs. Dev. |  |  |  |  |  |
|  | (Mean) |  |  |  |  |  |
1 | Moderate | glmmPen_FA | 99.50 | 3.47 | 99.80 | 0.53 | 0.59 | 0.26
|  | glmmPen | 100.00 | 37.40 | 100.00 | 0.00 | 2.62 | 0.27
| Large | glmmPen_FA | 97.20 | 3.80 | 99.90 | 0.80 | 0.49 | 0.33
|  | glmmPen | 99.00 | 60.60 | 100.00 | 0.00 | 3.18 | 0.34
2 | Moderate | glmmPen_FA | 100.00 | 2.27 | 98.40 | 0.47 | 0.47 | 0.27
|  | glmmPen | 100.00 | 13.67 | 99.20 | 0.00 | 2.84 | 0.43
| Large | glmmPen_FA | 99.80 | 3.53 | 99.80 | 0.73 | 0.48 | 0.35
|  | glmmPen | 100.00 | 30.93 | 100.00 | 0.00 | 2.53 | 0.50

Table 3: Results of the variable selection procedure for the $p=100$ logistics
mixed effects simulations, including true positive (TP) percentages for fixed
and random effects, false positive (FP) percentages for fixed and random
effects, the median time in hours for the algorithm to complete ($T^{med}$),
and the average of the mean absolute deviation (Abs. Dev. (Mean)) between the
coefficient estimates and the true $\beta$ values across all simulation
replicates. Column $\boldsymbol{B}$ describes the general size of both the
variances and eigenvalues of the resulting
$\boldsymbol{\Sigma}=\boldsymbol{BB}^{T}$ random effects covariance matrix.
All values of $r$ used in the glmmPen_FA method were from the Growth Ratio
estimates of $r$.

True $r$ | $\beta$ | $\boldsymbol{B}$ | Avg. $r$ |  $r$ Underestimated % |  $r$ Correct % |  $r$ Overestimated %
---|---|---|---|---|---|---
3 | 1 | Moderate | 3.00 | 0 | 100 | 0
|  | Large | 3.00 | 0 | 100 | 0
| 2 | Moderate | 2.76 | 24 | 76 | 0
|  | Large | 2.92 | 8 | 92 | 0
Table 4: Results of the Growth Ratio $r$ estimation procedure for glmmPen_FA
$p=25$ logistic mixed effects simulation results, including the average
estimate of $r$ across simulations and percent of times that the estimation
procedure underestimated $r$, gave the true $r$, or overestimated $r$. Column
$\boldsymbol{B}$ describes the general size of both the variances and
eigenvalues of the resulting $\boldsymbol{\Sigma}=\boldsymbol{BB}^{T}$ random
effects covariance matrix.

When comparing the glmmPen_FA and glmmPen results in these $p=25$ simulations,
we see that the median time for glmmPen to complete the variable selection
procedures ranged from 2.53 to 3.28 hours for all four simulation scenarios
considered. On the other hand, the glmmPen_FA method was able to fit these
variable selection procedures about 5-6 times faster, where the median running
time ranged from 0.47 to 0.59 hours.

Table 3.2 also shows that there is little difference in the true positives for
both the fixed and random effects between the two methods. However, glmmPen
tends to have more false positives in the fixed effects.

We also performed variable selection using glmmPen on the $r=3$, $p=100$
simulations described in Section 3.1. In these larger dimensions, we
simplified the glmmPen estimation procedure by assuming an independent
covariance matrix to reduce the number of random effects covariance
parameters. We let the glmmPen variable selection procedure proceed for 100
hours. In that time, glmmPen was able to complete the following number of
replicates out of the 100 total replicates: 83 for ($\beta=1,\boldsymbol{B}=$
Moderate), 71 for ($\beta=1,\boldsymbol{B}=$ Large), 100 for
($\beta=2,\boldsymbol{B}=$ Moderate), and 96 for for
($\beta=2,\boldsymbol{B}=$ Large). The minimum times needed to complete the
glmmPen variable selection procedures: 39.91, 57.60, 23.63, and 42.79 hours,
respectively; in summary, it took a day or more for the fastest simulation
replicates to complete when using the glmmPen method. In cases where we desire
to select true random effects from a large number of total predictors, it is
clear that the glmmPen_FA estimation procedure significantly reduces the
required time to perform variable selection.

###  3.3 Variable Selection in Binomial Data with 500 Predictors

In order to further illustrate the scalability of our method, we applied our
method to binary outcome simulations with $p=500$ covariates. We simulated the
binary responses from a logistic mixed effects model much like the procedure
described in Section 3.1, except the total number of predictors used in the
analyses was $p=500$ instead of $p=100$. All simulations assumed the true
number of latent factors $r$ was 3 and the Growth Ratio method was used to
estimate $r$. Just as in the $p=100$ binary outcome simulations, we specified
a full model for the algorithm such the random effect predictors equalled the
fixed effect predictors ( e.g. $q=p$), and our aim was to select the set of
true predictors and random effects. The variable selection results to these
simulations are given in Table 3.3. The median times needed to complete these
simulations took between 10.37 and 17.57 hours.

True $r$ | $\beta$ | $\boldsymbol{B}$ |  $r$ Avg. |  TP % |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
|  FP % |  |  | Fixef |  |  |  |  |
|  TP % |  |  |  |  |  |  |  |
|  FP % |  |  |  |  |  |  |  |
| $T^{med}$ |  Abs. Dev. |  |  |  |  |  |  |
|  | (Mean) |  |  |  |  |  |  |
3 | 1 | Moderate | 2.67 | 97.40 | 0.97 | 93.90 | 0.05 | 10.37 | 0.27
|  | Large | 2.85 | 88.50 | 1.73 | 94.30 | 0.33 | 17.57 | 0.37
| 2 | Moderate | 2.37 | 100.00 | 0.06 | 77.40 | 0.00 | 11.44 | 0.48
|  | Large | 2.41 | 99.60 | 0.19 | 88.10 | 0.03 | 13.22 | 0.55

Table 5: Results of the variable selection procedure for $p=500$ logistic
mixed effects simulation results, including true positive (TP) percentages for
fixed and random effects, false positive (FP) percentages for fixed and random
effects, the median time in hours for the algorithm to complete ($T^{med}$),
and the average of the mean absolute deviation (Abs. Dev. (Mean)) between the
coefficient estimates and the true $\beta$ values across all simulation
replicates. Column ‘$r$ Avg.’ gives the average Growth Ratio $r$ estimate used
within the algorithm. Column $\boldsymbol{B}$ describes the general size of
both the variances and eigenvalues of the resulting
$\boldsymbol{\Sigma}=\boldsymbol{BB}^{T}$ random effects covariance matrix.

###  3.4 Variable Selection in Poisson Data with 100 Predictors

While we have previously focused on binary outcome data in our simulations,
our proposed method also applies to other members of the generalized linear
model family, including the Poisson model for count outcome data. To
illustrate this, we simulated a Poisson mixed effects model with $p=100$
predictors, 5 of which had truly non-zero fixed and random effects, and the
other $p-5$ predictors had zero-valued fixed and random effects. As in the
previous Binomial simulations, we set the sample size to $N=2500$ and the
number of groups to $K=25$, with equal numbers of subjects per group. We set
$r$ to 3, assigned moderate predictor effects ($\beta=1$), and specified a
$\boldsymbol{B}$ matrix with 6 non-zero rows (for the 5 predictors plus an
intercept) that produced a ‘moderate’ covariance matrix (see the Appendix A
for details). Unlike in the previous Binomial simulations, we simulated
$x_{ki,j}\sim N(0,\sigma=0.10)$ for $j=1,...,p$ to reduce the overall spread
of the simulated $y_{ki}\sim Poisson(\mu_{ki})$ outcome values, where
$\mu_{ki}=\exp(\boldsymbol{x}_{ki}^{T}\boldsymbol{\beta}+\boldsymbol{z}_{ki}^{T}\boldsymbol{\gamma}_{k})$.

Using the Growth Ratio estimation procedure to estimate $r$, the average true
positive percentages were 83.40% for the fixed effects and for 77.20% the
random effects, and the average false positives were 5.92% for the fixed
effects and 3.23% for the random effects. The average estimate of $r$ across
the simulation replicates was 2.35.

##  4 Case study: Pancreatic Ductal Adenocarcinoma

Patients diagnosed with Pancreatic Ductal Adenocarcinoma (PDAC) generally face
a very poor prognosis, where the 5-year survival rate is 6% (Khorana et al.,
2016). The study by Moffitt et al. (2015) identified genes that are expressed
exclusively in pancreatic tumor cells. Using these tumor-specific genes,
Moffitt et al. (2015) was able to identify and validate two novel tumor
subtypes, termed ‘basal-like’ and ‘classical’. It was found that patients
diagnosed with basal-like tumors had significantly worse median survival than
those diagnosed with the classical tumors. Consequently, it is of clinical
interest to robustly predict this basal-like subtype in order to make and
improve tailored treatment recommendations.

A common problem observed in the application of gene signatures is the
inconsistency of gene signature selection across biomedical studies, where
gene signatures identified in one study may have little or no overlap with
ones identified in other studies (Waldron et al., 2014). This lack of
replicability of gene signatures also translates to a lack of replicability of
other findings, where models based upon these different gene signatures have
variable accuracy in predicting clinical outcomes in new studies (Sotiriou and
Piccart, 2007; Waldron et al., 2014) or provide contradictory effect estimates
relating genes to the outcome (Swisher et al., 2012). This lack of
replicability across studies comes from several sources, including small
sample size (Sotiriou and Piccart, 2007), low frequencies of one or more
molecular subtypes (Lusa et al., 2007), and differences in data pre-processing
steps (Lusa et al., 2007; Paquet and Hallett, 2015).

In order to improve replicability in the prediction of subtypes in PDAC, we
combine PDAC gene expression data from five different studies. The studies
used in these analyses are summarized in Web Appendix Table 3. In order to
account and adjust for between-study heterogeneity, we apply our new method
glmmPen_FA to fit a penalized logistic mixed effects model to our data to
select predictors with study-replicable effects, where we assume that
predictor effects may vary between studies.

Moffitt et al. (2015) identified a 500-member gene list relevant to
classifying the PDAC tumor subtypes. Of these 500 genes, the five studies
described in Web Appendix Table 3 had RNA-seq gene expression data for 432 of
these genes. There were some significant correlations between some of these
432 genes, as evaluated by Spearman correlations. In order to avoid having
very highly correlated covariates in the analyses, we decided to combine
highly correlated genes together into meta-genes. The clustering process used
to create these meta-genes is described in Web Appendix Section 1.6. The final
dataset included 110 cluster covariates. If clusters were composed of two or
more genes, the the raw RNA-seq gene expression for these genes was summed
together to create a meta-gene. The final covariates used in the analyses were
subject-level rank transformations of the gene expression corresponding to
these clusters.

Due to the presence of several pairwise Spearman correlation values greater
than 0.5 in this final dataset, we used the Elastic Net penalization procedure
(Friedman et al., 2010) to balance between ridge regression and the MCP
penalty. We let $\pi$ represent the balance between ridge regression and the
MCP penalty, where $\pi=0$ represents ridge regression and $\pi=1$ represents
the MCP penalty. In these analyses, we let $\pi=0.1$, and we estimated the
number of common factors $r$ using the Growth Ratio procedure. The same value
of $\pi$ was used for both the fixed effects and random effects penalization.
The sequence of $\lambda$ penalties used the the variable selection procedure
was the same as those used in the Binomial variable selection simulations for
$p=100$ as discussed in Web Appendix Section 1.2.

To see the general performance of glmmPen_FA under different covariate
correlation structures and under different values of $\pi$, see Web Appendix
Section 1.5.

The basal or classical subtype outcome was calculated using the clustering
algorithm specified in Moffitt et al. (2015). Further details are provided in
Web Appendix Section 1.6; the code for this procedure is provided in a GitHub
repository, see Supporting Information for more details.

In the final results, 8 of the original 110 cluster covariates had non-zero
fixed effect values in the best model, implying these covariates were
important for the prediction of the basal outcome. These 8 cluster covariates
represented 33 genes in total. Table 4 includes the cluster label for these 8
clusters, the sign of the associated fixed effect coefficient, and the gene
symbols of the genes that make up the cluster. Clusters with positive
coefficients indicate that having greater relative expression of these meta-
genes increases the odds of a subject being in the basal subtype, and vice
versa for negative coefficients. The best model contained a random intercept
(variance value 1.57) and no other random slopes. The Growth Ratio procedure
estimated $r=2$ latent common factors to model the random effects.

Cluster |  |
---|---|---
No. |  Coefficient |
| Gene Component Symbols |
21 | + | C16orf74, HES2, S100A2
25 | + | CCNG2, MBOAT2, MET, TGFB2
44 | + | FAM83A, MUC16, SCEL
45 | + | FANK1, HS3ST1, LEMD1
58 | + | KRT15, PAK6, PROM2, ST3GAL4, TRIM29
64 | + | MMP13, PMAIP1
70 | + | PSCA, VGLL1, WNT10A
91 | - | ANXA10, BTNL8, CLDN18, CYP2C18, GUCY2C, LRRC31,
|  | MYO1A, NR1I2, PIP5K1B, REG4

Table 6: Covariate cluster label within the case study dataset of the clusters
that had non-zero fixed effects in the final best model, the sign of the fixed
effect coefficient associated with the cluster, and the gene symbols of the
genes within the cluster.

Sensitivity analyses were conducted, which included additional values for
$\pi$ in the Elastic Net procedure and a larger value of the number of assumed
latent common factors $r$. Details are included in Web Appendix Section 1.7.

We also applied the glmmPen variable selection procedure to this data
(assuming an independent random effects covariance matrix). The 8 cluster
covariates selected by glmmPen_FA were also consistently selected by glmmPen.
One difference in the model results is that glmmPen also consistently selected
clusters 25 and 58 to have non-zero random effects. The main difference in
these variable selection procedures were the time to complete the variable
selection procedure, where glmmPen_FA completed the procedure within 0.4 hours
and glmmPen completed the procedure within 37.8 hours. More details about
glmmPen sensitivity analyses are provided in Web Appendix Section 1.7.

##  5 Discussion

By representing the random effects with a factor model, we reduce the latent
space from a large number of random effects to a smaller set of latent
factors. We have shown through simulations that by reducing the complexity of
the integral in the E-step, we can significantly improve the overall time
needed to perform variable selection in high dimensional generalized linear
mixed models.

The simulations in Section 3 also show how reducing the latent space increases
the feasible dimensionality of performing variable selection in generalized
linear mixed models. By using our novel formulation of the random effects, we
can perform variable selection on mixed models with hundreds of predictors
within a reasonable time-frame without any _a priori_ knowledge of which
predictors are relevant for the model, either in terms of fixed or random
effects. From the simulation results, we see that the glmmPen_FA method
results in accurate selection of the fixed and random effects across several
conditions.

In general, our method is limited by the need to provide an estimate for the
number of latent common factors. However, this limitation is tempered by
several observations from the simulation results. The simulation results show
that our data-driven estimation of the number of latent factors, based on the
Growth Ratio estimation procedure by Ahn and Horenstein (2013), provides
reasonable estimates. Even when it was estimated incorrectly by this
procedure, this mis-specification had very little impact on the general
variable selection performance or the fixed effects coefficient estimates.
Therefore, our method is not sensitive to the estimation of the number of
latent factors.

## Data Availability Statement

The data used for the case study analysis in Section 4 is provided as an RData
object in the Data folder of the following GitHub repository:
https://github.com/hheiling/paper_glmmPen_FA. The code used to create this
cleaned dataset from the individual studies is provided in the Code_CaseStudy
folder. The individual studies are also provided.

## References

  * Ahn and Horenstein (2013) Ahn, S. C. and Horenstein, A. R. (2013).  Eigenvalue ratio test for the number of factors.  Econometrica 81, 1203–1227.
  * Bai and Ng (2002) Bai, J. and Ng, S. (2002).  Determining the number of factors in approximate factor models.  Econometrica 70, 191–221.
  * 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, 1069–1077.
  * Breheny and Huang (2011) Breheny, P. and Huang, J. (2011).  Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection.  Annals of Applied Statistics 5, 232–253.
  * Breheny and Huang (2015) Breheny, P. and Huang, J. (2015).  Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors.  Statistics and computing 25, 173–187.
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of statistical software 76,.
  * Chen and Dunson (2003) Chen, Z. and Dunson, D. B. (2003).  Random effects selection in linear mixed models.  Biometrics 59, 762–769.
  * Fan et al. (2013) Fan, J., Liao, Y., and Mincheva, M. (2013).  Large covariance estimation by thresholding principal orthogonal complements.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75, 603–680.
  * Fan and Li (2012) Fan, Y. and Li, R. (2012).  Variable selection in linear mixed effects models.  Annals of statistics 40, 2043\.
  * Friedman et al. (2010) Friedman, J., Hastie, T., and Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of statistical software 33, 1–22.
  * Garcia et al. (2010) Garcia, R. I., Ibrahim, J. G., and Zhu, H. (2010).  Variable selection for regression models with missing data.  Statistica Sinica 20, 149\.
  * Gurka et al. (2011) Gurka, M. J., Edwards, L. J., and Muller, K. E. (2011).  Avoiding bias in mixed model inference for fixed effects.  Statistics in medicine 30, 2696–2707.
  * Heiling et al. (2023) Heiling, H., Rashid, N., Li, Q., and Ibrahim, J. (2023).  glmmPen: High Dimensional Penalized Generalized Linear Mixed Models (pGLMM).  R package version 1.5.3.4.
  * Hoffman and Gelman (2014) Hoffman, M. D. and Gelman, A. (2014).  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  J. Mach. Learn. Res. 15, 1593–1623.
  * 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, 495–503.
  * Kapetanios (2010) Kapetanios, G. (2010).  A testing procedure for determining the number of factors in approximate factor models with large datasets.  Journal of Business & Economic Statistics 28, 397–409.
  * Khorana et al. (2016) Khorana, A. A., Mangu, P. B., Berlin, J., Engebretson, A., Hong, T. S., Maitra, A., Mohile, S. G., Mumber, M., Schulick, R., Shapiro, M., et al. (2016).  Potentially curable pancreatic cancer: American society of clinical oncology clinical practice guideline.  Journal of Clinical Oncology 34, 2541–2556.
  * Lusa et al. (2007) Lusa, L., McShane, L. M., Reid, J. F., De Cecco, L., Ambrogi, F., Biganzoli, E., Gariboldi, M., and Pierotti, M. A. (2007).  Challenges in projecting clustering results across gene expression–profiling datasets.  JNCI: Journal of the National Cancer Institute 99, 1715–1723.
  * Moffitt et al. (2015) Moffitt, R. A., Marayati, R., Flate, E. L., Volmar, K. E., Loeza, S. G. H., Hoadley, K. A., Rashid, N. U., Williams, L. A., Eaton, S. C., Chung, A. H., et al. (2015).  Virtual microdissection identifies distinct tumor-and stroma-specific subtypes of pancreatic ductal adenocarcinoma.  Nature genetics 47, 1168\.
  * Onatski (2010) Onatski, A. (2010).  Determining the number of factors from empirical distribution of eigenvalues.  The Review of Economics and Statistics 92, 1004–1016.
  * Paquet and Hallett (2015) Paquet, E. R. and Hallett, M. T. (2015).  Absolute assignment of breast cancer intrinsic molecular subtype.  Journal of the National Cancer Institute 107, dju357.
  * Rashid et al. (2020) Rashid, N. U., Li, Q., Yeh, J. J., and Ibrahim, J. G. (2020).  Modeling between-study heterogeneity for improved replicability in gene signature selection and clinical prediction.  Journal of the American Statistical Association 115, 1125–1138.
  * Sotiriou and Piccart (2007) Sotiriou, C. and Piccart, M. J. (2007).  Taking gene-expression profiling to the clinic: when will molecular signatures become relevant to patient care?  Nature reviews cancer 7, 545–553.
  * Swisher et al. (2012) Swisher, E. M., Taniguchi, T., and Karlan, B. Y. (2012).  Molecular scores to predict ovarian cancer outcomes: a worthy goal, but not ready for prime time.
  * Thompson et al. (2017) Thompson, J. A., Fielding, K. L., Davey, C., Aiken, A. M., Hargreaves, J. R., and Hayes, R. J. (2017).  Bias and inference from misspecified mixed-effect models in stepped wedge trial analysis.  Statistics in medicine 36, 3670–3682.
  * Tran et al. (2020) Tran, M.-N., Nguyen, N., Nott, D., and Kohn, R. (2020).  Bayesian deep net glm and glmm.  Journal of Computational and Graphical Statistics 29, 97–113.
  * Waldron et al. (2014) Waldron, L., Haibe-Kains, B., Culhane, A. C., Riester, M., Ding, J., Wang, X. V., Ahmadifar, M., Tyekucheva, S., Bernau, C., Risch, T., et al. (2014).  Comparative meta-analysis of prognostic gene signatures for late-stage ovarian cancer.  JNCI: Journal of the National Cancer Institute 106,.

## Supporting Information

Web Appendices and Tables referenced in Sections 2, 3, and 4 are available
with this paper. The Web Appendix contains addition details about the variable
selection procedure, the simulations, and the case study. The glmmPen R
package, which contains the code for the original glmmPen formulation and the
new glmmPen_FA method, is available for download through CRAN at
https://cran.r-project.org/web/packages/glmmPen/index.html. The GitHub
repository https://github.com/hheiling/paper_glmmPen_FA contains the following
materials: the code to run the simulations, the code to create the cleaned
data for the case study, the code to analyze the case study data, the data
used in the case study, and the code to recreate the tables provided within
this paper. Please see the README files within this repository for more
details.

##  Appendix A A

###  A.1 B matrices used in simulations

The transpose of the first 11 rows of the deterministic ‘large’
$\boldsymbol{B}$ matrices used in the Binomial simulations in Section 3 are
given in equations (13) and (14), corresponding to $r=\\{3,5\\}$,
respectively. The deterministic ‘moderate’ $\boldsymbol{B}$ matrices are these
large $\boldsymbol{B}$ matrices multiplied by the constants 0.75 and 0.80 for
$r$ equal to 3 and 5, respectively. All other $p-10$ rows of the
$\boldsymbol{B}$ matrices were set to 0, where $p$ is the total number of
predictors used in the simulations.

$\scriptsize\boldsymbol{B}^{T}_{large,r=3}=\setcounter{MaxMatrixCols}{11}\begin{bmatrix}1&1&1&1&1&1&1&1&1&1&1\\\
0&-1&-1&-1&-1&-1&1&1&1&1&1\\\ -2&2&-1&0&1&-1&0&1&-1&0&1\\\ \end{bmatrix}$ (13)

$\scriptsize\boldsymbol{B}^{T}_{large,r=5}=\setcounter{MaxMatrixCols}{11}\begin{bmatrix}1&1&1&1&1&1&1&1&1&1&1\\\
0&-1&-1&-1&-1&-1&1&1&1&1&1\\\ -2&2&-1&0&1&-1&0&1&-1&0&1\\\
-1&1&1&-1&-1&1&1&-1&-1&1&-1\\\ -1&-1&0&1&1&-1&-1&0&1&1&-2\end{bmatrix}$ (14)

The transpose of the first 6 rows of the deterministic ‘moderate’
$\boldsymbol{B}$ matrix used in the Poisson simulations in Section 3.4 are
given in equation (15). All other $p-5$ rows of the $\boldsymbol{B}$ matrices
were set to 0, where $p=100$ in the Poisson simulations.

$\boldsymbol{B}^{T}_{poisson,r=3}=0.75\times\begin{bmatrix}1&1&1&1&1&1\\\
-1&-1&-1&1&1&1\\\ -1&0&1&-1&0&1\\\ \end{bmatrix}$ (15)

###  A.2 Algorithms

Algorithm 1 M-step of the MCECM algorithm

1\. The parameters $\boldsymbol{\theta}^{(s,0)}$ for M-step iteration $h=0$
are initialized using the results from the previous M-step,
$\boldsymbol{\theta}^{(s-1)}$.

2\. Conditional on $\boldsymbol{b}^{(s,h-1)}$ and $\tau^{(s-1)}$, each
$\beta_{j}^{(s,h)}$ for $j=1,...,p$ is given a single update using the
Majorization-Minimization algorithm specified by Breheny and Huang (2015).

3\. For each group $k$ in $k=1,...,K$, the augmented matrix
$\boldsymbol{\tilde{z}}_{ki}=(\boldsymbol{\tilde{\alpha}}_{k}^{(s)}\otimes\boldsymbol{z}_{ki})J$
is created for $i=1,...,n_{k}$ where
$\boldsymbol{\tilde{\alpha}}_{k}^{(s)}=((\boldsymbol{\alpha}_{k}^{(s,1)})^{T},...,(\boldsymbol{\alpha}_{k}^{(s,M)})^{T})^{T}$.

4\. Conditional on the $\tau^{(s-1)}$ and the recently updated
$\boldsymbol{\beta}^{(s,h+1)}$, each $\boldsymbol{b}_{t}^{(s,h)}$ for
$t=1,...,q$ is updated using the Majorization-Minimzation coordinate descent
grouped variable selection algorithm specified by Breheny and Huang (2015).

5\. Steps 2 through 4 are repeated until the M-step convergence criteria are
reached or until the M-step reaches its maximum number of iterations.

6\. Conditioning on the newly updated $\boldsymbol{\beta}^{(s)}$ and
$\boldsymbol{b}^{(s)}$, $\tau^{(s)}$ is updated (generically, using the
Newton-Raphson algorithm).

Algorithm 2 Full MCECM algorithm for single $(\lambda_{0},\lambda_{1})$
penalty combination

1\. Fixed and random effects $\boldsymbol{\beta}^{(0)}$ and
$\boldsymbol{b}^{(0)}$ are initialized as discussed the Web Appendix.

2\. In each E-step for EM iteration $s$, a burn-in sample from the posterior
distribution of the random effects is run and discarded. A sample of size
$M^{(s)}$ from the posterior is then drawn and retained for the M-step.

3\. Parameter estimates of $\boldsymbol{\beta}^{(s)}$, $\boldsymbol{b}^{(s)}$,
and $\tau^{(s)}$ are then updated as described in the M-step procedure given
above.

4\. Steps 2 and 3 are repeated until the convergence condition is met a pre-
specified consecutive number of times or until the maximum number of EM
iterations is reached.