stLMM: Extensions of Linear Mixed Models

• stLMM is a family of generalized linear mixed model extensions that relax standard Gaussian and continuity assumptions to address various data challenges.
• It includes methods like WarpedLMM for phenotype transformation, Probit-LMM for binary outcomes with sparsity, and Bayesian spatial models for ecological estimation.
• Each variant targets specific applications, offering improved heritability estimation, classification accuracy, and robustness in clustered, longitudinal, or spatial data.

stLMM is not a single universally standardized term. In the arXiv literature, it is used for several distinct linear mixed model extensions: a transformed-phenotype LMM in which the phenotype transformation is estimated from the data; a sparse LMM generalized to binary outcomes through a probit link; an R package for Bayesian spatial and space-time linear mixed models in ecological small-area estimation; and ST-LMM denoting skew-$t$ linear mixed models for clustered or longitudinal data (Fusi et al., 2014, Mandt et al., 2015, Finley, 3 Jul 2026, Schumacher et al., 2021). The common thread is the retention of mixed-model structure—fixed effects, random effects, and structured covariance—while relaxing a different modeling bottleneck: Gaussianity on the observed scale, continuity of the response, absence of spatial-temporal dependence, or symmetry and light tails.

1. Terminological scope and principal usages

The label “stLMM” is therefore context-dependent. In statistical genetics, it may denote a transformed-phenotype or sparse LMM workflow; in ecological statistics, it is the formal name of a software package; and in robust mixed-effects modeling, “ST-LMM” refers to a skew-$t$ distributional extension. This suggests that the term functions more as a family resemblance across LMM generalizations than as a single canonical model class.

Usage of “stLMM”	Core idea	Representative source
Transformed-phenotype LMM	Learn a monotone invertible transformation jointly with the LMM	WarpedLMM (Fusi et al., 2014)
Sparse binary-trait LMM	Probit mixed model with $\ell_1$-penalized fixed effects	Probit-LMM (Mandt et al., 2015)
Bayesian spatial and space-time LMM	Formula-based package with GP, NNGP, CAR, DAGAR, AR(1), and space-time effects	stLMM package (Finley, 3 Jul 2026)
Skew-$t$ linear mixed model	CFUST-based random effects and heavy-tailed errors	ST-LMM (Schumacher et al., 2021)

A useful organizing distinction is the object being generalized. WarpedLMM generalizes the phenotype scale; Probit-LMM generalizes the response likelihood; the stLMM package generalizes the latent effect structure for spatial, temporal, and space-time estimation; and the CFUST-based ST-LMM generalizes the joint distribution of random effects and errors.

2. Transformed-phenotype stLMM: WarpedLMM

In the transformed-phenotype sense, stLMM refers to a linear mixed model in which the observed phenotype $y$ is mapped to a latent Gaussian phenotype $z$ through an unknown monotone invertible transformation estimated jointly with the mixed-model parameters (Fusi et al., 2014). On the latent scale, the model is the standard LMM

$$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$

with

$$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$

where the genomic relationship matrix is

$K = \frac{1}{S}GG^\top.$

The observed phenotype is related to this latent variable by

$z_n = f(y_n;\psi),$

with $t$0 unknown, monotone, and invertible.

WarpedLMM uses the parametric warping family introduced following Snelson et al.: $t$1 under the constraints $t$2, $t$3, and $t$4. In all experiments, $t$5. Because $t$6 is increasing and all coefficients are nonnegative, the transformation is monotone increasing and hence invertible; $t$7 is obtained numerically by Newton–Raphson. The likelihood on the observed scale includes the Jacobian term,

$t$8

with $t$9. The Jacobian is essential because it prevents degenerate transformations and gives the correct likelihood on the original scale. Estimation is by maximizing a penalized likelihood over $\ell_1$0, combining standard (R)ML estimation of the variance components with numerical optimization over the transformation parameters.

This formulation differs sharply from ad hoc transformed-LMM practice. A standard transformed LMM chooses a fixed transformation $\ell_1$1—log, Box–Cox, or rank-based inverse normal—before fitting the mixed model. WarpedLMM instead estimates $\ell_1$2 in the presence of covariates and genetic covariance, thereby targeting Gaussian residuals after fitting the full model. The paper explicitly compares WarpedLMM with Box–CoxLMM, rank-based transformation followed by LMM, and the semi-parametric rank-based pipeline of Zhou and Stephens.

Once the transformation is estimated, association testing proceeds on the transformed phenotype $\ell_1$3 using standard LMM GWAS machinery,

$\ell_1$4

with Wald or likelihood-ratio testing for $\ell_1$5. Heritability is defined on the transformed scale as

$\ell_1$6

Prediction is performed by BLUP on the latent scale and then mapped back via the inverse warp: $\ell_1$7 A major practical consequence is that, unlike rank-based preprocessing, the method preserves a principled inverse map to the original phenotype units.

The empirical results are reported across simulations and real data from mouse, yeast, and human cohorts. In 50,000 simulations using HapMap3 chromosome 22 genotypes, WarpedLMM yields almost unbiased heritability estimates across varied $\ell_1$8, sample size, number of causal variants, covariate variance explained, and degree of nonlinearity, whereas the standard LMM systematically underestimates heritability and the bias worsens at higher true $\ell_1$9 and with more causal variants. In the Valdar et al. mouse data, 18 of 47 traits show significant heritability differences between WarpedLMM and the untransformed LMM, and in 17 of those 18 the WarpedLMM estimate is higher; 10-fold cross-validation shows consistently higher out-of-sample $t$0. In the Bloom et al. yeast data, 17 of 45 phenotypes show significant $t$1 differences, mostly increases, with corresponding improvements in prediction. In the Northern Finnish birth cohort analyzed by Sabatti et al., all methods retain genomic control $t$2, but WarpedLMM finds 6 distinct loci for LDL whereas the comparison methods find 3; it also finds 3 QTLs for HDL, with one missed by the alternatives. The paper further reports p-value correlation $t$3 with the semi-parametric rank-based transformation approach of Zhou and Stephens, while emphasizing the invertibility and joint likelihood-based estimation of WarpedLMM.

3. Sparse and binary-trait stLMM: Probit-LMM

A second usage of stLMM appears in sparse or regularized LMMs for feature selection, especially when extended to binary phenotypes. In this setting, the relevant model is the Sparse Probit Linear Mixed Model, or Probit-LMM, which generalizes the LMM paradigm to binary responses while preserving kinship-type covariance correction for confounding (Mandt et al., 2015). The generative model is

$t$4

with

$t$5

and in the main linear-kernel form

$t$6

For some experiments, an additional side-information kernel is included: $t$7

The sparse component enters through an $t$8 penalty on the fixed-effect coefficient vector $t$9. After absorbing the label signs into $y$0 and $y$1, the objective becomes

$y$2

The paper proves that this objective is convex. It also gives an equivalent two-weight-vector representation in which a dense Gaussian component $y$3 captures background population structure and the sparse vector $y$4 captures interpretable feature-specific effects: $y$5 This mirrors the sparse-plus-polygenic decomposition commonly associated with LMM-Lasso, but with a probit likelihood rather than Gaussian regression. The authors explicitly position Probit-LMM as a natural extension of LMM-Lasso to binary outcomes and note that naively applying LMM-Lasso to binary labels leads to lower predictive accuracies.

The technical difficulty is that the correlated probit likelihood has no closed form. Two scalable approximate inference schemes are proposed. Both use ADMM to handle the smooth-plus-$y$6 structure,

$y$7

with alternating updates in $y$8, $y$9, and the dual variable $z$0. The faster approximation is a MAP treatment of the dense random component $z$1. The main method is an approximate EM algorithm in which the E-step requires moments of a truncated multivariate normal,

$z$2

approximated by Expectation Propagation. The resulting gradient and Hessian are expressed through the posterior moments $z$3 and $z$4, after which Newton steps are embedded in the ADMM outer loop. The paper typically performs only one Newton step per ADMM iteration.

Confounder correction remains conceptually identical to standard LMM practice: confounding variation is represented through the covariance $z$5. In the tuberculosis experiment, age is incorporated through an RBF kernel $z$6; in other settings, the linear kernel $z$7 captures similarity or kinship structure. The random-effect term is therefore handled as correlated Gaussian noise, while sparsity is reserved for the fixed effects that are meant to have a more specific or causal interpretation. The authors emphasize that, under the MAP approximation, $z$8 and $z$9 enter symmetrically in the likelihood and are distinguished only by their $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$0 and $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$1 regularizers.

Empirically, the method is evaluated on synthetic data and three real datasets. On tuberculosis gene-expression data with 40 TB cases, 103 controls, 48,803 features, and an age-based side kernel, Probit-LMM improves AUC by up to 12 percentage points over sparse Probit, by up to 3 points over GP classification, by up to 7 points over LMM-Lasso, and by up to 7 points over the MAP approximation. In malware detection with 200 Android apps and 545,333 binary features, the reported normalized low-FPR metric $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$2 is $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$3 for Probit-LMM, compared with $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$4 for the MAP approximation, $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$5 for sparse Probit, $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$6 for GP classification, and $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$7 for LMM-Lasso. In Arabidopsis flowering time, binarized from 199 accessions and 216,130 SNPs, Probit-LMM attains $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$8 AUC, modestly above the other binary methods and clearly above LMM-Lasso at $$\mathbf{z} = X\beta + u + \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$9. The paper also reports that selected features under Probit-LMM are less correlated with the first principal component of the kernel than those from sparse Probit, and that bootstrap feature sets are more stable.

4. stLMM as a Bayesian spatial and space-time modeling package

In ecological statistics, stLMM is the explicit name of an R package for Bayesian linear mixed models with spatial, temporal, and space-time latent effects, designed for small-area ecological estimation (Finley, 3 Jul 2026). The motivating use case is programs such as the US Forest Inventory and Analysis, where analysts require estimates for counties, years, ecoregions, or other reporting domains that may have few or no observations in a given period. The package is intended to support both area-level and unit-level small-area estimation within one framework.

For Gaussian responses, the package uses the latent model

$$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$0

where $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$1 denotes explicit i.i.d. grouped random effects, $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$2 stacks structured latent processes, and $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$3 has either homoskedastic Gaussian variance or a diagonal precision structure

$$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$4

This residual formulation is central for direct-estimate models: $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$5 may be fixed as $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$6 for known sampling variances, or modeled through resid() when the variances themselves are uncertain. The structured latent processes are combined as

$$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$7

and the key computational matrix is

$$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$8

The package provides a shared formula interface across a substantial range of latent-effect classes: iid() for grouped random effects, gp() for dense point-referenced Gaussian processes, nngp() for nearest-neighbor Gaussian processes, ar1() for temporal autoregression, car() for areal CAR models, dagar() for DAGAR models, and the separable space-time terms car_time() and dagar_time(). Structured varying coefficients are specified through interactions such as x:car(...) or x:car_time(...), with the mapping matrix $$\mathbf{z} \sim \mathcal{N}\bigl(X\beta,\; \sigma_g^2 K + \sigma_e^2 I\bigr),$$9 multiplying latent effects by covariates. For non-Gaussian responses, binomial and fixed-size negative binomial likelihoods are handled through Pólya–Gamma data augmentation, yielding a conditionally Gaussian form with a diagonal working precision updated within MCMC.

The computational design is based on sparse precision matrices and collapsed sampling. Most structured processes—CAR, DAGAR, CAR-time, DAGAR-time, AR(1), and NNGP—supply sparse precision blocks. During fitting, the structured latent processes $K = \frac{1}{S}GG^\top.$0 are integrated out, reducing the MCMC state to fixed effects, iid random effects, residual-variance parameters, and latent-process hyperparameters. After fitting, recover() samples $K = \frac{1}{S}GG^\top.$1 from its posterior conditional distribution, and the recovered draws are used for fitted(), diagnostics, prediction, aggregation, and mapping. The implementation uses CHOLMOD for sparse Cholesky factorization and BLAS/LAPACK for dense linear algebra. The package characterizes this as a “one posterior-draw workflow,” because the same posterior draws are propagated through prediction and subsequent aggregation.

A central modeling feature is the retention of rows with missing responses. Such rows can remain in the data with their covariates and support indices intact; they contribute to $K = \frac{1}{S}GG^\top.$2 and $K = \frac{1}{S}GG^\top.$3 but not to the likelihood through $K = \frac{1}{S}GG^\top.$4, and later become prediction targets. This is particularly important for small-area estimation because entire county-year combinations with no direct estimate can still receive model-based posterior predictions with uncertainty.

The main application in the paper is a Washington county-year biomass example using FIA data and tree canopy cover. The fitted model is

$K = \frac{1}{S}GG^\top.$5

where $K = \frac{1}{S}GG^\top.$6 is a county-level CAR spatially varying coefficient on tree canopy cover, $K = \frac{1}{S}GG^\top.$7 is a separable county-year CAR-time process, and the residual variance is modeled through the scaled-variance specification

$K = \frac{1}{S}GG^\top.$8

The corresponding stLMM() call combines county_mean_tcc_scaled:car(...), car_time(county_fips, year, ...), and resid(model = "scaled", ...). Posterior draws of

$K = \frac{1}{S}GG^\top.$9

are then recovered for all county-year combinations, including those without direct estimates. The package article series states that reproducible source code, data, diagnostics, and related model variants are provided.

5. ST-LMM as a canonical fundamental skew-$z_n = f(y_n;\psi),$0 linear mixed model

A fourth usage is ST-LMM in the sense of the canonical fundamental skew-$z_n = f(y_n;\psi),$1 linear mixed model developed by Matos, Schumacher, and Cabral for longitudinal and clustered data (Schumacher et al., 2021). The starting point is the standard mixed model

$z_n = f(y_n;\psi),$2

but the Gaussian assumptions on random effects and errors are replaced by a canonical fundamental skew-$z_n = f(y_n;\psi),$3 distribution. The motivation is robustness to skewness and heavy tails, particularly when Gaussian LMMs are sensitive to outliers or asymmetric subject-specific effects.

The construction begins with the CFUSN and CFUST families. If $z_n = f(y_n;\psi),$4, $z_n = f(y_n;\psi),$5, and $z_n = f(y_n;\psi),$6 is a $z_n = f(y_n;\psi),$7 skewness matrix, then

$z_n = f(y_n;\psi),$8

defines a CFUSN random vector. Introducing a Gamma scale mixture,

$z_n = f(y_n;\psi),$9

yields the CFUST family. In the mixed-model application, the joint vector of random effects and within-subject errors is assumed to satisfy

$t$00

where $t$01 is the random-effects covariance matrix, $t$02 is the within-subject covariance, $t$03 is the skewness matrix, $t$04 is the skew dimension, and $t$05 is the degrees of freedom parameter. The scalar $t$06 is chosen so that $t$07.

Under this construction, $t$08 is skew-$t$09, $t$10 is symmetric Student-$t$11, and $t$12 has a CFUST marginal density. The hierarchical representation is central: $t$13

$t$14

$t$15

This closed-form hierarchy enables an ECME algorithm for maximum likelihood estimation. The E-step requires conditional expectations such as $t$16, $t$17, and truncated-$t$18 moments for the latent skewness variables $t$19. The M-step then updates $t$20, $t$21, $t$22, and $t$23, while $t$24 is updated by directly maximizing the observed log-likelihood. The paper also discusses standard-error estimation via Louis’ method and empirical Bayes prediction of random effects.

An important property of this ST-LMM is nesting. By parameter restriction it contains the Gaussian LMM, the $t$25-LMM, skew-normal LMMs, Sahu–Dey–Branco skew-$t$26 cases, and the Azzalini–Capitanio skew-$t$27 case. This permits likelihood-based model comparison. In the schizophrenia application, six variants are fitted: SN-LMM with $t$28, SN-LMM with SDB skewness, SN-LMM with $t$29, ST-LMM with $t$30, ST-LMM with SDB skewness, and ST-LMM with $t$31. The reported AIC values are 1579.12, 1605.26, 1581.44, 1499.84, 1532.56, and 1499.02, respectively, with the full ST-LMM $t$32 variant giving the smallest AIC. In that application, $t$33, $t$34, and $t$35 are significant, whereas $t$36 and $t$37 are not significant at 5%, supporting the stated equivalence conclusion for the therapy comparison.

The paper’s simulation studies further indicate that the ECME estimator behaves well in finite samples, with bias shrinking as $t$38 increases and Louis-based standard errors broadly tracking Monte Carlo variability. A practical point emphasized by the authors is initialization: different starting strategies can materially affect convergence, and a “best of” multi-start strategy is recommended.

6. Cross-cutting distinctions and interpretive cautions

The principal source of confusion around stLMM is nomenclatural rather than mathematical. The four cited usages address different statistical problems. WarpedLMM is a transformed-phenotype model for non-Gaussian continuous traits; Probit-LMM is a sparse binary-response model for confounded classification; the stLMM package is a Bayesian spatial and space-time framework for ecological small-area estimation; and ST-LMM is a robust skew-$t$39 mixed model for clustered or longitudinal responses (Fusi et al., 2014, Mandt et al., 2015, Finley, 3 Jul 2026, Schumacher et al., 2021).

Their latent structures are correspondingly different. WarpedLMM uses a standard genetic LMM after a learned monotone warp and relies on REML-style variance-component estimation plus optimization over transformation parameters. Probit-LMM retains a kinship-style covariance but replaces Gaussian regression with a latent-Gaussian probit likelihood and combines ADMM, Newton updates, and EP. The ecological stLMM package builds a block-diagonal precision representation for structured spatial, temporal, and space-time effects, collapses those effects during MCMC, and recovers them afterward for prediction and aggregation. The CFUST-based ST-LMM changes the joint distribution of $t$40 and $t$41, introducing latent half-normal and Gamma variables and estimating parameters via ECME.

A second distinction concerns the inferential target. WarpedLMM emphasizes GWAS power, heritability estimation, and phenotype prediction on the original scale. Probit-LMM emphasizes classification accuracy, confounder-robust feature selection, and sparse interpretability. The stLMM package emphasizes posterior-draw propagation, domain prediction with missing-response rows retained as targets, and direct-estimate versus unit-level modeling within one interface. ST-LMM emphasizes robustness to skewness and heavy tails, likelihood-based model selection among nested LMM families, and improved fit in longitudinal settings with asymmetric subject heterogeneity.

Accordingly, “stLMM” should be read only after the surrounding context is fixed. In statistical genetics it may refer either to a transformed LMM or to a sparse mixed model extension; in ecological Bayes modeling it is a package name; and in robust mixed-effects theory ST-LMM denotes a skew-$t$42 distributional specification. The term therefore identifies a broad pattern of LMM generalization rather than a unique model, software system, or inferential protocol.