Predict Fine-Support Quantities — cos_predict

Predicts fine-support effects and latent means. For spatial fits this uses the conditional Gaussian distribution of fine-support effects given observed-support effects. For non-spatial fits the spatial effect is fixed at zero and predictions use recovered regression coefficients only.

Usage

cos_predict_fine(fit, rec_B, pred_coords = NULL, X_pred = NULL,
  offset_pred = NULL,
  spatial_uncertainty = c("conditional_mean", "marginal", "joint"),
  target = c("latent", "observed"), keep_samples = FALSE,
  n_threads = NULL, verbose = TRUE)

Arguments

fit: A "cos_fit" object.
rec_B: A "cos_recovery_B" object.
pred_coords: Optional matrix of prediction coordinates with columns x and y. If both this and X_pred are omitted, uses the complete raster cells intersecting observed polygons stored in fit. Supply both arguments to predict over the full raster domain or elsewhere.
X_pred: Optional matrix of prediction covariates with column names matching the fitted covariates.
offset_pred: Optional prediction offset for Polya-Gamma response families. Supply NULL or a finite numeric vector of length one or nrow(X_pred). For negative-binomial prediction this is commonly a log exposure. NULL means zero; fitted offsets are not automatically reused.
spatial_uncertainty: How to handle unobserved fine-support spatial effects. "conditional_mean" is faster but omits conditional spatial prediction uncertainty. "marginal" draws independent marginal conditional effects for each prediction cell. "joint" draws the joint conditional spatial effect across fine prediction cells and can require a large dense covariance matrix.
target: Prediction target. For Gaussian fits, "latent" returns latent fine-support means eta_A = X_A beta + omega_A, while "observed" also adds independent fine-support nugget draws with variance tau_sq. For negative-binomial fits, "latent" summarizes the log mean eta and mean count mu = exp(eta), while "observed" also simulates count draws from the fitted negative-binomial likelihood.
keep_samples: Logical; if TRUE, store posterior sample matrices and add 2.5%, 50%, and 97.5% quantiles to the per-cell summary. If FALSE, return means and standard deviations without quantiles.
n_threads: Optional number of OpenMP threads used by compiled covariance code. If NULL, uses fit$prep$n_threads.
verbose: Logical; if TRUE, print progress messages.

Details

By default the function reports the latent fine-support mean target. This is the usual target for mapping the underlying forest condition. Use target = "observed" when the desired posterior predictive target is a fine-support response with finite-cell nugget variation included.

With spatial_uncertainty = "conditional_mean", posterior summaries reflect uncertainty in recovered parameters and observed-support effects, but they omit conditional uncertainty in unobserved fine-support spatial effects. Use spatial_uncertainty = "marginal" for independent per-cell marginal posterior predictive uncertainty or spatial_uncertainty = "joint" when joint fine-surface samples are needed.

With spatial_uncertainty = "joint", the conditional covariance can be positive semidefinite when predicting fine cells constrained by observed-support effects. In that case the function uses eigen-based sampling for the valid degenerate Gaussian distribution.

For binomial fits, the function reports link-scale eta_* summaries and response probability p_* summaries. For negative-binomial fits, it reports link-scale eta_* summaries and mean-count mu_* summaries. When target = "observed" for a negative-binomial fit, it also reports simulated count y_* summaries.

Value

A list with class "cos_prediction_fine" containing prediction summaries and, optionally, posterior samples. Existing eta_* summaries are always returned. When target = "observed", y_* summaries and, if requested, y_samples are also returned. Binomial fits return p_* summaries and optional p_samples; negative-binomial fits return mu_* summaries and optional mu_samples. The object also contains timing, a named numeric vector with user, system, and elapsed seconds.

Examples

if (FALSE) { # \dontrun{
library(sf)
library(raster)
library(ggplot2)

chm <- raster(system.file("extdata", "example_chm.tif", package = "cosTapered"))
names(chm) <- "chm"
B_sf <- st_read(system.file("extdata", "example_B.gpkg",
                            package = "cosTapered"), quiet = TRUE)
forest_sf <- st_read(system.file("extdata", "example_forest.gpkg",
                                 package = "cosTapered"), quiet = TRUE)
truth <- readRDS(system.file("extdata", "example_truth.rds",
                             package = "cosTapered"))

intercept <- chm
intercept[] <- ifelse(is.na(getValues(chm)), NA_real_, 1)
names(intercept) <- "intercept"
X_rast <- stack(intercept, chm)
names(X_rast) <- c("intercept", "chm")

prep <- cos_prepare(X_rast, B_sf, response_col = "y_B",
                    gamma = truth$gamma, taper_code = truth$taper_code,
                    verbose = FALSE)
priors <- cos_default_priors(
  prep,
  beta_mu = c(intercept = 40, chm = 1),
  beta_sd = c(intercept = 20, chm = 2),
  tau_B_shape = 2,
  tau_B_scale = truth$tau_B_sq_true,
  sigma_shape = 2,
  sigma_scale = truth$sigma_sq_true,
  phi_lower = 3 / 600,
  phi_upper = 3 / 80
)
fit <- cos_fit(
  prep, priors,
  starting = data.frame(tau_B_sq = truth$tau_B_sq_true,
                        sigma_sq = truth$sigma_sq_true,
                        phi = truth$phi_true),
  tuning = c(log_tau_B_sq = 0.25, log_sigma_sq = 0.25, z_phi = 0.25),
  n_batch = 50,
  batch_length = 8,
  seed = 1,
  verbose = FALSE
)
rec_B <- cos_recover_B(fit, burn_in = 100, n_samples = 100,
                       seed = 1, verbose = FALSE)

pred <- cos_predict_fine(fit, rec_B,
                         spatial_uncertainty = "conditional_mean",
                         verbose = FALSE)
summary(pred)

fine_df <- pred$summary[, c("x", "y", "eta_mean", "eta_sd")]

ggplot() +
  geom_tile(data = fine_df, aes(x = x, y = y, fill = eta_mean),
            width = xres(chm), height = yres(chm)) +
  geom_sf(data = B_sf, fill = NA, color = "black", linewidth = 0.25) +
  geom_sf(data = forest_sf, fill = NA, color = "white", linewidth = 0.35) +
  scale_fill_viridis_c(option = "magma") +
  coord_sf(expand = FALSE) +
  labs(title = "Fine-support posterior mean",
       fill = "eta mean", x = NULL, y = NULL) +
  theme_bw()

ggplot() +
  geom_tile(data = fine_df, aes(x = x, y = y, fill = eta_sd),
            width = xres(chm), height = yres(chm)) +
  geom_sf(data = B_sf, fill = NA, color = "black", linewidth = 0.25) +
  geom_sf(data = forest_sf, fill = NA, color = "white", linewidth = 0.35) +
  scale_fill_viridis_c(option = "cividis") +
  coord_sf(expand = FALSE) +
  labs(title = "Fine-support posterior SD",
       fill = "eta SD", x = NULL, y = NULL) +
  theme_bw()
} # }