BayesMRA (version 1.0.0)

mcmc_mra: 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(
  y,
  X,
  locs,
  params,
  priors = NULL,
  M = 4,
  n_neighbors = 68,
  n_coarse_grid = 10,
  n_padding = 5L,
  n_cores = 1L,
  inits = NULL,
  config = NULL,
  verbose = FALSE,
  use_spam = TRUE,
  n_chain = 1
)

Arguments

y

is a \(n\) vector of Gaussian data.

X

is a \(n \times p\) matrix of fixed effects (like latitude, elevation, etc)

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.

priors

is the list of prior settings.

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).

n_cores

is the number of cores for parallel computation using openMP.

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.

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.

verbose

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

use_spam

is a boolean flag to determine whether the output is a list of spam matrix objects (use_spam = TRUE) or a an \(n \times n\) sparse Matrix of class "dgCMatrix" use_spam = FALSE (see spam and Matrix packages for details).

n_chain

is the MCMC chain id. The default is 1.

Examples

# NOT RUN {
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 <- do.call(cbind, MRA$W)

n_dims   <- rep(NA, length(MRA$W))
dims_idx <- NULL
for (i in 1:M) {
    n_dims[i] <- ncol(MRA$W[[i]])
    dims_idx  <- c(dims_idx, rep(i, n_dims[i]))
}
## 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
A_constraint <- sapply(1:M, function(i){
    tmp = rep(0, sum(n_dims))
    tmp[dims_idx == i] <- 1
    return(tmp)
})
a_constraint <- rep(0, M)
alpha <- as.vector(spam::rmvnorm.prec.const(
    n = 1,
    mu = rep(0, nrow(W)),
    Q = Q_alpha_tau2,
    A = t(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(
    y             = y,
    X             = X,
    locs          = locs,
    params        = params,
    priors        = priors,
    M             = 2,
    n_coarse_grid = 4,
    n_cores       = 1L,
    verbose       = FALSE
)

# }