---
title: Penalized Linear Mixed Models
url: https://www.emergentmind.com/topics/penalized-linear-mixed-models-plmms
type: topic
---

# Penalized Linear Mixed Models

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 [2502.01577]. 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 [1211.0457, 1306.2427].

## 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 [2502.01577]. 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 [2502.01577]. 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
\[
\tilde y\sim N(\tilde X\beta,I).
\]
On the transformed scale, the errors are independent and standard penalized regression machinery becomes applicable [2502.01577]. 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 [2502.01577]. 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
\[
\hat\beta
=
\arg\min_\beta
\left\{
\frac{1}{2}\|\tilde y-\tilde X\beta\|_2^2
+
\sum_{j=1}^p P_\lambda(|\beta_j|)
\right\},
\]
with the lasso as the special case
\[
\hat\beta
=
\arg\min_\beta
\left\{
\frac{1}{2}\|\tilde y-\tilde X\beta\|_2^2
+
\lambda\|\beta\|_1
\right\}.
\]
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 [2502.01577, 1607.02883].

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" [1607.02883], the objective is the penalized negative marginal log-likelihood
\[
Q_{n,\lambda}(\beta,\eta)
=
-l_n(\beta,\eta)+\sum_{j=1}^p P_{n,\lambda}(|\beta_j|),
\]
where \(\eta=(\sigma,\theta)\) 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 \(\beta\).

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 [1211.0457]. 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" [1306.2427] 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 [2305.08204]. "Scalable Subset Selection in Linear Mixed Models" [2506.20425] goes further by replacing continuous shrinkage with a hierarchical \(\ell_0\) criterion,
\[
l(\beta,\gamma)+\lambda\alpha\|\beta\|_0+\lambda(1-\alpha)\|\gamma\|_0,
\]
subject to
\[
\beta_k=0\Rightarrow \gamma_k=0,
\]
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 [1306.2427]. 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
\[
P_z=(I+\sigma^{-2}Z\mathcal G Z^\top)^{-1}
\]
is replaced by a proxy-based matrix
\[
\widetilde P_z=(I+ZMZ^\top)^{-1},
\]
and they prove model-selection consistency for both fixed and random effects under conditions on the proxy matrix [1211.0457]. 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 [1607.02883]. 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 \(s\log p=o(n)\), 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 \(K\), a **fit** step that uses coordinate descent on the transformed problem, and a **format** step that returns coefficients on the original scale [2502.01577]. If \(K=UDU^\top\), then
\[
\Sigma=U(\sigma_s^2D+\sigma_\epsilon^2I)U^\top
\quad\text{and}\quad
\Sigma^{-1/2}=U(\sigma_s^2D+\sigma_\epsilon^2I)^{-1/2}U^\top,
\]
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" [1704.08891] 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
\[
S_{n+1}^{sa}=(1-\delta_{n+1})S_n^{sa}+\frac{\delta_{n+1}}{m_{n+1}}\sum_{j=1}^{m_{n+1}}S(Z_{j,n}),
\]
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 [2305.08204]. The latent-factor extension [2305.08201] replaces a high-dimensional random-effects covariance by a low-rank factorization \(\Sigma=BB^\top\), reducing the latent dimension from \(q\) to \(r\ll q\). 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 \(\ell_0\) mixed-model procedure of [2506.20425] applies block coordinate descent to pairs \((\beta_k,\gamma_k)\), 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 \(O(n_i^3)\) to \(O(n_i^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 \((X_1,y_1)\), test design \(X_2\), and block covariance
\[
\hat\Sigma=
\begin{bmatrix}
\hat\Sigma_{11} & \hat\Sigma_{12}\\
\hat\Sigma_{21} & \hat\Sigma_{22}
\end{bmatrix},
\]
the recommended validation prediction is
\[
\hat y_2=X_2\hat\beta+\hat\Sigma_{21}\hat\Sigma_{11}^{-1}(y_1-X_1\hat\beta),
\]
which is the BLUP-adjusted predictor emphasized in the recent cross-validation literature [2503.14374]. In GWAS-oriented software, BLUP-based prediction is also exposed directly through `predict()` and cross-validation wrappers [2502.01577].

The cross-validation problem is now recognized as methodologically nontrivial. "Cross-Validation in Penalized Linear Mixed Models: Addressing Common Implementation Pitfalls" [2503.14374] 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 \(0.92\), FDR \(0.63\), NVAR \(16\), and RSEE \(0.81\); **inner CV** had TDR \(0.92\), FDR \(0.80\), NVAR \(29\), and RSEE \(0.84\); **outer CV** selected \(648\) variables with FDR \(0.99\) and RSEE \(4.30\) [2503.14374]. In the PennCath real-data example, the same paper reports that outer CV selected \(556\) variables at \(\lambda_{\min}\), whereas full CV selected \(5\). 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 [2503.14374].

## 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 | [1507.04777], [2206.12267] |
| Multi-trait quantitative traits | Kronecker covariance \(K\otimes\Sigma_g + I_n\otimes\Sigma_e\) with grouped nonconvex penalties | [1305.4413] |
| Longitudinal outcomes | GRM random intercept plus subject-specific random intercepts/slopes | [2501.11083] |
| Functional predictors | Penalized spline population slope and subject-specific random slope functions | [1601.01039] |
| Joint fixed/random selection in GLMMs | Group penalties on covariance-factor rows, MCECM estimation | [2305.08204] |
| Low-rank random-effect structure | Latent factor decomposition \( \Sigma=BB^\top \) | [2305.08201] |
| Exact subset selection | Hierarchical \(\ell_0\) penalties on fixed effects and random-effect variances | [2506.20425] |

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 \(\ell_1\) term [1507.04777]. 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
\[
H=K\otimes\Sigma_g+I_n\otimes\Sigma_e
\]
and then applies group MCP or sparse group MCP to trait-grouped SNP effects [1305.4413]. 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 [2501.11083]. Functional extensions go further by treating the slope itself as a function \(\beta(t)\) and adding subject-specific random slope functions \(b_i(t)\), estimated by penalized splines and REML-EM updates [1601.01039].

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" [2510.03531] 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 [2502.01577]. 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 \(1{,}401\) participants and \(696{,}644\) 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** [2502.01577]. 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 \(\lambda\), and the lasso selected **53 SNPs**, including hits in **NTN1**, **PAX7**, **IRF6**, and **FOXE1** [2502.01577].

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 [2305.08204]. The latent-factor extension is implemented in the same package and is intended to make high-dimensional random-effect selection more scalable [2305.08201]. **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 \(\ell_0\)-regularized subset selection for linear mixed models, with paths of 100 models and around \(10^4\) predictors reported in seconds to minutes [2506.20425].

The empirical use regimes differ. **plmmr** is targeted at joint sparse modeling in correlated high-dimensional Gaussian settings, especially GWAS [2502.01577]. **glmmPen** and **PenalizedGLMM.jl** are aimed at generalized outcomes [2305.08204, 2206.12267]. **glmmsel** is aimed at exact subset selection in grouped high-dimensional data [2506.20425]. 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 [1306.2427]. 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 [1306.2427].

Computational bottlenecks remain important. In preconditioning-based GWAS PLMMs, the eigendecomposition of \(K\) is the dominant cost as \(n\) grows, and file-backing does not remove the expense of storing and diagonalizing an \(n\times n\) relationship matrix [2502.01577]. 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 [2501.11083]. 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 [2506.20425].

Theoretical coverage is uneven. Some papers establish model-selection consistency or oracle properties for fixed-effect or fixed/random selection under explicit regularity conditions [1211.0457, 1607.02883]. Others are primarily optimization papers and do not address selection consistency, de-biasing, or asymptotic distributions of penalized estimators [1704.08891]. The new \(\ell_0\) subset-selection work provides a finite-sample KL-divergence bound, but not support-recovery consistency [2506.20425]. 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, 1507.04777]. 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 [2502.01577].

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 [2503.14374]. A plausible implication is that future PLMM development will continue to balance statistical validity against increasingly aggressive computational shortcuts.

Source: https://www.emergentmind.com/topics/penalized-linear-mixed-models-plmms