Posterior prediction from recovered stLMM fits

Returns posterior samples for the mean response at fitted rows or at newdata rows. For binomial models, the default mean response is the posterior fitted probability. For negative-binomial models, the default mean response is the posterior fitted mean count.

Usage

# S3 method for class 'stLMM'
predict(object, newdata = NULL, y_samples = FALSE,
  n_omp_threads = NULL, verbose = FALSE, sub_sample = list(start = 1L,
  thin = 1L), scale = c("response", "link"), ...)

# S3 method for class 'stLMM_chains'
predict(object, newdata = NULL, y_samples = FALSE,
  n_omp_threads = NULL, verbose = FALSE, sub_sample = list(start = 1L,
  thin = 1L), scale = c("response", "link"), ...)

# S3 method for class 'stLMM_recovery'
predict(object, newdata = NULL, y_samples = FALSE,
  joint = FALSE, joint_method = c("full", "vecchia"), pred_m = NULL,
  pred_ordering = "maxmin", st_scale = NULL,
  return_w_samples = TRUE, n_omp_threads = NULL, verbose = FALSE,
  sub_sample = list(start = 1L, thin = 1L), scale = c("response", "link"),
  ...)

# S3 method for class 'stLMM_recovery_chains'
predict(object, newdata = NULL,
  y_samples = FALSE, joint = FALSE, joint_method = c("full", "vecchia"),
  pred_m = NULL, pred_ordering = "maxmin", st_scale = NULL,
  return_w_samples = TRUE, n_omp_threads = NULL, verbose = FALSE,
  sub_sample = list(start = 1L, thin = 1L), scale = c("response", "link"),
  ...)

Arguments

object

An stLMM fit without structured process terms, a Polya-Gamma process fit with saved in-chain process draws, a recovered object returned by recover(), or the corresponding multi-chain object.

newdata

Optional data frame. Grouped random-effect levels must already exist in the fitted model support. For recovered process models, process coordinates may be existing support locations or new AR1, GP, or NNGP nodes. CAR and CAR-time prediction rows must use existing graph areas; CAR-time rows must also use fitted time values.

y_samples

Logical; if TRUE, also simulate observation-level posterior predictive samples. Gaussian fits use the fitted residual variance. Binomial fits draw counts from the fitted probabilities and prediction trial counts.

joint

Logical; for new GP or NNGP nodes, whether to draw unique new nodes with cross-location dependence. The default FALSE simulates unique new nodes independently conditional on the fitted support or fitted-support neighbors and is recommended when only marginal prediction summaries are needed. For dense gp() terms, joint = TRUE uses the exact dense GP conditional covariance among unique new nodes.

joint_method

Method used for NNGP terms when joint = TRUE. "full" keeps the historical behavior: fitted-support neighbors are used for each new node and the new-node residual covariance is drawn from a dense matrix. "vecchia" orders the unique new prediction nodes and simulates them sequentially from a Vecchia/NNGP-style history containing all fitted support nodes plus previously simulated prediction nodes. Dense gp() terms ignore joint_method and use exact joint Gaussian conditioning when joint = TRUE.

pred_m

Optional positive integer neighbor count for joint_method = "vecchia". The default NULL uses the fitted NNGP term's m. The first prediction nodes always have the fitted support available; if fewer than pred_m history nodes are available, all available history nodes are used.

pred_ordering

Ordering for unique new prediction nodes under joint_method = "vecchia". May be one of the NNGP ordering options "coord", "default", "maxmin", "hilbert", "random", or a numeric permutation of the unique new prediction nodes. For space-time NNGP terms, "hilbert" is currently rejected because the implemented Hilbert ordering is two-dimensional.

st_scale

NULL, a positive scalar, or a named positive vector/list by process term. For space-time NNGP new-node prediction, NULL inherits the st_scale stored on the fitted NNGP graph, falling back to 1 for older objects. A supplied value overrides the fitted scale. The final coordinate is multiplied by st_scale before prediction ordering and nearest-neighbor ranking. This is not automatic standardization: spatial coordinates are left on their original scale and the time coordinate is left on its original scale except for this single multiplier during graph construction. The covariance calculation still uses the original coordinates. Non-default st_scale values are valid only for space-time NNGP terms.

return_w_samples

Logical; for recovered process prediction, whether to return raw latent process samples for each process term in w_samples. The default TRUE is useful for diagnostics and process maps. Set to FALSE for large prediction problems when memory use matters. For covariate-scaled process terms, returned values are the raw latent process draws, not the covariate-multiplied row-level contribution.

n_omp_threads

Optional positive integer number of OpenMP threads for threaded prediction kernels. The default NULL uses the thread count stored in the fitted or recovered object. This affects threaded NNGP prediction neighbor searches and compiled NNGP prediction simulation kernels; R-level dense joint GP and dense joint NNGP prediction are not OpenMP-threaded.

verbose

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

sub_sample

A list with optional integer entries start and thin. For predict.stLMM() without process terms, these select posterior draw indices. For Polya-Gamma process fits with saved process draws, and for predict.stLMM_recovery(), these select saved or recovered draw indices through object$recover_iter.

scale

For family = "binomial", return fitted probabilities on the "response" scale or the linear predictor on the "link" scale. For family = "negative_binomial", return fitted mean counts on the "response" scale or the log mean on the "link" scale. For Gaussian models these scales are identical.

...

Currently unused.

Details

Prediction always represents the model as fit. Prediction-time dropping of model components is not supported. New grouped random-effect levels are not supported.

By default, prediction uses the OpenMP thread count stored in the fitted or recovered object. Use n_omp_threads to override that value for a prediction call without refitting or re-running recover(). Because the long-running prediction kernels can be OpenMP threaded, verbose = TRUE reports progress at the R phase level rather than printing from worker threads.

With newdata = NULL, prediction returns fitted-row posterior mean samples for the original data rows. With newdata, fixed-effect and iid random-effect design matrices are rebuilt from the fitted formula metadata. If the fitted formula contains offset(), the same offset expression is evaluated in newdata and added to the prediction linear predictor. New AR1, GP, and NNGP process nodes are supported. Repeated prediction rows sharing a new process node reuse the same simulated latent draw within a posterior sample.

When return_w_samples = TRUE, recovered-process prediction objects also contain w_samples, a named list with one matrix per process term. Each matrix has posterior draws in rows and prediction rows in columns, aligned with mu_samples. These are raw latent process values on the prediction rows. For spatially varying coefficient terms, w_samples does not include the covariate scaling used when adding the process contribution to mu_samples.

For family = "binomial", prediction first constructs posterior samples of the linear predictor and then applies the inverse-logit transform when scale = "response". If y_samples = TRUE, prediction simulates binomial counts. With newdata = NULL, the fitted trial counts are used. With newdata, a column named trials is used when present; otherwise prediction defaults to one trial per row.

For family = "negative_binomial", prediction first constructs posterior samples of the log mean and then exponentiates them when scale = "response". If y_samples = TRUE, prediction simulates negative-binomial counts using the fitted fixed size.

For dense gp() terms, joint = FALSE draws each unique new node from its marginal conditional distribution given the fitted support. With joint = TRUE, unique new GP nodes are drawn jointly from the exact conditional Gaussian distribution using the dense covariance among prediction nodes and fitted nodes. This is appropriate for small prediction sets and is dense in the number of unique new GP nodes.

For NNGP terms, joint = FALSE uses the independent new-node prediction approach based on nearest fitted-support neighbors, following the conditional NNGP prediction machinery described by Finley et al. (2019). Each unique new node is simulated independently conditional on its fitted-support neighbors. With joint = TRUE, joint_method = "full", nearest neighbors are still chosen only from the fitted support, but unique new NNGP nodes are drawn jointly from a dense residual covariance. This can require substantial memory for large prediction sets.

With joint = TRUE, joint_method = "vecchia", unique new NNGP nodes are ordered by pred_ordering. Neighbors for each ordered prediction node are chosen from the combined history made of all fitted support nodes and earlier prediction nodes in that ordering. The prediction draw is then generated sequentially, so simulated earlier prediction nodes can enter the conditional mean of later prediction nodes. This avoids a dense new-node covariance matrix and is usually the more scalable joint prediction method. This option follows the observed-prediction ordering idea and latent-first full-conditioning construction described by Katzfuss et al. (2020, Sects. 3.1, 4.2, and 4.3.4), adapted to the package workflow in which fitted latent process draws are saved or recovered first and prediction nodes are then simulated sequentially conditional on those draws. The method is ordering dependent, as in other Vecchia/NNGP approximations.

The independent, dense-joint, and Vecchia-joint approaches should have similar posterior means and medians for each prediction location when the neighbor sets are adequate. They differ mainly in the cross-location dependence of the simulated prediction draws.

For space-time NNGP prediction, st_scale controls the search geometry used to order new nodes and choose nearest neighbors. The search distance is equivalent to Euclidean distance after replacing the final coordinate by st_scale * time. Thus st_scale = 1 uses the coordinates as supplied, while st_scale = s treats one unit of the final coordinate as s spatial-coordinate units for graph construction. Larger values make temporal separation matter more and tend to choose neighbors closer in time; smaller values make temporal separation matter less. This choice does not change the covariance parameters or covariance evaluations once the neighbors have been selected.

References

Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019). Efficient Algorithms for Bayesian Nearest Neighbor Gaussian Processes. Journal of Computational and Graphical Statistics, 28(2), 401–414. doi:10.1080/10618600.2018.1537924 .

Katzfuss, M., Guinness, J., Gong, W., and Zilber, D. (2020). Vecchia Approximations of Gaussian-Process Predictions. Journal of Agricultural, Biological and Environmental Statistics, 25(3), 383–414. doi:10.1007/s13253-020-00401-7 .

Value

An object of class stLMM_prediction, a list with mu_samples, optional y_samples, optional w_samples, draw_index, and metadata. Binomial prediction objects also include scale. For multi-chain input, the returned object has class stLMM_prediction_chains and contains chains, a list of ordinary stLMM_prediction objects. The sub_sample rule is applied within each chain.

See also

recover, fitted.stLMM, summary.stLMM_prediction

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)
pred_fit <- predict(rec)
newdat <- data.frame(x = c(-0.5, 0.5), time = c(9, 10))
pred_new <- predict(rec, newdata = newdat, y_samples = TRUE)
summary(pred_new)
#> stLMM prediction summary
#>   draws: 8
#>   rows: 2
#>   newdata: TRUE
#>   process samples: ar1_1
#> 
#> mu:
#>      mean     sd    q2.5   q50.0  q97.5
#> 1 -0.1472 0.4973 -0.5674 -0.3293 0.6803
#> 2 -0.2663 0.9015 -1.4613 -0.1272 0.8418
#> 
#> y:
#>      mean     sd    q2.5   q50.0  q97.5
#> 1 -0.6677 1.1030 -1.8275 -0.9740 1.1862
#> 2 -0.6035 1.0926 -1.5325 -0.8202 1.5008