Penalized Linear Mixed Models

• Penalized Linear Mixed Models (PLMMs) are high-dimensional regression models that blend sparse fixed effects with mixed-model correlation adjustments.
• They use preconditioning and penalties like lasso, SCAD, and MCP to enhance variable selection and address biases from correlated observations.
• PLMMs are applied in genetics, clinical prediction, and batch effect correction, with robust software implementations for large-scale data analysis.

Searching arXiv for recent and foundational papers on penalized linear mixed models and related extensions. Searching arXiv for recent work on PLMM software, cross-validation, and scalable mixed-model selection. Penalized linear mixed models (PLMMs) are mixed-effects regression models for high-dimensional data that combine sparse penalization of fixed effects with explicit modeling of correlation among observations. In the canonical Gaussian formulation, $y=X\beta+u+\epsilon$ with $u\sim N(0,\sigma_s^2K)$, $\epsilon\sim N(0,\sigma_\epsilon^2I)$, and $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$; penalization acts on the fixed effects, while $K$ captures relatedness, population structure, batch effects, or other sample-level dependence (Peter et al., 3 Feb 2025). The term now covers a family of methods rather than a single estimator: penalized profile-likelihood approaches, grouped penalties for random-effect structure, rotation-based GWAS solvers, and extensions to binary, multivariate, longitudinal, functional, and generalized mixed-model settings (Fan et al., 2012, Müller et al., 2013).

1. Statistical role and core formulation

PLMMs arise when two features of modern data occur simultaneously: the predictor dimension is large, and the observations are correlated. In GWAS and related omics settings, the sources of dependence explicitly emphasized in the literature include family relationships, cryptic relatedness, population structure or ancestry differences, geographic recruitment effects, and batch effects (Peter et al., 3 Feb 2025). Standard sparse regression methods such as the lasso assume independent observations; if that assumption fails, correlation structure can be mistaken for signal, producing biased estimates or spurious selection.

A standard mixed-model representation writes

$$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$

In GWAS-oriented PLMMs, a common choice is

$K=\frac{1}{p}XX^\top,$

the genomic relatedness matrix computed from the standardized design matrix (Peter et al., 3 Feb 2025). Marginally,

$$y\sim N(X\beta,\Sigma), \qquad \Sigma=\sigma_s^2K+\sigma_\epsilon^2I.$$

A central computational idea is preconditioning or decorrelation. Multiplying by $\Sigma^{-1/2}$ yields

$$\tilde y=\Sigma^{-1/2}y,\qquad \tilde X=\Sigma^{-1/2}X,$$

so that

$u\sim N(0,\sigma_s^2K)$0

On the transformed scale, the errors are independent and standard penalized regression machinery becomes applicable (Peter et al., 3 Feb 2025). This reduction is one of the main reasons PLMMs became attractive in structured high-dimensional genetics.

The same literature also makes clear that PLMMs are not limited to genetics. The motivating use cases explicitly include predicting clinical outcomes, plant and animal breeding, polygenic risk score construction, and Mendelian randomization contexts (Peter et al., 3 Feb 2025). A plausible implication is that PLMMs are best viewed as a general sparse correlated-regression framework whose GWAS formulations happen to be especially mature.

2. Penalization regimes and selection targets

The fixed-effect component is typically estimated by minimizing a penalized quadratic or likelihood-based criterion. In preconditioned Gaussian formulations, the generic objective is

$u\sim N(0,\sigma_s^2K)$1

with the lasso as the special case

$u\sim N(0,\sigma_s^2K)$2

The penalty repertoire explicitly documented across the literature includes lasso, SCAD, and MCP, with unpenalized covariates such as age, sex, and recruitment site often retained in the fixed-effects design but excluded from shrinkage (Peter et al., 3 Feb 2025, Ghosh et al., 2016).

PLMMs differ substantially in what they penalize. Some methods regularize only the fixed effects while estimating variance parameters jointly. In "Non-Concave Penalization in Linear Mixed-Effects Models and Regularized Selection of Fixed Effects" (Ghosh et al., 2016), the objective is the penalized negative marginal log-likelihood

$u\sim N(0,\sigma_s^2K)$3

where $u\sim N(0,\sigma_s^2K)$4 denotes the variance parameters. This formulation assumes the random-effects structure is low-dimensional and pre-specified, and it places the regularization burden entirely on $u\sim N(0,\sigma_s^2K)$5.

Other methods target both fixed and random effects. Fan and Li formulate fixed-effect selection through a penalized profile likelihood and random-effect selection through a group variable selection criterion in which the unknown random-effects covariance is replaced by a proxy matrix (Fan et al., 2012). In their framework, fixed effects are penalized coordinatewise, whereas random effects are penalized at the group level through norms of subject-stacked random-effect components. The 2013 review "Model Selection in Linear Mixed Models" (Müller et al., 2013) synthesizes related approaches that penalize Cholesky or alternative Cholesky parameterizations, because arbitrary elementwise penalization of covariance matrices can violate positive semidefiniteness.

More recent work extends the selection target further. The glmmPen framework penalizes fixed effects individually and penalizes rows of a covariance-factor matrix groupwise, thereby encoding random-effect inclusion or exclusion through grouped penalties on covariance proxies (Heiling et al., 2023). "Scalable Subset Selection in Linear Mixed Models" (Thompson et al., 25 Jun 2025) goes further by replacing continuous shrinkage with a hierarchical $u\sim N(0,\sigma_s^2K)$6 criterion,

$u\sim N(0,\sigma_s^2K)$7

subject to

$u\sim N(0,\sigma_s^2K)$8

so that a random effect can only appear if the corresponding fixed effect is present. This suggests that the PLMM literature now includes both shrinkage-based and exact-subset-selection formulations.

A recurrent theme across these formulations is that random-effect selection is structurally grouped. The review literature is explicit that selecting a random effect usually corresponds to removing a whole row of a Cholesky or covariance-factor matrix rather than thresholding isolated covariance entries (Müller et al., 2013). That distinction is central to the mixed-model version of sparsity.

3. Estimation and optimization strategies

The oldest high-dimensional PLMM formulations are likelihood-based. Fan and Li derive fixed-effect selection from a penalized profile likelihood in which the mixed-model weight matrix

$u\sim N(0,\sigma_s^2K)$9

is replaced by a proxy-based matrix

$\epsilon\sim N(0,\sigma_\epsilon^2I)$0

and they prove model-selection consistency for both fixed and random effects under conditions on the proxy matrix (Fan et al., 2012). This line of work emphasizes theoretical support recovery while avoiding direct estimation of a high-dimensional random-effects covariance.

Joint likelihood optimization with nonconcave penalties is developed in (Ghosh et al., 2016). That paper combines local quadratic approximation for SCAD-type penalties with Newton-style updates in moderate dimension and a coordinate gradient descent algorithm in high dimension. It proves low-dimensional oracle properties and high-dimensional selection consistency under assumptions including $\epsilon\sim N(0,\sigma_\epsilon^2I)$1, while also stressing that the objective is nonconvex and the algorithm can converge to local rather than global optima.

Preconditioning-based GWAS methods separate the problem into covariance estimation and sparse regression. In plmmr, the plmm() workflow consists of a prep step that eigendecomposes $\epsilon\sim N(0,\sigma_\epsilon^2I)$2, a fit step that uses coordinate descent on the transformed problem, and a format step that returns coefficients on the original scale (Peter et al., 3 Feb 2025). If $\epsilon\sim N(0,\sigma_\epsilon^2I)$3, then

$\epsilon\sim N(0,\sigma_\epsilon^2I)$4

which makes the decorrelation step computationally tractable once the eigendecomposition is available.

When the likelihood or score is not available in closed form, PLMM-relevant optimization becomes stochastic. "Stochastic Proximal Gradient Algorithms for Penalized Mixed Models" (Fort et al., 2017) studies penalized maximum likelihood problems in latent-variable models and develops MCPG and SAPG recursions based on proximal-gradient updates with Monte Carlo approximations of latent expectations. In the paper’s linear mixed-effects toy example, the fixed-effect coefficients are lasso-penalized while the latent random effects are integrated through Gaussian expectations. The main algorithmic point is that SAPG uses a stochastic approximation

$\epsilon\sim N(0,\sigma_\epsilon^2I)$5

which smooths Monte Carlo noise across iterations; with fixed batch size, the paper shows SAPG’s approximation error vanishes while MCPG’s does not.

For generalized mixed models, glmmPen uses a Monte Carlo Expectation Conditional Minimization algorithm with grouped MM coordinate descent for random-effect blocks and MCMC-based posterior sampling of random effects (Heiling et al., 2023). The latent-factor extension (Heiling et al., 2023) replaces a high-dimensional random-effects covariance by a low-rank factorization $\epsilon\sim N(0,\sigma_\epsilon^2I)$6, reducing the latent dimension from $\epsilon\sim N(0,\sigma_\epsilon^2I)$7 to $\epsilon\sim N(0,\sigma_\epsilon^2I)$8. That reformulation is broader than classical Gaussian PLMMs, but it is directly relevant to high-dimensional random-slope selection in mixed models.

Exact subset-selection methods use different computational ideas. The $\epsilon\sim N(0,\sigma_\epsilon^2I)$9 mixed-model procedure of (Thompson et al., 25 Jun 2025) applies block coordinate descent to pairs $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$0, majorizes the smooth likelihood by a quadratic upper bound, and uses a threshold map together with local search. Sherman–Morrison updates and the matrix determinant lemma reduce the cost of changing one random-effect variance at a time from $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$1 to $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$2 per cluster. The paper proves monotone descent and finite active-set stabilization, but not global optimality.

4. Correlation adjustment, prediction, and cross-validation

PLMMs are often used not only for selection but also for prediction. In mixed-model form, prediction naturally incorporates BLUP adjustments. With training data $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$3, test design $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$4, and block covariance

$\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$5

the recommended validation prediction is

$\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$6

which is the BLUP-adjusted predictor emphasized in the recent cross-validation literature (Peter et al., 18 Mar 2025). In GWAS-oriented software, BLUP-based prediction is also exposed directly through predict() and cross-validation wrappers (Peter et al., 3 Feb 2025).

The cross-validation problem is now recognized as methodologically nontrivial. "Cross-Validation in Penalized Linear Mixed Models: Addressing Common Implementation Pitfalls" (Peter et al., 18 Mar 2025) argues that the entire fitting pipeline must be nested within each fold: standardization, covariance estimation, eigendecomposition, preconditioning, penalized fitting, and BLUP prediction. The paper identifies four implementation pitfalls: constructing the preconditioner incorrectly, failing to re-standardize within each fold, computing validation BLUPs on an inconsistent scale, and performing rotation outside cross-validation.

Its most important conceptual claim is that outer rotation is invalid because preconditioning destroys exchangeability. The paper compares full, inner, and outer CV. In the small-signal simulation, full CV achieved TDR $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$7, FDR $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$8, NVAR $\Sigma=\sigma_s^2K+\sigma_\epsilon^2I$9, and RSEE $K$0; inner CV had TDR $K$1, FDR $K$2, NVAR $K$3, and RSEE $K$4; outer CV selected $K$5 variables with FDR $K$6 and RSEE $K$7 (Peter et al., 18 Mar 2025). In the PennCath real-data example, the same paper reports that outer CV selected $K$8 variables at $K$9, whereas full CV selected $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$0. This is one of the clearest demonstrations that ordinary lasso-style CV logic does not transfer unchanged to PLMMs.

A common misconception is that covariance preconditioning can be performed once on the full dataset and then treated as harmless preprocessing. The recent CV literature rejects that view explicitly. Another misconception is that BLUP prediction can reuse covariance blocks from a globally standardized matrix even when fitting used fold-specific scaling; the same paper shows that this inconsistency inflates estimation error (Peter et al., 18 Mar 2025).

5. Major extensions and neighboring model classes

The Gaussian single-trait PLMM is only one branch of the broader literature. Several nearby models preserve the PLMM logic—sparse fixed effects plus mixed-model dependence correction—while changing the outcome type, covariance structure, or predictor type.

Setting	Core modification	Representative paper
Binary outcomes	Probit orthant-likelihood or penalized GLMM/PQL	(Mandt et al., 2015, 2206.12267)
Multi-trait quantitative traits	Kronecker covariance $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$1 with grouped nonconvex penalties	(Liu et al., 2013)
Longitudinal outcomes	GRM random intercept plus subject-specific random intercepts/slopes	(St-Pierre et al., 19 Jan 2025)
Functional predictors	Penalized spline population slope and subject-specific random slope functions	(Liu et al., 2016)
Joint fixed/random selection in GLMMs	Group penalties on covariance-factor rows, MCECM estimation	(Heiling et al., 2023)
Low-rank random-effect structure	Latent factor decomposition $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$2	(Heiling et al., 2023)
Exact subset selection	Hierarchical $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$3 penalties on fixed effects and random-effect variances	(Thompson et al., 25 Jun 2025)

Binary-response extensions are especially important because they clarify a limitation of Gaussian PLMMs. The Sparse Probit Linear Mixed Model replaces the Gaussian likelihood by a correlated probit orthant probability and penalizes the sparse fixed effect with an $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$4 term (Mandt et al., 2015). The pGLMM framework argues that penalized LMMs are misspecified for binary traits because the working residual variance depends on the current mean, so a single Gaussian rotation is no longer valid; it replaces the PLMM likelihood by a penalized PQL-based GLMM while retaining a kinship-driven random effect (2206.12267). Taken together, these papers show that “PLMM for binary data” is usually shorthand for a generalized or probit extension rather than a literal Gaussian mixed model.

The multivariate direction is represented by the penalized-MTMM, which models stacked traits under

$$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$5

and then applies group MCP or sparse group MCP to trait-grouped SNP effects (Liu et al., 2013). The longitudinal direction is represented by penalized mixed models with multiple random effects, including a GRM-structured polygenic random intercept and subject-specific random intercepts and slopes for repeated measurements (St-Pierre et al., 19 Jan 2025). Functional extensions go further by treating the slope itself as a function $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$6 and adding subject-specific random slope functions $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$7, estimated by penalized splines and REML-EM updates (Liu et al., 2016).

A further recent use case is latent confounding and batch correction. "Penalized mixed models to adjust for batch effects and unobserved confounding in high dimensional regression" (Lu et al., 3 Oct 2025) decomposes hidden confounding into a bias term and an exogenous confounding noise term, and compares LASSO, PC-LASSO, and PLMM. Its main empirical conclusion is that PLMM is the most robust overall as confounding complexity increases. This suggests that covariance-based correction can be preferable to fixed-PC adjustment when the confounding structure is not effectively low-dimensional.

6. Software, scalability, and empirical regimes

A defining feature of recent PLMM work is the shift from methodological proposals to end-to-end software. The plmmr package is an R implementation specifically for genome-wide data with complex correlation structure (Peter et al., 3 Feb 2025). Its workflow includes reading PLINK or text input, constructing a standardized design with create_design(), fitting via plmm(), and optional tuning with cv_plmm(). The implementation builds on bigmemory, bigalgebra, and biglasso, stores coefficient paths in sparse Matrix format, and keeps the genotype matrix file-backed rather than loading it entirely into RAM.

The package paper reports two GWAS examples. In PennCath, with $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$8 participants and $$y=X\beta+u+\epsilon, \qquad u\sim N(0,\sigma_s^2K), \qquad \epsilon\sim N(0,\sigma_\epsilon^2I), \qquad u\perp \epsilon.$$9 autosomal SNPs after quality control, the full pipeline runtime on a single-core laptop ranged from about 1.5 minutes for the smallest subset to about 22 minutes for the full dataset, with preprocessing never exceeding about 5 minutes (Peter et al., 3 Feb 2025). In the much larger Pittsburgh Orofacial Cleft study, the dataset contained 10,545 participants, over 2,500 families, recruitment across 14 global sites on 5 continents, and 469,577 SNPs after QC; create_design() took about 9 minutes, eigendecomposition took 15.4 hours, model fitting took 16.2 hours, 5-fold CV selected $K=\frac{1}{p}XX^\top,$0, and the lasso selected 53 SNPs, including hits in NTN1, PAX7, IRF6, and FOXE1 (Peter et al., 3 Feb 2025).

Several software ecosystems now cover neighboring parts of the design space. glmmPen supports Gaussian, Binomial, and Poisson families and performs simultaneous fixed- and random-effect selection by penalizing fixed coefficients and grouped covariance-factor rows; its Gaussian family is directly relevant as a penalized Gaussian mixed model (Heiling et al., 2023). The latent-factor extension is implemented in the same package and is intended to make high-dimensional random-effect selection more scalable (Heiling et al., 2023). PenalizedGLMM.jl brings the PLMM logic to binary GWAS through a penalized GLMM that can scale to approximately 300,000 SNPs by null-model variance estimation, PQL, and a single spectral decomposition trick (2206.12267). glmmsel provides hierarchical $K=\frac{1}{p}XX^\top,$1-regularized subset selection for linear mixed models, with paths of 100 models and around $K=\frac{1}{p}XX^\top,$2 predictors reported in seconds to minutes (Thompson et al., 25 Jun 2025).

The empirical use regimes differ. plmmr is targeted at joint sparse modeling in correlated high-dimensional Gaussian settings, especially GWAS (Peter et al., 3 Feb 2025). glmmPen and PenalizedGLMM.jl are aimed at generalized outcomes (Heiling et al., 2023, 2206.12267). glmmsel is aimed at exact subset selection in grouped high-dimensional data (Thompson et al., 25 Jun 2025). This suggests that “PLMM software” is now better understood as a toolkit family indexed by outcome type, correlation structure, and the desired form of sparsity.

7. Limitations, controversies, and open problems

The literature is explicit that covariance selection in mixed models is harder than variable selection in ordinary regression. Positive semidefinite constraints and boundary points at zero variance components make covariance-model selection numerically and theoretically delicate (Müller et al., 2013). Cholesky-based parameterizations help preserve legality, but the review literature also notes their order dependence: different orderings of random effects can induce different selected models in finite samples (Müller et al., 2013).

Computational bottlenecks remain important. In preconditioning-based GWAS PLMMs, the eigendecomposition of $K=\frac{1}{p}XX^\top,$3 is the dominant cost as $K=\frac{1}{p}XX^\top,$4 grows, and file-backing does not remove the expense of storing and diagonalizing an $K=\frac{1}{p}XX^\top,$5 relationship matrix (Peter et al., 3 Feb 2025). In longitudinal penalized mixed models, replacing the full GRM by a sparse GRM yields major computational gains but biases variance-component estimates, typically shifting variance from weak genetic links into non-genetic random intercept terms (St-Pierre et al., 19 Jan 2025). Exact subset-selection methods scale much better than older joint-selection mixed-model procedures, but they remain nonconvex and rely on local search rather than global optimality certificates (Thompson et al., 25 Jun 2025).

Theoretical coverage is uneven. Some papers establish model-selection consistency or oracle properties for fixed-effect or fixed/random selection under explicit regularity conditions (Fan et al., 2012, Ghosh et al., 2016). Others are primarily optimization papers and do not address selection consistency, de-biasing, or asymptotic distributions of penalized estimators (Fort et al., 2017). The new $K=\frac{1}{p}XX^\top,$6 subset-selection work provides a finite-sample KL-divergence bound, but not support-recovery consistency (Thompson et al., 25 Jun 2025). This suggests that “theory for PLMMs” is fragmented across distinct inferential goals.

A separate controversy concerns outcome type. The generalized-model literature argues that applying Gaussian penalized LMM machinery to binary outcomes is fundamentally misspecified because the working covariance depends on fitted means; binary GWAS therefore motivate penalized GLMM or probit extensions rather than direct Gaussian PLMM use (2206.12267, Mandt et al., 2015). In practice, some software still allows binary outcomes to be treated numerically in a linear model, but the package papers identify that as a current limitation rather than a settled solution (Peter et al., 3 Feb 2025).

Finally, cross-validation itself is now recognized as a first-order methodological issue. Full CV is computationally more expensive because every fold must recompute covariance estimation, eigendecomposition, rotation, and prediction, yet the recent literature argues that this cost is the price of valid tuning (Peter et al., 18 Mar 2025). A plausible implication is that future PLMM development will continue to balance statistical validity against increasingly aggressive computational shortcuts.