spBayes (version 0.0-1)

sp.effect: Recovers the random spatial effects from a ggt.sp object

Description

The function sp.effect recovers the vector of random spatial effects from a given marginalized model contained in ggt.sp object.

Usage

sp.effect(ggt.sp.obj, start=1, end, thin=1, verbose=TRUE, ...)

Arguments

ggt.sp.obj

an object returned by ggt.sp (i.e., object of class ggt.sp).

start

specifies the first sample included in spatial process recovery. This is useful for those who choose to acknowledge chain burn-in.

end

specifies the last sample included in spatial process recovery. The default is to include from start to nrow(ggt.sp.obj$p.samples).

thin

a sample thinning factor. The default of 1 considers all samples between start and end. For example, if thin = 10 then 1 in 10 samples are considered between start and end.

verbose

if TRUE calculation progress is printed to the screen. Otherwise, nothing is printed to the screen.

...

currently no additional arguments.

Value

• coords
  the $n \times 2$ matrix of the observation coordinates from ggt.sp.obj$coords.
• pp.samples
  an $n \times m$ matrix that holds $m$ samples from the posterior predictive distribution of the spatial process. See below for details.

Details

The Gaussian model fit by ggt.sp 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.

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

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

##############################################
##recover random spatial effect
##Y = X^tB + w + e, where w is the random
##spatial effect and e is pure error
##############################################

##get the spatial random effect w
sp.effect.out <- sp.effect(ggt.sp.out, start=500, thin=2)
w <- rowMeans(sp.effect.out$pp.samples)

par(mfrow=c(1,2))
int.obs <- interp(BEF.subset$XUTM, BEF.subset$YUTM, BEF.subset$BE_BA_AREA)
image(int.obs, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Observed density
of American beech (Y)")
contour(int.obs, add=TRUE)

int.w <- interp(BEF.subset$XUTM, BEF.subset$YUTM, w)
image(int.w, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Spatial random effect density (w)")
contour(int.w, add=TRUE)