ggt.sp: General Gaussian Template for common spatial covariance functions

Description

The function ggt.sp performs Bayesian analysis of Gaussian models with spatially dependent error structure.

Usage

ggt.sp(formula, data = parent.frame(), coords, run.control, beta.control, cov.model, var.update.control,
                     DIC=TRUE, DIC.start = 1, verbose=TRUE, ...)

Arguments

formula

a symbolic description of the model to be fit. The details of model specification are given below.

data

an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which ggt.sp is called.

coords

an $n \times 2$ matrix of the observation coordinates in $R^2$ (e.g., easting and northing).

run.control

a list with tag n.samples the value of which is the number of MCMC iterations.

beta.control

a list with tags beta.update.method, tuning.matrix, beta.starting, beta.prior, and sample.order. The value of beta.update.method is either "gibbs" or

cov.model

a quoted key word that specifies the covariance function used to model the spatial dependence structure among the observations. Supported covariance model key words are: "exponential", "matern", "spherical"

var.update.control

a list with each tag corresponding to a parameter's name in the specified cov.model. The value portion of each of these tags is another list with tags sample.order, starting, tuning, and

DIC

if TRUE, DIC and associated statistics will be calculated to assess model fit.

DIC.start

the sampler iteration to start the DIC calculation. This is useful for those who choose to acknowledge chain burn-in.

verbose

if TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.

...

currently no additional arguments.

Value

An object of class ggt.sp, which is a list with the following tags:
cov.modelthe covariance model specified by cov.model.
no.tausqa logical value indicating if the parameter $\tau^2$ was used.
coordsthe $n \times 2$ matrix specified by coords.
Xthe $n \times p$ matrix of regressors specified by formula.
(Y){the the $n \times 1$ response variable vector specified by formula. } p.samples{a matrix of posterior sample for the regressors and parameters used in the specified cov.model. } acceptance{a matrix of the Metropolis-Hastings sampling acceptance rate. Row names correspond to the parameters used in the var.function and if the beta.update.method is set to "mh" the acceptance rate for Metropolis-Hastings block proposal of the regressors is included. } DIC{a matrix that holds $DIC$ and associated values $\bar{D}$, $D[\bar{\theta}]$, and $pD$ (Banerjee et al. 2004 and Spiegelhalter et al. 2002). }
The Gaussian model assumes that $y \sim N(X \beta, \Sigma_y)$, where $y$ is the $n \times 1$ response vector, $X$ is the $n \times p$ matrix of regressors, and $\Sigma_y = (\sigma^2 R(\phi)+\tau^2 I)$, where $R$ is a valid correlation function specified by cov.model. Under this marginalized model, we assume $\Sigma_y \sim MVN(0, \sigma^2 R(\phi)+\tau^2 I)$. The unmarginalized model partitions the variance terms into spatial and non-spatial components which follow $MVN(0, \sigma^2 R(\phi))$ and $MVN(0, \tau^2 I)$, respectively. Once the marginalized model is fit with ggt.sp the vector of random spatial effects can be recovered with sp.effect function (see Banerjee et al. 2004, Chapter 5 for further details).
The parameters associated with the $X$ are regressors and defined in the beta.control list. The cov.model defines the structure of $n \times n$ $\Sigma_y$. Covariance models and associated parameters are:
- Exponential:
{$\sigma^2 exp(-\phi d)+\tau^2$, for $d > 0$, where $sigma^2$ is the partial sill, $\tau^2$ is the nugget, $\phi$ is the decay parameter, and $d$ is the distance between any two sites. } Gaussian:{$\sigma^2 exp(-\phi^2 d^2)+\tau^2$, for $d > 0$, where $sigma^2$ is the partial sill, $\tau^2$ is the nugget, $\phi$ is the decay parameter, and $d$ is the distance between any two sites. } Spherical:{$\sigma^2 (1 - 3/2 \phi d + 1/2(\phi d)^3)+\tau^2$, for $0 < d <= 1="" \phi$,="" where="" $sigma^2$="" is="" the="" partial="" sill,="" $\tau^2$="" nugget,="" $\phi$="" decay="" parameter,="" and="" $d$="" distance="" between="" any="" two="" sites.="" }="" Matern:{$\sigma^2/(2^(\nu - 1) \Gamma(\nu)) (2 \nu^0.5 d \phi)^\nu K_\nu (2 \nu^0.5 d \phi)$, for $d > 0$, where $sigma^2$ is the partial sill, $\tau^2$ is the nugget, $\phi$ is the decay parameter, $\nu$ is decay smoothness, $\Gamma$ is the gamma function, and $K_\nu$ is the modified Bessel of the second kind. }
For all covariance models the $\tau^2$ parameter is optional. To omit $\tau^2$ from the model, do not include it in the var.update.control The update method for the regressors is either the Metropolis-Hastings algorithm or Gibbs. If the Metropolis-Hastings algorithm is chosen, then updates come from a single block draw from a multi-variate normal distribution with step size determined by the tuning matrix. Otherwise, a Gibbs sampler is used to sample from the conditional posterior distribution of the regressors given the error terms.
The Metropolis-Hastings algorithm is used to estimate the joint posterior distribution of the model's variance and, if specified, regression parameters. In the simplest case, for each iteration of the algorithm a single a single vector of candidate values is drawn from a multivariate normal distribution with variance matrix specified by the tuning values in the var.update.control and beta.control lists. This single block draw is specified by setting the sample.order values equal in the var.update.control and beta.control lists. In so doing, a single acceptance rate is monitored for all parameters. An acceptance rate of about 30 percent is generally recommended (see Gelman et al., 1996). Hand tuning can be tricky, and it is often useful to update parameters independently or in smaller sets (i.e., blocks). This is accomplished by setting the sample.order accordingly. For instance, to update variance parameters as a single block separate from the regressors, set the all sample.order values in the var.update.control list equal and different than the sample.order in the beta.control list. In this way, two separate acceptance rates are monitored and reported. The extreme is setting all sample.order values different, which allows each parameter to be monitored independently. However, the likelihood is evaluated for each each parameter block, this requires an inversion of the covariance matrix, which if $n$ is large could take a significant amount of time.

References

Banerjee, S., Carlin, B.P., and Gelfand, A.E. (2004). Hierarchical modeling and analysis for spatial data. Chapman and Hall/CRC Press, Boca Raton, Fla.

Gelman, A., Roberts, G.O., and Gilks, W.R. (1996). Efficient Metropolis jumping rules. In Bayesian Statistics 5, eds. J.M. Bernardo, J.O. Berger, A.P. Dawid, and A.F.M. Smith. Oxford: Oxfod University Press, pp. 599-607. Spiegelhalter, D.J., Best, N., Carlin, B.P., and van der Linde, A. (2002). Bayesian measures of model complexity and fit (with discussion). J. Roy. Statist. Soc., Ser. B, 64, 583-639.

Further information on the package spBayes can be found at: http://blue.fr.umn.edu/spatialBayes.

Examples

Run this code

###############################################
##Subset data
###############################################
set.seed(1234)

data(BEF)

n.subset <- 100
BEF.subset <- BEF[sample(1:nrow(BEF),n.subset),]

##############################################
##General ggt.sp setup and model fit
##############################################

##Specify the priors, hyperparameters, and variance parameter starting values.
sigmasq.prior <- prior(dist="IG", shape=2, scale=30)
tausq.prior <- prior(dist="IG", shape=2, scale=30)
##The prior on phi corresponds to a prior of 500-2000 meters
##for the effective range (i.e., -log(0.05)/0.0015, -log(0.05)/0.006, when
##using the exponential covariance function).
phi.prior <- prior(dist="LOGUNIF", a=0.0015, b=0.006)

var.update.control <- list("sigmasq"=list(sample.order=0, starting=30, tuning=0.5, prior=sigmasq.prior),
                           "tausq"=list(sample.order=0, starting=30, tuning=0.5, prior=tausq.prior),
                           "phi"=list(sample.order=0, starting=0.006, tuning=0.8, prior=phi.prior))

##Specify the number of samples.
run.control <- list("n.samples" = 1500)

##Specify the model, assign the prior, and starting values for the regressors.
model <- BE_BA_AREA~ELEV+SPR_02_TC1+SPR_02_TC2+SPT_02_TC3

##Note, precision of 0 is same as flat prior alternatively if you want a flat prior
##do not include the beta.prior tag in the beta.control list.
beta.prior <- prior(dist="NORMAL", mu=rep(0, 5), precision=diag(0, 5))

beta.control <- list(beta.update.method="gibbs", beta.starting=rep(0, 5), beta.prior=beta.prior)

model.coords <- cbind(BEF.subset$XUTM, BEF.subset$YUTM)

ggt.sp.out <-ggt.sp(formula=model, data=BEF.subset,
                    coords=model.coords,
                    var.update.control=var.update.control,
                    beta.control=beta.control,
                    cov.model="exponential",
                    run.control=run.control, DIC=TRUE, verbose=TRUE)

print(ggt.sp.out$accept)
print(ggt.sp.out$DIC)

##Get the effective range of spatial dependence.
eff.range <- -log(0.05)/ggt.sp.out$p.samples[,"phi"]

##Plot the chain with coda.
mcmc.obj <- mcmc(cbind(ggt.sp.out$p.samples, eff.range))
plot(mcmc.obj)

##############################################
##ggt.sp setup and model fit with a
##fixed variance parameter
##############################################

var.update.control <- list("sigmasq"=list(sample.order=0, starting=30, tuning=0.5, prior=sigmasq.prior),
                           "tausq"=list(sample.order=0, starting=30, tuning=0.5, prior=tausq.prior),
                           "phi"=list(fixed=0.001))

ggt.sp.out <-ggt.sp(formula=model, data=BEF.subset,
                    coords=model.coords,
                    var.update.control=var.update.control,
                    beta.control=beta.control,
                    cov.model="exponential",
                    run.control=run.control, DIC=TRUE, verbose=TRUE)

print(ggt.sp.out$accept)
print(ggt.sp.out$DIC)

##Summarize and plot chain with coda
mcmc.obj <- mcmc(ggt.sp.out$p.samples)
summary(mcmc.obj)
##plot(mcmc.obj)



##############################################
##ggt.sp setup and model fit with single
##block update of the variance and regression
##parameters
##############################################

var.update.control <- list("sigmasq"=list(sample.order=0, starting=30, tuning=0.5, prior=sigmasq.prior),
                           "tausq"=list(sample.order=0, starting=30, tuning=0.5, prior=tausq.prior),
                           "phi"=list(sample.order=0, starting=0.006, tuning=0.8, prior=phi.prior))

#Using Metropolis-Hastings and default flat prior for regression parameters.
beta.control <- list(beta.update.method="mh", sample.order=0)

ggt.sp.out <-ggt.sp(formula=model, data=BEF.subset,
                    coords=model.coords,
                    var.update.control=var.update.control,
                    beta.control=beta.control,
                    cov.model="exponential",
                    run.control=run.control, DIC=TRUE, verbose=TRUE)

##Note single acceptance rate.
print(ggt.sp.out$accept)

Run the code above in your browser using DataLab