Recover latent process draws

Recover posterior draws of latent process terms from a fitted collapsed stLMM model.

Usage

recover(object, ...)

# S3 method for class 'stLMM'
recover(object, n_omp_threads = NULL, verbose = FALSE,
  sub_sample = list(start = 1L, thin = 1L), ...)

# S3 method for class 'stLMM_chains'
recover(object, n_omp_threads = NULL, verbose = FALSE,
  sub_sample = list(start = 1L, thin = 1L), ...)

Arguments

object

An object returned by stLMM(). Multi-chain fits are handled chain by chain.

n_omp_threads

Optional positive integer number of OpenMP threads for threaded recovery setup kernels. The default NULL uses the thread count stored in the fitted object. This mainly affects NNGP covariance-coefficient updates needed while reconstructing latent process draws; sparse CHOLMOD factorization and random-number generation are not OpenMP-threaded by recover().

verbose

Logical; if TRUE, print phase-level recovery progress messages. Progress is reported from R before and after major recovery phases rather than from individual OpenMP worker threads.

sub_sample

A list with optional integer entries start and thin selecting posterior draws to recover or select from saved process draws.

...

Currently unused.

Details

For Gaussian process models, the collapsed sampler integrates latent process values out during fitting. This function reconstructs the posterior conditional distribution of the latent process for selected posterior draws and samples from it using the stored backend, posterior parameter samples, and CHOLMOD sparse factorization.

By default, recovery uses the OpenMP thread count stored in the fitted object. Use n_omp_threads to override that value for a recovery call without refitting. The override is most relevant for NNGP terms because their covariance-coefficient updates are threaded. Other major recovery work, including sparse factorization and latent Gaussian random draws, is not parallelized by this argument.

For Polya-Gamma process models, process draws are saved during stLMM() by default as save_process = list(start = 1, thin = 1). In that case, recover() only selects and labels the saved in-chain process draws. If a Polya-Gamma process model was fit with save_process = FALSE, recover() errors and the model must be refit with save_process = TRUE or save_process = list(start = ..., thin = ...).

Recovery is only needed for models with structured process terms such as ar1(), gp(), nngp(), car(), dagar(), or car_time(), or dagar_time(). Explicit iid grouped random effects are sampled during fitting and do not need this recovery step. When the fitted data include rows with missing numeric responses, recovered process samples are on the full latent support, including latent nodes associated only with those missing-response rows.

For Polya-Gamma process models, the sampler uses warm-running internal Polya-Gamma latent-process draws to update the augmentation variables. Those same in-chain process draws are the process samples used by recover(). Post-fit Polya-Gamma reconstruction is intentionally not supported.

Value

The fitted object with recovered latent process draws added. User-facing w_samples are returned by process term; NNGP and DAGAR terms are in original support order. The internal ordered copy used by fitted values and prediction is stored in w_samples_ordered. The object also includes w_samples_stacked and recover_iter. The returned object has class stLMM_recovery and can be passed to fitted() or predict(). For multi-chain input, the returned object has class stLMM_recovery_chains and contains chains, a list of ordinary stLMM_recovery objects. The sub_sample rule is applied within each chain.

See also

stLMM, fitted.stLMM, predict.stLMM_recovery, summary.stLMM_recovery

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