stLMM: Bayesian Spatial and Space-Time Linear Mixed Models for Small-Area Ecological Estimation

Abstract: stLMM is an R package for Bayesian linear mixed models with spatial, temporal, and space-time latent effects. It provides a common formula interface for independent and identically distributed (iid) grouped effects, autoregressive (AR) temporal effects, Gaussian process (GP) and nearest-neighbor Gaussian process (NNGP) point-referenced effects, conditional autoregressive (CAR) and directed acyclic graph autoregressive (DAGAR) areal effects, separable areal space-time effects, and structured varying coefficients. The package is designed for ecological small-area estimation workflows in which analysts must move between direct-estimate and unit-level models, combine sampling variances or residual-variance models with spatial and temporal borrowing, and retain missing response rows as prediction targets. A shared sparse-precision implementation underlies the model terms. Structured latent effects are collapsed during fitting, then recovered or retained for fitted values, diagnostics, prediction, and posterior summaries. This gives users one posterior-draw workflow for propagating uncertainty from model fitting through prediction and aggregation. The package is demonstrated with a Washington county biomass example from the package article series, using Forest Inventory and Analysis (FIA) data and tree canopy cover to estimate county-year biomass means. The full article series provides reproducible source code, data, diagnostics, and related model variants.

• The paper proposes stLMM as a unified Bayesian framework integrating spatial, temporal, and space-time random effects for ecological small-area estimation.
• It employs a collapsed MCMC sampler and sparse precision matrix architecture to efficiently manage Gaussian, binomial, and count data models.
• The approach enables robust uncertainty propagation and performance benchmarking, demonstrated via county-year biomass estimation in Washington.

Bayesian Spatial and Space-Time Linear Mixed Models for Small-Area Ecological Estimation: The stLMM R Package

Overview and Motivation

The stLMM R package provides a unified Bayesian framework for spatial, temporal, and space-time ecological small-area estimation (SAE) using linear mixed models with structured latent effects. SAE is a persistent challenge in ecological monitoring where reporting domains (e.g., counties, years) often have insufficient direct sample support. Classical area-level (Fay-Herriot) and unit-level (Battese-Harter-Fuller) models are foundational but lack the flexibility necessary to combine direct and unit-level modeling, allow complex spatial or temporal random effects, and propagate uncertainty through aggregation on arbitrary domains. stLMM addresses this gap by supporting a diverse array of linear mixed model terms, unified under a sparse precision matrix architecture, and by providing a robust interface for model specification, fitting, and prediction.

Modeling Framework and Implementation

stLMM enables modeling of Gaussian, binomial, and count data via a general linear mixed model expressed as:

$y = X\beta + Zb + A w + \varepsilon,$

where $X\beta$ are fixed effects; $Zb$ are iid random effects; $w$ is a stacked vector of structured latent processes (e.g., spatial, temporal, or space-time random effects); $A$ is the process mapping matrix; and $\varepsilon$ encodes residual variance, supporting both homoskedasticity and complex, row-specific structures (e.g., direct-estimate sampling variance). Posterior inference leverages a collapsed MCMC scheme, integrating out structured latent processes during sampling, leading to efficiency gains in memory and computation by reducing the primary MCMC state space.

Structured latent effects are modular, with terms including:

• iid(): for grouped random effects,
• gp()/nngp(): for Gaussian (dense) and nearest-neighbor (sparse, scalable) point-referenced spatial or space-time processes,
• car()/dagar(): for areal spatial dependence (CAR, DAGAR, including space-time via car_time/dagar_time),
• ar1(): for temporal AR(1) effects,
• Structured varying coefficients: via interactions with process terms.

All process terms, except dense gp(), are assembled as sparse precision matrices, with computation built atop R's Matrix/CHOLMOD APIs, leveraging BLAS/LAPACK for dense sub-tasks. Data augmentation (e.g., Pólya-Gamma for logistic/negative binomial models) unifies the algorithmic paradigm for non-Gaussian outcomes, maintaining conditional normality for latent process inference.

User Interface and Workflow

User interaction centers on R's formula syntax, extended to encode both standard and structured latent terms. The package provides a strictly minimal, yet expressive, workflow:

fit <- stLMM(formula, data = dat, ...) rec <- recover(fit, ...) pred <- predict(rec, newdata = newdat, ...)

Model specification supports both area-level and unit-level responses, explicit modeled or known sampling variances (direct-estimate), structured process recovery for posterior summaries, and full uncertainty propagation through posterior sampling. Posterior inferences are returned as draws, facilitating downstream aggregation and transformation.

Comparative Software Perspective

stLMM is specialized relative to general probabilistic programming frameworks (e.g., Stan, NIMBLE), INLA, and domain-specific spatial packages (sdmTMB, FRK, spBayes, spNNGP, CARBayes, geostan, hSDM, Hmsc, spOccupancy, spAbundance). Its unique focus is on seamless movement between model classes (point-referenced, areal, space-time), integrated handling of direct and unit-level models, and a robust pipeline for spatially/temporally indexed prediction and summary—all with a unified posterior prediction machinery suitable for ecological small-area estimation.

While it is less general in terms of model expressivity, it fills an unaddressed methodological niche: the ability to compare and combine multiple classical and modern SAE modeling strategies in a shared Bayesian interface, propagating posterior uncertainty end-to-end, and efficiently utilizing sparse precision computations.

Application: County-Year Biomass Estimation in Washington State

The practical utility of stLMM is demonstrated via a Washington county-year biomass SAE analysis using US Forest Service FIA data and county tree canopy cover (TCC) as an auxiliary covariate. Here, the response is county-year direct estimates of biomass and associated sampling variances.

A representative fitted model is:

$$\widehat{y}_i = \beta_0 + \{\beta_1 + u(a_i)\} x_{a_i, t_i} + w(a_i, t_i) + e_i, \quad e_i \sim N(0, \tau^2_i),$$

where the TCC coefficient is allowed to vary spatially via a county-level CAR process on $u(a)$, and a separable CAR-time process $w(a, t)$ induces spatial-temporal smoothing. Residual variance modeling incorporates both the estimated sampling variance and a calibrated shrinkage parameter, enabling both sharp and diffuse uncertainty quantification for ill-supported domains.

Direct estimates show high variance and data sparsity across both spatial and temporal domains, while the posterior model means (smoothed via space-time CAR process and TCC covariate) provide robust, spatially and temporally coherent inferences:

Figure 1: Washington county-year direct estimates and posterior means from the scaled-variance CAR-time model with a spatially varying TCC coefficient. Grey counties have no model-ready direct estimate in that year.

The spatially varying TCC coefficient surface summarizes the spatial heterogeneity in the canopy cover–biomass association:

Figure 2: Posterior mean tree canopy cover (TCC) coefficient, in Mg/ha per one standard deviation increase in county mean TCC.

Posterior summaries for selected counties across time illustrate improved uncertainty quantification; direct estimates are plotted with point sizes scaled by the available FIA plot count:

Figure 3: Posterior county-year biomass means and 95% credible intervals for four counties, with direct estimates shown as points sized by FIA plot count.

This case study highlights stLMM’s ability to (1) accommodate spatially varying coefficient models, (2) integrate sampling variance modeling with scalable spatial and space-time smoothing, and (3) support missing response imputation via posterior predictive draws. The end-to-end posterior propagation for both domain-level estimates and uncertainty is retained.

Implications, Limitations, and Future Directions

stLMM provides a practical and highly extensible toolkit for ecological SAE, enabling flexible model-based borrowing across spatial, temporal, and biophysical domains. The approach is especially pertinent for national forest inventory applications, species distributions, and carbon stock monitoring, where direct estimation in small areas is underpowered and out-of-domain prediction with uncertainty is essential.

The main trade-off is computational. The collapsed sampler and sparse matrix operations are efficient in moderate-size analyses but may not always outperform highly specialized model-specific samplers in very large-scale settings or when a single model class suffices. Nonetheless, the balance of transparency, extensibility, and unified workflow—especially for analysts performing sensitivity analyses or cross-model comparisons—is a strong asset.

Potential future developments include broader support for non-Gaussian likelihoods, further scalability improvements for massive space-time datasets (e.g., improved NNGP and spatiotemporal process implementations), and extensions to high-dimensional, multi-response scenarios (such as stacking or joint species distribution models). Algorithmic advances in sparse linear algebra, parallelization, and development of specialized block samplers can further increase efficiency and model complexity.

Conclusion

stLMM advances ecological SAE by operationalizing a comprehensive Bayesian linear mixed modeling workflow over spatial, temporal, and space-time supports, with support for complex direct and unit-level models, structured residuals, and posterior sampling. Its formula-based interface and collapsed sparse-precision machinery make explicit, reproducible modeling accessible to practitioners needing credible inference and reporting in under-sampled ecological domains. The package’s contributions are particularly significant for workflows involving sensitivity analysis, benchmarking, and reporting policy-relevant ecological quantities with uncertainty.

Citation: "stLMM: Bayesian Spatial and Space-Time Linear Mixed Models for Small-Area Ecological Estimation" (2607.02836)