linear.model.Bayes: Bayesian estimation for the geostatistical linear Gaussian model

Description

This function performs Bayesian estimation for the geostatistical linear Gaussian model.

Usage

linear.model.Bayes(
  formula,
  coords,
  data,
  kappa,
  control.mcmc,
  control.prior,
  low.rank = FALSE,
  knots = NULL,
  messages = TRUE
)

Arguments

formula

an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted.

coords

an object of class formula indicating the geographic coordinates.

data

a data frame containing the variables in the model.

kappa

shape parameter of the Matern covariance function.

control.mcmc

output from control.mcmc.Bayes.

control.prior

output from control.prior.

low.rank

logical; if low.rank=TRUE a low-rank approximation is fitted.

knots

if low.rank=TRUE, knots is a matrix of spatial knots used in the low-rank approximation. Default is knots=NULL.

messages

logical; if messages=TRUE then status messages are printed on the screen (or output device) while the function is running. Default is messages=TRUE.

Value

An object of class "Bayes.PrevMap". The function summary.Bayes.PrevMap is used to print a summary of the fitted model. The object is a list with the following components:

estimate: matrix of the posterior samples for each of the model parameters.

S: matrix of the posterior samplesfor each component of the random effect. This is only returned for the low-rank approximation.

y: response variable.

D: matrix of covariarates.

coords: matrix of the observed sampling locations.

kappa: vaues of the shape parameter of the Matern function.

knots: matrix of spatial knots used in the low-rank approximation.

const.sigma2: vector of the values of the multiplicative factor used to adjust the sigma2 in the low-rank approximation.

h1: vector of values taken by the tuning parameter h.theta1 at each iteration.

h2: vector of values taken by the tuning parameter h.theta2 at each iteration.

h3: vector of values taken by the tuning parameter h.theta3 at each iteration.

call: the matched call.

Details

This function performs Bayesian estimation for the geostatistical linear Gaussian model, specified as $$Y = d'\beta + S(x) + Z,$$ where $Y$ is the measured outcome, $d$ is a vector of coavariates, $\beta$ is a vector of regression coefficients, $S(x)$ is a stationary Gaussian spatial process and $Z$ are independent zero-mean Gaussian variables with variance tau2. More specifically, $S(x)$ has an isotropic Matern covariance function with variance sigma2, scale parameter phi and shape parameter kappa. The shape parameter kappa is treated as fixed.

Priors definition. Priors can be defined through the function control.prior. The hierarchical structure of the priors is the following. Let $\theta$ be the vector of the covariance parameters $(\sigma^2,\phi,\tau^2)$; then each component of $\theta$ can have independent priors freely defined by the user. However, uniform and log-normal priors are also available as default priors for each of the covariance parameters. To remove the nugget effect $Z$, no prior should be defined for tau2. Conditionally on sigma2, the vector of regression coefficients beta has a multivariate Gaussian prior with mean beta.mean and covariance matrix sigma2*beta.covar, while in the low-rank approximation the covariance matrix is simply beta.covar.

Updating the covariance parameters using a Metropolis-Hastings algorithm. In the MCMC algorithm implemented in linear.model.Bayes, the transformed parameters $$(\theta_{1}, \theta_{2}, \theta_{3})=(\log(\sigma^2)/2,\log(\sigma^2/\phi^{2 \kappa}), \log(\tau^2))$$ are independently updated using a Metropolis Hastings algorithm. At the $i$-th iteration, a new value is proposed for each from a univariate Gaussian distrubion with variance, say $h_{i}^2$, tuned according the following adaptive scheme $$h_{i} = h_{i-1}+c_{1}i^{-c_{2}}(\alpha_{i}-0.45),$$ where $\alpha_{i}$ is the acceptance rate at the $i$-th iteration (0.45 is the optimal acceptance rate for a univariate Gaussian distribution) whilst $c_{1} > 0$ and $0 < c_{2} < 1$ are pre-defined constants. The starting values $h_{1}$ for each of the parameters $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ can be set using the function control.mcmc.Bayes through the arguments h.theta1, h.theta2 and h.theta3. To define values for $c_{1}$ and $c_{2}$, see the documentation of control.mcmc.Bayes.

Low-rank approximation. In the case of very large spatial data-sets, a low-rank approximation of the Gaussian spatial process $S(x)$ might be computationally beneficial. Let $(x_{1},\dots,x_{m})$ and $(t_{1},\dots,t_{m})$ denote the set of sampling locations and a grid of spatial knots covering the area of interest, respectively. Then $S(x)$ is approximated as $\sum_{i=1}^m K(\|x-t_{i}\|; \phi, \kappa)U_{i}$, where $U_{i}$ are zero-mean mutually independent Gaussian variables with variance sigma2 and $K(.;\phi, \kappa)$ is the isotropic Matern kernel (see matern.kernel). Since the resulting approximation is no longer a stationary process (but only approximately), sigma2 may take very different values from the actual variance of the Gaussian process to approximate. The function adjust.sigma2 can then be used to (approximately) explore the range for sigma2. For example if the variance of the Gaussian process is 0.5, then an approximate value for sigma2 is 0.5/const.sigma2, where const.sigma2 is the value obtained with adjust.sigma2.

References

Diggle, P.J., Giorgi, E. (2019). Model-based Geostatistics for Global Public Health. CRC/Chapman & Hall.

Giorgi, E., Diggle, P.J. (2017). PrevMap: an R package for prevalence mapping. Journal of Statistical Software. 78(8), 1-29. doi: 10.18637/jss.v078.i08

Higdon, D. (1998). A process-convolution approach to modeling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics 5, 173-190.