Fit a Bayesian mixed model with structured latent processes

Fits Bayesian linear mixed models with fixed effects, optional explicit iid grouped random effects, and optional structured latent process terms. Structured process terms are integrated out during parameter updates. For Gaussian process models, their latent draws are recovered afterward with recover. For Polya-Gamma process models, the in-chain process draws already used by the augmentation step are saved and reused by fitted-value, recovery, and prediction methods.

For family = "binomial", the sampler uses Polya-Gamma latent variables. Conditional on the Polya-Gamma weights, the model is fit as a Gaussian working model with diagonal observation precision equal to the current weights. Polya-Gamma random variates are generated through BayesLogit's hybrid sampler. Binomial trials default to one per row and can be supplied through trials. The Gaussian resid() residual variance term and tau_sq controls are not used for binomial fits. For structured binomial process models, the sampler saves in-chain process draws by default; recover selects from those saved draws.

For family = "negative_binomial", the response must be non-negative integer counts and size must be supplied as a positive fixed scalar. The sampler uses the log-mean parameterization $\eta_i=\log(\mu_i)$ and an internal logit tilt $\psi_i=\eta_i-\log(\code{size})$. Conditional on the Polya-Gamma weights, the model again becomes a Gaussian working model with diagonal observation precision equal to the current weights. Residual resid() terms and tau_sq controls are not used for this likelihood. Larger size values make the negative-binomial distribution closer to Poisson, but can slow the sampler because the Polya-Gamma shape is the observed count plus size.

Usage

stLMM(formula, data = parent.frame(), starting = list(), tuning = NULL,
  priors = NULL, family = "gaussian", trials = NULL,
  size = NULL, n_samples, n_omp_threads = 1, nngp_search = c("fast", "brute"),
  verbose = TRUE, n_report = 100, warmup = TRUE, metropolis = list(),
  cholmod_control = list(), save_process = NULL, chains = 1L,
  chain_control = list(), describe_terms = FALSE, ...)

Arguments

formula: A model formula with fixed effects, optional offset() terms, optional iid grouped random-effect terms such as iid(group) or x:iid(group), and optional structured process terms such as ar1(), gp(), nngp(), car(), dagar(), car_time(), or dagar_time().
data: A data frame or environment containing variables in formula.
starting: Optional named list of starting values. The public convention uses underscores. Default residual variance is supplied as starting = list(resid = list(tau_sq = value)). Process starts use term-specific blocks, for example starting = list(nngp_1 = c(sigma_sq = 1, phi = 1)). IID random-effect variance starts use the same term-block convention, for example starting = list(iid_1 = list(sigma_sq = 1)). Covariance and residual variance parameters can be fixed with fixed, for example starting = list(resid = list(tau_sq = fixed(0.5))).
tuning: Optional named list of initial Metropolis proposal scales for covariance and residual variance parameters. A zero tuning value fixes the corresponding parameter. Parameters supplied with fixed are assigned zero tuning automatically.
priors: Named list of prior objects created with constructors such as flat, normal, ig, uniform, log_normal, gamma_dist, half_normal, half_t, and beta_dist. Fixed effects use priors = list(beta = flat()) or priors = list(beta = normal(mean, sd)). Default residual variance uses resid$tau_sq; for example, priors = list(resid = list(tau_sq = ig(shape, scale))) puts an inverse-gamma prior directly on tau_sq, while priors = list(resid = list(tau_sq = half_t(df, scale))) defines a half-$t$ prior on the residual standard deviation $\sqrt{\tau^2}$ and transforms it internally to the variance scale. Process terms require term-specific blocks, for example priors = list(nngp_1 = list(sigma_sq = ig(shape, scale), phi = uniform(lower, upper))). IID random-effect variance priors use the same term-block convention, for example priors = list(iid_1 = list(sigma_sq = ig(shape, scale))). Priors are not required for covariance or residual variance parameters supplied with fixed.
family: Model likelihood family. Use "gaussian" for Gaussian linear mixed models or "binomial" for binomial logistic models using Polya-Gamma Gibbs updates. Use "negative_binomial" for fixed-size negative-binomial log-mean count models using Polya-Gamma Gibbs updates. The default is "gaussian". For compatibility with familiar R model syntax, gaussian() and binomial() are also accepted.
trials: Optional binomial trial counts for family = "binomial". The default NULL uses one trial per row, i.e. Bernoulli data. A single character string names a column in data; otherwise supply a positive integer-valued vector with length matching the model frame. The response should contain successes and must be between zero and trials.
size: Positive fixed negative-binomial size parameter for family = "negative_binomial". The fitted linear predictor is $\eta_i = \log(\mu_i)$, where $\mu_i$ is the mean count.
n_samples: Number of MCMC samples.
n_omp_threads: Number of OpenMP threads for supported kernels.
nngp_search: Nearest-neighbor search used when constructing NNGP graphs. "fast" uses the exact history-restricted k-d tree implementation. "brute" checks every eligible historical node and is mainly useful for testing and debugging.
verbose: Logical; print setup, sampler-start, and progress reports.
n_report: Progress report interval for retained sampling and warmup. For warmup, reports are rounded to completed warmup batches. Use 0 to suppress progress reports while keeping other verbose setup messages.
warmup: Logical or list controlling discarded Metropolis proposal-scale calibration before retained sampling. TRUE uses default calibration; FALSE disables it. A list may contain enabled, batch_length, min_batches, max_batches, target, and near_zero. min_batches forces at least that many warmup batches before early stopping is allowed; max_batches remains the ceiling. Defaults are enabled = TRUE, batch_length = 25, min_batches = 0, max_batches = 20, and near_zero = 0.02. When target is omitted, the default acceptance interval is c(0.15, 0.45) for block updates and c(0.30, 0.60) for scalar updates.
metropolis: Optional sampler-control list for Metropolis updates of covariance and residual variance parameters. The default list(blocking = "joint") keeps all active parameters in one adaptive covariance block. Built-in alternatives are "by_term", "residual_process", "variance_theta", "process_variance", and "scalar". A single string such as metropolis = "by_term" is also accepted. A list may also contain target_accept for retained-sampling adaptation and batch_length for retained adaptation batches. Defaults are target_accept = 0.234 for block updates, 0.44 for scalar updates, and batch_length = 25.
cholmod_control: Optional sparse factorization control list. The ordering element controls the CHOLMOD fill-reducing ordering used for latent sparse matrix factorizations. Supported values are "auto" (the default CHOLMOD strategy), "best" (try CHOLMOD's full built-in method set), "natural", "amd", "metis", "nesdis", and "colamd". A single string such as cholmod_control = "amd" is also accepted. The postorder element is logical and defaults to TRUE. METIS-based choices require the CHOLMOD library exposed by Matrix to have been built with METIS/partitioning support.
save_process: Controls whether structured latent process draws are saved during MCMC. Use NULL for the default: process draws are saved automatically for Polya-Gamma process models as list(start = 1, thin = 1) and not saved otherwise. Use TRUE as shorthand for list(start = 1, thin = 1), use FALSE to turn saving off, or use a list with start and thin, for example save_process = list(start = 501, thin = 5). For Polya-Gamma process models, saved in-chain process draws are required for recover(), fitted(), and process prediction. Enabling save_process is only supported for Polya-Gamma models with structured process terms.
chains: Positive integer number of MCMC chains to run. The default 1L returns an ordinary stLMM object. Values greater than one run sequential chains and return an stLMM_chains object.
chain_control: Optional list controlling multi-chain initialization. seed initializes reproducible chain-specific random number streams when supplied. dispersion controls log-scale jitter for omitted positive starting values.
describe_terms: Logical; print model term diagnostics before sampling.
...: Currently unused.

Value

With chains = 1, an object of class stLMM. With chains > 1, an object of class stLMM_chains with component chains, a list of ordinary stLMM objects, one per chain.

For each ordinary fit, posterior samples are available both as top-level entries such as beta_samples, tau_sq_samples, sigma_sq_samples, theta_samples, and alpha_samples, and through the active-component list samples. Empty sample blocks for absent model components are NULL. Covariance Metropolis diagnostics are returned as covariance_acceptance and adaptive_metropolis. The default is one global covariance block; optional blocking and scalar-update diagnostics are stored in adaptive_metropolis$blocks. The object also contains timing, with CPU and elapsed wall time for the compiled sampler call, term_description, and the stored backend used by recovery and prediction. Multi-chain fits also include top-level sampler timing by chain and total sampler timing across chains.

Latent process draws are saved during fitting only for Polya-Gamma likelihoods with structured process terms. The default in that case is save_process = list(start = 1, thin = 1), because those in-chain process draws are already needed by the Polya-Gamma augmentation and post-fit Polya-Gamma process reconstruction is not used. For Gaussian process models, the sampler remains fully collapsed with respect to the process terms and recover() is the intended post-fit step for drawing latent process samples. Multi-chain fits can be converted to a coda::mcmc.list with as_mcmc and summarized with summary(), which includes chain diagnostics when possible.

Details

Predictor variables must be complete. Missing numeric responses are allowed: only observed response rows enter the likelihood, while the full row support is retained for fitted values, recovery, and prediction metadata. The number of missing responses is reported by print() and summary(). Formula offsets are supported through standard R offset() terms. The offset is a known additive component of the linear predictor. For negative-binomial count models, offset(log_exposure) gives $\log(\mu_i) = \log(\mathrm{exposure}_i) + \cdots$; for binomial models it is a known log-odds contribution. Fitted values and predictions include the offset on both link and response scales. When predicting with newdata, variables used inside offset() must be present in newdata.

The intended workflow is:


fit <- stLMM(...)
summary(fit)
rec <- recover(fit)
fitted(rec)
predict(rec, newdata = ...)

For models without structured process terms, fitted() and predict() can be called directly on the fit object.

Prior lists use explicit prior constructors. Raw numeric vectors such as c(2, 1) are not accepted as priors. Common patterns are:


## residual variance
priors = list(resid = list(tau_sq = ig(shape, scale)))
priors = list(resid = list(tau_sq = half_t(df = 3, scale = 0.5)))

## fixed effects
priors = list(beta = flat())
priors = list(beta = normal(mean = 0, sd = 10))

## known additive offset
y ~ x + offset(log_exposure) + car(area, graph = g)

## fixed direct-estimate sampling variances
y ~ x + resid(model = "fixed", variance = vhat)

## sampled group-specific residual variances
y ~ x + resid(model = "group", group = area)
priors = list(resid = list(tau_sq = half_t(df = 3, scale = 1)))

## direct-estimate-informed group-specific residual variances
y ~ x + resid(model = "group", group = area, variance = vhat,
              prior = "ig", shape = 6)

## scaled direct-estimate sampling variances
y ~ x + resid(model = "scaled", variance = vhat, n = n_eff)

## iid random-effect variance for iid(group)
priors = list(iid_1 = list(sigma_sq = ig(shape, scale)))

## fixed covariance parameters
starting = list(
  resid = list(tau_sq = fixed(0.5)),
  nngp_1 = list(sigma_sq = fixed(1), phi = fixed(3))
)

## process covariance controls
priors = list(nngp_1 = list(sigma_sq = ig(shape, scale),
                            phi = uniform(lower, upper)))

## CAR controls
priors = list(car_1 = list(sigma_sq = ig(shape, scale),
                           rho = uniform(lower, upper)))

## DAGAR controls
priors = list(dagar_1 = list(sigma_sq = ig(shape, scale),
                             rho = uniform(lower, upper)))

## separable CAR-time controls
priors = list(car_time_1 = list(sigma_sq = ig(shape, scale),
                                rho = uniform(lower, upper),
                                phi = uniform(lower, upper)))

## separable CAR-time controls with time_model = "exp"
priors = list(car_time_1 = list(sigma_sq = ig(shape, scale),
                                rho = uniform(lower, upper),
                                lambda = uniform(lower, upper)))

## separable DAGAR-time controls
priors = list(dagar_time_1 = list(sigma_sq = ig(shape, scale),
                                  rho = uniform(lower, upper),
                                  phi = uniform(lower, upper)))

Starting values, tuning, and fixed parameters. Omitting starting uses data-based or conservative defaults. Omitting tuning uses default initial Metropolis proposal scales. Setting a tuning value to zero fixes that covariance or residual variance parameter. Parameters marked with fixed() use their supplied values and are omitted from Metropolis proposal blocks.
Priors and parameter transforms. Process terms require priors for parameters that are sampled. The Gaussian likelihood defaults to a flat fixed effect prior, while Polya-Gamma likelihoods default to normal(mean = 0, sd = 10) for fixed effects. The flat() and normal() constructors are accepted only for priors$beta. Variance parameters are positive. The half_normal() and half_t() variance-prior constructors are specified on the corresponding standard deviation scale and transformed internally to the variance scale. Covariance theta parameters use a bounded-logit proposal transform and therefore require finite support. For positive theta parameters, log_normal() and gamma_dist() must include support = c(lower, upper). See stLMM_priors and stLMM_terms for details.
Multiple chains. For chains > 1, omitted starting values are dispersed automatically from the defaults. With chain_control$seed, chain-specific random number streams are generated from the supplied seed, so dispersed starts are reproducible across different Metropolis blocking choices. If a starting value is supplied by the user, it controls that parameter: scalar values are recycled to all chains, and vectors with length chains provide chain-specific starts. Other lengths are rejected. For process terms, chain-specific starts should be supplied in list blocks, for example starting = list(nngp_1 = list(sigma_sq = c(1, 2, 3), phi = c(0.2, 0.4, 0.6))).
Warmup and retained-sampling adaptation. By default, the Metropolis sampler runs a short discarded warmup phase that calibrates proposal scales before retained samples are collected. Warmup draws are not returned in posterior sample matrices; warmup diagnostics are stored in fit$adaptive_metropolis$warmup. Retained-sampling adaptation is controlled separately through metropolis$target_accept and metropolis$batch_length.
Metropolis blocking. Metropolis blocking changes only the proposal partition used for covariance and residual variance parameters; it does not change the model. Blocking can improve mixing for multi-term or residual-variance models, but each block may require a separate collapsed likelihood evaluation and sparse factorization. The default "joint" block is usually the fastest choice. Treat the other blocking schemes as diagnostic or model-specific alternatives: compare trace plots, R-hat, effective sample sizes, and runtime before changing the default for routine analyses. metropolis = "scalar" updates one active covariance or residual variance parameter at a time using scalar random-walk adaptation. It can be useful when block proposals mix poorly, but it may be substantially slower for models with expensive collapsed likelihood evaluations.
Sparse factorization controls. The cholmod_control argument is intended for advanced diagnostics and performance benchmarking. The default ordering = "auto" uses CHOLMOD's default ordering strategy. Alternative orderings can change sparse Cholesky fill-in, floating-point work, and runtime without changing the statistical model. Use describe_terms = TRUE or inspect fit$term_description to compare the selected CHOLMOD ordering, fill ratio, lnz, and flops. The "best" option can make symbolic analysis slower because CHOLMOD tries more orderings before factorization.
Fixed residual variances. When resid(model = "fixed", variance = vhat) is used, tau_sq is not sampled and should not be supplied in starting, tuning, or priors.

See stLMM_terms for the package-native formula terms recognized by stLMM().

Examples

set.seed(1)
dat <- data.frame(
  y = rnorm(8),
  x = rnorm(8),
  time = seq_len(8)
)

fit <- stLMM(
  y ~ x + ar1(time),
  data = dat,
  n_samples = 8,
  priors = list(
    resid = list(tau_sq = ig(2, 1)),
    ar1_1 = list(sigma_sq = ig(2, 1), phi = uniform(-0.8, 0.8))
  ),
  warmup = FALSE,
  verbose = FALSE
)
rec <- recover(fit, sub_sample = list(start = 5, thin = 1))
fitted(rec)
#> [1] -0.66560060  0.36628335 -0.42578945  1.11258894  0.06457469 -0.28426435
#> [7]  0.85185568  1.03213562