ggt: General Gaussian Template

Description

The function ggt performs Bayesian analysis of Gaussian models with potentially complex hierarchical error structures.

Usage

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

Arguments

formula

a symbolic description of the regression portion of the model. 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 is called.

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

var.function

an R function that defines the error structure of the model. The output of this function must be a $n \times n$ positive definite matrix, where $n$ is the number of observations defined by the formula. The function must not have

var.update.control

a list with each tag corresponding to a parameter's name in the var.function. 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, which is a list with the following tags:
p.samplesa matrix of posterior sample for the regressors and parameters used in the var.function.
acceptancea 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.
var.functiona copy of the var.function.
DICa matrix that holds $DIC$ and associated values $\bar{D}$, $D[\bar{\theta}]$, and $pD$ (Banerjee et al. 2004 and Spiegelhalter et al. 2002).

Details

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 = Q_y\sigma^2$. The parameters associated with the $X$ are defined in the beta.control list. The var.function defines the structure of $n \times n$ $\Sigma_y$. For example, under the simple linear model, $Q_y$ is the $n \times n$ identity matrix and $\sigma^2$ is the only parameter defined in the var.update.control list. 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 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 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))

##Get the distance matrix.
D <- as.matrix(dist(cbind(BEF.subset$XUTM, BEF.subset$YUTM)))

##Define the variance function (e.g., exponential covariance)
var.function <- function(){
  sigmasq*exp(-phi*D)+tausq*diag(n.subset)
}

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


ggt.out <- ggt(formula=model, data=BEF.subset,
               var.update.control=var.update.control,
               beta.control=beta.control,
               var.function = var.function,
               run.control=run.control, DIC=TRUE, verbose=TRUE)

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

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

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

##############################################
##ggt 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.out <- ggt(formula=model, data=BEF.subset,
               var.update.control=var.update.control,
               beta.control=beta.control,
               var.function = var.function,
               run.control=run.control, DIC=TRUE, verbose=TRUE)

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

Run the code above in your browser using DataLab