| Title: | Spatial Analysis with Misaligned Data Using Atom-Based Regression Models |
|---|---|
| Description: | Implements atom-based regression models (ABRM) for analyzing spatially misaligned data. Provides functions for simulating misaligned spatial data, preparing NIMBLE model inputs, running MCMC diagnostics, and providing results. All main functions return S3 objects with print(), summary(), and plot() methods for intuitive result exploration. Methods originally described in Mugglin et al. (2000) <doi:10.1080/01621459.2000.10474279>, further investigated in Trevisani & Gelfand (2013), and applied in Nethery et al. (2023) <doi:10.1101/2023.01.10.23284410>. |
| Authors: | Yunzhe Qian [aut, cre], Rachel Nethery [aut], Nancy Krieger [ctb] (Contributed to the project conceptualization and manuscript), Nykesha Johnson [ctb] (Contributed to the project conceptualization and manuscript) |
| Maintainer: | Yunzhe Qian <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.2.8 |
| Built: | 2026-05-23 08:55:32 UTC |
| Source: | https://github.com/bellayqian/spatialatomizer |
Returns the posterior mean estimates for all regression coefficients as a named numeric vector.
## S3 method for class 'abrm' coef(object, ...)## S3 method for class 'abrm' coef(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
A named numeric vector of posterior mean estimates. Names correspond
to parameter labels (e.g., "beta_x[1]", "beta_y[1]",
"beta_0_y").
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) coef(results)sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) coef(results)
Returns the 95% posterior credible intervals for all estimated parameters. Note: these are Bayesian credible intervals from MCMC quantiles, not frequentist confidence intervals.
## S3 method for class 'abrm' confint(object, parm = NULL, level = 0.95, ...)## S3 method for class 'abrm' confint(object, parm = NULL, level = 0.95, ...)
object |
An object of class |
parm |
Optional character vector of parameter names to subset. If
|
level |
Confidence level for the credible interval (default |
... |
Additional arguments (currently unused). |
A matrix with two columns, CI.Lower and CI.Upper, and
one row per parameter.
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) confint(results) confint(results, parm = "beta_x[1]")sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) confint(results) confint(results, parm = "beta_x[1]")
Computes posterior-mean fitted values at the Y-grid level by applying the estimated regression coefficients to the supplied covariate data.
## S3 method for class 'abrm' fitted(object, ...)## S3 method for class 'abrm' fitted(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
A data frame with one row per Y-grid cell and four columns:
y_grid_idThe original Y-grid cell ID.
observedThe observed outcome value y for each
Y-grid cell. Note: the outcome is observed at the Y-grid level for
all cells. This is distinct from covariates, which may be
latent (unobserved at the Y-grid level) for cells that span
multiple atoms.
fittedThe model-predicted outcome on the original outcome scale (counts for Poisson/binomial; continuous for normal).
residualObserved minus fitted.
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) fitted(results)sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) fitted(results)
Returns the NIMBLE code for the Atom-Based Regression Model with mixed-type variables. Automatically registers custom distributions if not already registered.
get_abrm_model()get_abrm_model()
A nimbleCode object containing the NIMBLE model specification
for the ABRM. Pass this directly to run_abrm.
model_code <- get_abrm_model() print(model_code)model_code <- get_abrm_model() print(model_code)
Displays MCMC diagnostic plots — trace plots and posterior density plots — for the parameters monitored during model fitting. If no diagnostic plots are stored in the object, a message is issued instead.
## S3 method for class 'abrm' plot(x, ...)## S3 method for class 'abrm' plot(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x. Called for its side effect of rendering
the MCMC diagnostic plots.
## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Display MCMC trace and posterior density plots plot(results)## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Display MCMC trace and posterior density plots plot(results)
Visualizes the spatial layout of a misaligned_data object. By default,
both grids are overlaid to illustrate the misalignment between them.
## S3 method for class 'misaligned_data' plot( x, which = c("both", "gridy", "gridx", "atoms"), color_y = "blue", color_x = "orange", title = NULL, ... )## S3 method for class 'misaligned_data' plot( x, which = c("both", "gridy", "gridx", "atoms"), color_y = "blue", color_x = "orange", title = NULL, ... )
x |
A |
which |
Character string specifying what to plot. One of |
color_y |
Color for the Y-grid boundaries. Default |
color_x |
Color for the X-grid boundaries. Default |
title |
Optional character string for the plot title. If |
... |
Additional arguments passed to |
The input x, invisibly.
## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) # Default: overlay both grids to visualise misalignment plot(sim_data) # Plot a single component plot(sim_data, which = "gridy") plot(sim_data, which = "gridx") plot(sim_data, which = "atoms") # Custom colours plot(sim_data, color_y = "steelblue", color_x = "firebrick") # Custom title plot(sim_data, title = "Utah Spatial Misalignment")## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) # Default: overlay both grids to visualise misalignment plot(sim_data) # Plot a single component plot(sim_data, which = "gridy") plot(sim_data, which = "gridx") plot(sim_data, which = "atoms") # Custom colours plot(sim_data, color_y = "steelblue", color_x = "firebrick") # Custom title plot(sim_data, title = "Utah Spatial Misalignment")
Prints a concise summary of a fitted ABRM, including the number of parameters estimated and, when true parameter values are available, the mean absolute bias and 95% credible-interval coverage rate.
## S3 method for class 'abrm' print(x, ...)## S3 method for class 'abrm' print(x, ...)
x |
An abrm object |
... |
Additional arguments (unused) |
Invisibly returns x. Called for its side effect of printing
a concise ABRM model summary including number of parameters, mean
absolute bias, and credible-interval coverage rate (when true parameter
values are available).
## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM ## (niter and nburnin are intentionally small for illustration only; ## use larger values, e.g. niter = 50000, nburnin = 30000, in practice) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Print a concise model summary print(results) # Reports: number of parameters, mean absolute bias, coverage rate## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM ## (niter and nburnin are intentionally small for illustration only; ## use larger values, e.g. niter = 50000, nburnin = 30000, in practice) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Print a concise model summary print(results) # Reports: number of parameters, mean absolute bias, coverage rate
Print Method for Misaligned Data Objects
## S3 method for class 'misaligned_data' print(x, ...)## S3 method for class 'misaligned_data' print(x, ...)
x |
A misaligned_data object |
... |
Additional arguments (unused) |
Invisibly returns the input object x. The function is called for its side effect of printing a summary of the simulated misaligned spatial data including grid dimensions and number of atoms.
## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) print(sim_data)## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) print(sim_data)
Prints the full parameter table from a summary.abrm object,
including posterior means, standard deviations, credible intervals, and
(when available) bias and coverage.
## S3 method for class 'summary.abrm' print(x, digits = 4, ...)## S3 method for class 'summary.abrm' print(x, digits = 4, ...)
x |
An object of class |
digits |
Integer; number of significant digits. Default |
... |
Additional arguments (currently unused). |
Invisibly returns x.
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) s <- summary(results) print(s) # same as typing summary(results) interactively print(s, digits = 2) # fewer decimal placessim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) s <- summary(results) print(s) # same as typing summary(results) interactively print(s, digits = 2) # fewer decimal places
Prints all variance-covariance matrices stored in a "vcov.abrm"
object returned by vcov.abrm.
## S3 method for class 'vcov.abrm' print(x, ...)## S3 method for class 'vcov.abrm' print(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x. Called for its side effect of printing
the intercept variance and the X-grid and Y-grid covariance matrices.
## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Extract and print the variance-covariance matrices vcov_mats <- vcov(results) print(vcov_mats) # Prints: intercept variance, X-grid and Y-grid covariance matrices## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Extract and print the variance-covariance matrices vcov_mats <- vcov(results) print(vcov_mats) # Prints: intercept variance, X-grid and Y-grid covariance matrices
Displays a formatted WAIC summary for a fitted ABRM.
## S3 method for class 'waic_abrm' print(x, digits = 3, ...)## S3 method for class 'waic_abrm' print(x, digits = 3, ...)
x |
An object of class |
digits |
Integer; number of decimal places. Default |
... |
Additional arguments (currently unused). |
Invisibly returns x.
Fits an Atom-Based Regression Model (ABRM) to spatially misaligned data using Bayesian MCMC via NIMBLE. Works with both simulated and real-world data.
run_abrm( gridx, gridy, atoms, model_code, true_params = NULL, norm_idx_x = NULL, pois_idx_x = NULL, binom_idx_x = NULL, norm_idx_y = NULL, pois_idx_y = NULL, binom_idx_y = NULL, dist_y = 2, niter = 50000, nburnin = 30000, nchains = 2, thin = 10, seed = NULL, sim_metadata = NULL, save_plots = FALSE, output_dir = NULL, compute_waic = FALSE )run_abrm( gridx, gridy, atoms, model_code, true_params = NULL, norm_idx_x = NULL, pois_idx_x = NULL, binom_idx_x = NULL, norm_idx_y = NULL, pois_idx_y = NULL, binom_idx_y = NULL, dist_y = 2, niter = 50000, nburnin = 30000, nchains = 2, thin = 10, seed = NULL, sim_metadata = NULL, save_plots = FALSE, output_dir = NULL, compute_waic = FALSE )
gridx |
The X-grid sf dataframe, containing a numeric area ID variable named 'ID' and covariates named 'covariate_x_1','covariate_x_2',... |
gridy |
The Y-grid sf dataframe, containing a numeric area ID variable named 'ID', covariates named 'covariate_y_1','covariate_y_2',...., and an outcome named 'y'. |
atoms |
The atom sf dataframe, which should contain numeric variables named 'ID_x' and 'ID_y' holding the X-grid and Y-grid cell IDs for each atom, as well as an atom-level population count named 'population'. |
model_code |
NIMBLE model code from get_abrm_model() |
true_params |
The true outcome model regression coefficient parameters, if known (e.g., from simulate_misaligned_data()) |
norm_idx_x |
Vector of numeric indices of X-grid covariates (ordered as 'covariate_x_1','covariate_x_2',...) that should be treated as normally-distributed |
pois_idx_x |
Vector of numeric indices of X-grid covariates (ordered as 'covariate_x_1','covariate_x_2',...) that should be treated as Poisson-distributed |
binom_idx_x |
Vector of numeric indices of X-grid covariates (ordered as 'covariate_x_1','covariate_x_2',...) that should be treated as binomial-distributed |
norm_idx_y |
Vector of numeric indices of Y-grid covariates (ordered as 'covariate_y_1','covariate_y_2',...) that should be treated as normally-distributed |
pois_idx_y |
Vector of numeric indices of Y-grid covariates (ordered as 'covariate_y_1','covariate_y_2',...) that should be treated as Poisson-distributed |
binom_idx_y |
Vector of numeric indices of Y-grid covariates (ordered as 'covariate_y_1','covariate_y_2',...) that should be treated as binomial-distributed |
dist_y |
Distribution type for outcome (1=normal, 2=poisson, 3=binomial) |
niter |
Number of MCMC iterations (default: 50000) |
nburnin |
Number of burn-in iterations (default: 30000) |
nchains |
Number of MCMC chains (default: 2) |
thin |
Thinning interval (default: 10) |
seed |
Integer seed for reproducibility. Each chain uses seed+(chain_number-1) (default: NULL) |
sim_metadata |
Optional named list of simulation metadata (e.g.,
|
save_plots |
Logical, whether to save diagnostic plots (default: FALSE) |
output_dir |
Directory for saving outputs (default: NULL) |
compute_waic |
Logical; if |
An object of class "abrm": a named list with components
mcmc_results, parameter_estimates, all_parameters,
fitted_values (numeric vector of Y-grid-level predicted outcome
values on the original outcome scale),
y_grid_ids (integer vector of Y-grid cell IDs in model order), and
y_observed (numeric vector of observed outcome values for directly
observed Y-grid cells). Use fitted() to extract a formatted
comparison table of observed vs. fitted values, and waic() to
retrieve model fit criteria when compute_waic = TRUE.
# Simulate misaligned spatial data with one normal covariate per grid sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) # Retrieve the pre-compiled NIMBLE model code model_code <- get_abrm_model() # Fit the ABRM (use small niter/nburnin for illustration only) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) print(results) # concise model summary summary(results) # full parameter table plot(results) # MCMC trace and density plots vcov(results) # posterior variance-covariance matrices# Simulate misaligned spatial data with one normal covariate per grid sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) # Retrieve the pre-compiled NIMBLE model code model_code <- get_abrm_model() # Fit the ABRM (use small niter/nburnin for illustration only) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) print(results) # concise model summary summary(results) # full parameter table plot(results) # MCMC trace and density plots vcov(results) # posterior variance-covariance matrices
Run NIMBLE Model with Diagnostics
run_nimble_model( constants, data, inits, sim_metadata = NULL, model_code, niter = 50000, nburnin = 30000, nchains = 2, thin = 10, seed = NULL, save_plots = FALSE, output_dir = NULL, compute_waic = FALSE )run_nimble_model( constants, data, inits, sim_metadata = NULL, model_code, niter = 50000, nburnin = 30000, nchains = 2, thin = 10, seed = NULL, save_plots = FALSE, output_dir = NULL, compute_waic = FALSE )
constants |
List of model constants |
data |
List of data |
inits |
List of initial values |
sim_metadata |
List with simulation metadata (optional) |
model_code |
NIMBLE code object |
niter |
Number of MCMC iterations (default: 50000) |
nburnin |
Number of burn-in iterations (default: 30000) |
nchains |
Number of MCMC chains (default: 2) |
thin |
Thinning interval (default: 10) |
seed |
Integer seed for reproducibility. Each chain uses seed+(chain_number-1) (default: NULL) |
save_plots |
Logical, whether to save diagnostic plots (default: FALSE) |
output_dir |
Directory for saving plots (default: NULL) |
compute_waic |
Logical; if |
A named list with components samples (per-chain MCMC sample matrices), summary (posterior summary statistics), and
convergence (Gelman-Rubin statistics, effective sample sizes, and optional diagnostic plots).
# The recommended workflow is to call run_abrm(), which calls # run_nimble_model() internally after preparing all inputs. sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) model_code <- get_abrm_model() results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) summary(results)# The recommended workflow is to call run_abrm(), which calls # run_nimble_model() internally after preparing all inputs. sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) model_code <- get_abrm_model() results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) summary(results)
Simulate Misaligned Spatial Data
simulate_misaligned_data( seed = 2, dist_covariates_x = c("normal", "poisson", "binomial"), dist_covariates_y = c("normal", "poisson", "binomial"), dist_y = "poisson", x_intercepts = rep(0, 3), y_intercepts = rep(0, 3), rho_x = 0.6, rho_y = 0.6, x_correlation = 0.5, y_correlation = 0.5, beta0_y = NULL, beta_x = NULL, beta_y = NULL, diff_pops = TRUE, xy_cov_cor = FALSE )simulate_misaligned_data( seed = 2, dist_covariates_x = c("normal", "poisson", "binomial"), dist_covariates_y = c("normal", "poisson", "binomial"), dist_y = "poisson", x_intercepts = rep(0, 3), y_intercepts = rep(0, 3), rho_x = 0.6, rho_y = 0.6, x_correlation = 0.5, y_correlation = 0.5, beta0_y = NULL, beta_x = NULL, beta_y = NULL, diff_pops = TRUE, xy_cov_cor = FALSE )
seed |
Random seed (default = 2) |
dist_covariates_x |
Vector specifying distribution type for each synthetic X-grid covariate ('poisson', 'binomial', or 'normal') |
dist_covariates_y |
Vector specifying distribution type for each synthetic Y-grid covariate ('poisson', 'binomial', or 'normal') |
dist_y |
Distribution type for synthetic outcome variable (one of 'poisson', 'binomial', or 'normal') |
x_intercepts |
Intercepts for X covariates |
y_intercepts |
Intercepts for Y covariates |
rho_x |
Spatial correlation parameter for X-grid covariates (0 to 1 with higher values yielding more spatial correlation, default = 0.6) |
rho_y |
Spatial correlation parameter for Y-grid covariates and outcome (0 to 1 with higher values yielding more spatial correlation, default = 0.6) |
x_correlation |
Between-variable correlation for all pairs of X-grid covariates (default = 0.5) |
y_correlation |
Between-variable correlation for all pairs of Y-grid covariates (default = 0.5) |
beta0_y |
Intercept for outcome model |
beta_x |
Outcome model coefficients for X-grid covariates |
beta_y |
Outcome model coefficients for Y-grid covariates |
diff_pops |
Logical, indicating whether the atoms should be generated with different population sizes (diff_pops = TRUE) or a common population size (diff_pops = FALSE) |
xy_cov_cor |
Logical, indicating whether the atom-level spatial random effects for X-grid and Y-grid covariates should be correlated (xy_cov_cor = TRUE) or not. When set to TRUE, the x_correlation and rho_x parameters are used to generate all covariates (separate correlation parameters are not allowed for X-grid and Y-grid covariates). |
An object of class "misaligned_data": a named list with components gridy (Y-grid sf data frame with outcome and Y
covariates), gridx (X-grid sf data frame with X covariates), atoms (atom-level sf data frame, the spatial intersection of the
two grids), and true_params (list of the true regression coefficients used for data generation).
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c('normal', 'poisson', 'binomial'), dist_covariates_y = c('normal', 'poisson', 'binomial'), dist_y = "poisson", x_intercepts = c(4, -1, -1), y_intercepts = c(4, -1, -1), beta0_y = -1, x_correlation = 0.5, y_correlation = 0.5, beta_x = c(-0.03, 0.1, -0.2), beta_y = c(0.03, -0.1, 0.2) ) class(sim_data) # "misaligned_data" print(sim_data) # grid dimensions and atom count summary(sim_data) # includes true beta values names(sim_data$atoms) # atom-level variables sim_data$true_params$beta_x # true X-grid coefficientssim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c('normal', 'poisson', 'binomial'), dist_covariates_y = c('normal', 'poisson', 'binomial'), dist_y = "poisson", x_intercepts = c(4, -1, -1), y_intercepts = c(4, -1, -1), beta0_y = -1, x_correlation = 0.5, y_correlation = 0.5, beta_x = c(-0.03, 0.1, -0.2), beta_y = c(0.03, -0.1, 0.2) ) class(sim_data) # "misaligned_data" print(sim_data) # grid dimensions and atom count summary(sim_data) # includes true beta values names(sim_data$atoms) # atom-level variables sim_data$true_params$beta_x # true X-grid coefficients
Returns a structured summary of a fitted ABRM, including posterior means,
standard deviations, and 95
coefficients. When true parameter values are available (i.e., the model
was fitted on simulated data via simulate_misaligned_data),
mean absolute bias and credible-interval coverage are also reported.
## S3 method for class 'abrm' summary(object, ...)## S3 method for class 'abrm' summary(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
An object of class "summary.abrm" (a list) with components:
n_paramsNumber of regression parameters estimated.
mean_abs_biasMean absolute bias across parameters, or
NA if true values were not supplied.
coveragePercentage of parameters whose true value falls
within its 95% credible interval, or NA if unavailable.
estimatesData frame of posterior summaries.
The object is printed via print.summary.abrm.
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1 ) summary(results) # full parameter table with bias and coveragesim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1 ) summary(results) # full parameter table with bias and coverage
Prints grid dimensions, atom counts, and the true regression coefficients
(beta_x and beta_y) used to generate the data.
## S3 method for class 'misaligned_data' summary(object, ...)## S3 method for class 'misaligned_data' summary(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns object. Called for its side effect of
printing the data summary including the true parameter values.
## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) summary(sim_data)## Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) summary(sim_data)
Extracts variance-covariance matrices for regression coefficients from MCMC posterior samples. Returns separate matrices for X-grid and Y-grid coefficients.
## S3 method for class 'abrm' vcov(object, ...)## S3 method for class 'abrm' vcov(object, ...)
object |
An object of class |
... |
Additional arguments (unused) |
The variance-covariance matrices are computed from the posterior samples of the MCMC chains. If multiple chains were run, samples are combined across chains before computing covariances.
A list with class "vcov.abrm" containing:
vcov_beta_x |
Variance-covariance matrix for X-grid coefficients |
vcov_beta_y |
Variance-covariance matrix for Y-grid coefficients |
vcov_beta_0 |
Variance of the intercept (scalar) |
vcov_all |
Combined variance-covariance matrix for all parameters |
## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Extract posterior variance-covariance matrices vcov_mats <- vcov(results) # Posterior standard errors for X-grid and Y-grid coefficients sqrt(diag(vcov_mats$vcov_beta_x)) sqrt(diag(vcov_mats$vcov_beta_y)) # Posterior correlation matrix across all beta parameters cov2cor(vcov_mats$vcov_all)## Step 1: Simulate misaligned spatial data sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = c("normal", "poisson"), dist_covariates_y = c("normal", "poisson"), dist_y = "poisson", x_intercepts = c(0, -1), y_intercepts = c(0, -1), beta0_y = -1, beta_x = c(0.1, -0.05), beta_y = c(-0.1, 0.05) ) ## Step 2: Retrieve the NIMBLE model code model_code <- get_abrm_model() ## Step 3: Fit the ABRM results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = model_code, true_params = sim_data$true_params, norm_idx_x = 1, pois_idx_x = 2, norm_idx_y = 1, pois_idx_y = 2, dist_y = 2, niter = 1000, nburnin = 500, nchains = 2, seed = 1, save_plots = FALSE ) ## Step 4: Extract posterior variance-covariance matrices vcov_mats <- vcov(results) # Posterior standard errors for X-grid and Y-grid coefficients sqrt(diag(vcov_mats$vcov_beta_x)) sqrt(diag(vcov_mats$vcov_beta_y)) # Posterior correlation matrix across all beta parameters cov2cor(vcov_mats$vcov_all)
Extracts the Widely Applicable Information Criterion (WAIC) from a fitted model object.
waic(object, ...)waic(object, ...)
object |
A fitted model object. |
... |
Additional arguments passed to methods. |
A "waic_abrm" object for abrm models. See
waic.abrm for details.
Returns the Widely Applicable Information Criterion (WAIC) for a fitted
ABRM. WAIC is the recommended information criterion for Bayesian models and
is computed by NIMBLE when compute_waic = TRUE is passed to
run_abrm.
## S3 method for class 'abrm' waic(object, ...)## S3 method for class 'abrm' waic(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
WAIC (Watanabe 2010) is defined as:
where is the log pointwise predictive density and
is the effective number of parameters. Unlike AIC,
WAIC is fully Bayesian and averages over the posterior distribution rather
than evaluating at a point estimate.
WAIC values are only meaningful when compared between models fitted on
identical data with the same outcome distribution. A true marginal
log-likelihood is not available for ABRM, so AIC and
logLik are not supported; use waic() instead.
An object of class "waic_abrm", which is a named list with
components:
waicThe scalar WAIC value. Lower values indicate better predictive accuracy penalised for model complexity.
lppdLog pointwise predictive density — the goodness-of-fit component of WAIC.
penaltyThe effective number of parameters
() — the complexity penalty component.
n_paramsThe number of regression parameters estimated.
Returns an error if the model was not fitted with
compute_waic = TRUE.
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571–3594.
waic, print.waic_abrm,
run_abrm
sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, compute_waic = TRUE # required; re-fit if this was FALSE ) w <- waic(results) print(w) # formatted table: WAIC, lppd, pWAIC, n_params w$waic # raw scalar for comparing models w$penalty # effective number of parameters (pWAIC)sim_data <- simulate_misaligned_data( seed = 1, dist_covariates_x = "normal", dist_covariates_y = "normal", dist_y = "normal", x_intercepts = 0, y_intercepts = 0, beta0_y = 0, beta_x = 0.1, beta_y = -0.1 ) results <- run_abrm( gridx = sim_data$gridx, gridy = sim_data$gridy, atoms = sim_data$atoms, model_code = get_abrm_model(), norm_idx_x = 1, norm_idx_y = 1, dist_y = 1, niter = 1000, nburnin = 500, nchains = 2, seed = 1, compute_waic = TRUE # required; re-fit if this was FALSE ) w <- waic(results) print(w) # formatted table: WAIC, lppd, pWAIC, n_params w$waic # raw scalar for comparing models w$penalty # effective number of parameters (pWAIC)