ref.MCMC {ref.ICAR}

MCMC for Reference Prior on an Intrinsic Conditional Autoregressive Random Effects Model for Areal Data

Description

Implements the Metropolis-within-Gibbs sampling algorithm proposed by Ferreira et al. (2021), to perform posterior inference for the intrinsic conditional autoregressive model with spatial random effects. This algorithm uses the spectral domain for the hierarchical model to create the Spectral Gibbs Sampler (SGS), which provides notable speedups to the MCMC algorithm proposed by Keefe et al (2019).

Usage

ref.MCMC(
  y,
  X,
  H,
  iters = 10000,
  burnin = 5000,
  verbose = TRUE,
  tauc.start = 1,
  beta.start = 1,
  sigma2.start = 1,
  step.tauc = 0.5,
  step.sigma2 = 0.5
)

Arguments

y

Vector of responses.

X

Matrix of covariates. This should include a column of 1's for models with a non-zero intercept.

H

The neighborhood matrix.

iters

Number of MCMC iterations to perform. Defaults to 10,000.

burnin

Number of MCMC iterations to discard as burn-in. Defaults to 5,000.

verbose

If FALSE, MCMC progress is not printed.

tauc.start

Starting value for the spatial dependence parameter.

beta.start

Starting value for the vector of fixed effect regression coefficients.

sigma2.start

Starting value for the variance of the unstructured random effects.

step.tauc

Step size for the spatial dependence parameter

step.sigma2

Step size for the variance of the unstructured random effects.

Value

A list containing MCMC chains and parameter summaries.

MCMCchain

Matrix of MCMC chains.

tauc.MCMC

MCMC chains for the spatial dependence parameter.

sigma2.MCMC

MCMC chains for the variance of the unstructured random effects.

phi.MCMC

MCMC chains for the spatial random effects.

beta.MCMC

MCMC chains for the fixed effect regression coefficients.

accept.sigma2

Final acceptance number for variance of the unstructured random effects.

accept.tauc

Final acceptance number for spatial dependence parameter.

accept.phi

Final acceptance number for spatial random effects.

Author(s)

Erica M. Porter, Matthew J. Keefe, Christopher T. Franck, and Marco A.R. Ferreira

References

Keefe MJ, Ferreira MAR, Franck CT (2019). “Objective Bayesian analysis for Gaussian hierarchical models with intrinsic conditional autoregressive priors.” Bayesian Analysis, 14(1), 181 – 209. doi:10.1214/18-BA1107.

Keefe MJ, Ferreira MAR, Franck CT (2018). “On the formal specification of sum-zero constrained intrinsic conditional autoregressive models.” Spatial Statistics, 24, 54–65. doi:10.1016/j.spasta.2018.03.007.

Ferreira MAR, Porter EM, Franck CT (2021). “Fast and scalable computations for Gaussian hierarchical models with intrinsic conditional autoregressive spatial random effects.” Computational Statistics and Data Analysis, 162, 107264. ISSN 0167-9473, doi:10.1016/j.csda.2021.107264, https://www.sciencedirect.com/science/article/pii/S0167947321000980.

Examples

#### Fit the model for simulated areal data on a grid ####

### Load extra libraries
library(sp)
library(methods)
library(spdep)
library(mvtnorm)

### Generate areal data on a grid
rows=5; cols=5
tauc=1
sigma2=2; beta=c(1,5)

### Create grid
grid <- GridTopology(c(1,1), c(1,1), c(cols,rows))
polys <- as(grid, "SpatialPolygons")
spgrid <- SpatialPolygonsDataFrame(polys,data=data.frame(row.names=row.names(polys)))

### Create neighborhood matrix
grid.nb <- poly2nb(spgrid,queen=FALSE)
W <- nb2mat(grid.nb, style="B")

### Put spatially correlated data in grid
p <- length(beta)
num.reg <- (rows*cols)
if(p>1){x1<-rmvnorm(n=num.reg,mean=rep(0,p-1),sigma=diag(p-1))} else{x1<-NULL}
X <- cbind(rep(1,num.reg),x1)
Dmat <- diag(apply(W,1,sum))
H <- Dmat - W
row.names(H) <- NULL

### Obtain true response vector
theta_true <- rnorm(num.reg,mean=0,sd=sqrt(sigma2))
Q <- eigen(H,symmetric=TRUE)$vectors
eigH <- eigen(H,symmetric=TRUE)$values
D <- diag(eigH)
Qmat <- Q[,1:(num.reg-1)]
phimat <- diag(1/sqrt(eigH[1:(num.reg-1)]))
z <- t(rmvnorm(1,mean=rep(0,num.reg-1),sigma=diag(num.reg-1)))
phi_true <- sqrt((1/tauc)*sigma2)*(Qmat%*%phimat%*%z)
Y <- X%*%beta + theta_true + phi_true

### Fit the model
set.seed(5432)
model <- ref.MCMC(y=Y,X=X,H=H,iters=15000,burnin=5000,verbose=TRUE,tauc.start=.1,beta.start=-1,
sigma2.start=.1,step.tauc=0.5,
      step.sigma2=0.5)
      

#### Small example for checking
model <- ref.MCMC(y=Y,X=X,H=H,iters=1000,burnin=50,verbose=TRUE,tauc.start=.1,beta.start=-1,
sigma2.start=.1,step.tauc=0.5,
      step.sigma2=0.5)