mcmc_mra_vecchia: Bayesian Multi-resolution Spatial Regression

View source: R/mcmc_mra_vecchia.R

mcmc_mra_vecchiaR Documentation

Bayesian Multi-resolution Spatial Regression

Description

this function runs Markov Chain Monte Carlo to estimate the Bayesian multi-resolution spatial regression model.

Usage

mcmc_mra_vecchia(
  y,
  locs,
  params,
  M = 3,
  n_neighbors = 68,
  n_coarse_grid = 10,
  n_padding = 5L,
  inits = NULL,
  normalize = TRUE,
  verbose = FALSE,
  config = NULL,
  constraint = "overall",
  n_chain = 1
)

Arguments

y

is a n vector of Gaussian data.

locs

is a n \times 2 matrix of observation locations.

params

is a list of parameter settings. The list params must contain the following values:

  • n_adapt: A positive integer number of adaptive MCMC iterations.

  • n_mcmc: A positive integer number of total MCMC iterations post adaptation.

  • n_thin: A positive integer number of MCMC iterations per saved sample.

  • n_message: A positive integer number of frequency of iterations to output a progress message. For example, n_message = 50 outputs progress messages every 50 iterations.

M

The number of resolutions.

n_neighbors

The expected number of neighbors for each interior basis function. This determines the basis radius parameter.

n_coarse_grid

The number of basis functions in one direction (e.g. n_coarse_grid = 10 results in a 10 \times 10 course grid which is further extended by the number of additional padding basis functions given by n_padding.

n_padding

The number of additional boundary points to add on each boundary. For example, n_padding = 5 will add 5 boundary knots to the both the left and right side of the grid).

inits

is the list of initial values if the user wishes to specify initial values. If these values are not specified, then the initial values will be randomly sampled from the prior.

normalize

A logical value of whether to normalize the basis functions or not.

verbose

Should verbose output be printed? Typically this is only useful for troubleshooting.

config

is the list of configuration values if the user wishes to specify initial values. If these values are not specified, then default a configuration will be used.

constraint

What constraint should be applied to the spatial process? Options include no constraint (constraint = "unconstrained"), a constraint so the entire MRA process sums to 0 (constraint = "overall"), a constraint so that each of the M levels of the MRA process sum to 0 (constraint = "resolution"), or whether the predicted process must sum to 0 (constraint = "predicted"). Note: constraint = "predicted" is NOT currently supported.

n_chain

is the MCMC chain id. The default is 1.

Examples

set.seed(111)
## genereate the locations
locs <- matrix(runif(20), 10, 2)
## generate some covariates and regression coefficients
X <- cbind(1, matrix(rnorm(30), 10, 3))
beta <- rnorm(ncol(X))

## simulate the MRA process
M <- 2
MRA <- mra_wendland_2d(locs, M = M, n_coarse_grid = 4)
W <- MRA$W
n_dims <- MRA$n_dims
dims_idx <- MRA$dims_idx

## set up the process precision matrices
Q_alpha <- make_Q_alpha_2d(sqrt(n_dims), c(0.9, 0.8))
Q_alpha_tau2 <- make_Q_alpha_tau2(Q_alpha, tau2 = c(2, 4))

## add in constraints so each resolution has random effects that sum to 0
constraints <- make_constraint(MRA, constraint = "resolution", joint = TRUE)
A_constraint <- constraints$A_constraint
a_constraint <- constraints$a_constraint
alpha <- as.vector(spam::rmvnorm.prec.const(
    n = 1,
    mu = rep(0, nrow(W)),
    Q = Q_alpha_tau2,
    A = A_constraint,
    a = a_constraint))
## define the data
y <-  as.vector(X %*% beta + W %*% alpha + rnorm(10))

## define the params for MCMC fitting
params <- list(
    n_mcmc    = 5,
    n_adapt   = 5,
    n_thin    = 1,
    n_message = 5)

## define the model priors
priors <- list(
    alpha_tau2   = 1,
    beta_tau2    = 1,
    alpha_sigma2 = 1,
    beta_sigma2  = 1,
    mu_beta      = rep(0, ncol(X)),
    Sigma_beta   = 5 * diag(ncol(X)))

## Fit the MRA model using MCMC
# out     <- mcmc_mra_vecchia(
#     y             = y,
#     X             = X,
#     locs          = locs,
#     params        = params,
#     priors        = priors,
#     M             = 2,
#     n_coarse_grid = 4
# )


jtipton25/BayesMRA documentation built on Feb. 28, 2024, 1:27 p.m.