sp.predict: Prediction for new points given a ggt.sp object

Description

The function sp.predict collects a posterior predictive sample for a set of new points given a ggt.sp object.

Usage

sp.predict(ggt.sp.obj, pred.joint = TRUE, pred.coords, pred.covars, use.covar.names = TRUE, start=1, end, thin=1, verbose=TRUE, ...)

Arguments

ggt.sp.obj

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

pred.joint

if TRUE the spatial dependence among the prediction points is preserved; otherwise, predicted points are considered independent. See below for details.

pred.coords

an $m \times 2$ matrix of $m$ prediction point coordinates in $R^2$ (e.g., easting and northing). The first column is assumed to be easting coordinates and the second column northing coordinates (i.e., x,y coordinates).

pred.covars

an $m \times q$ matrix or data frame containing the regressors associated with pred.coords. Several scenarios exist which determine the specification of this matrix or data frame. If the model contained in ggt.sp.obj

use.covar.names

if TRUE an attempt is made to match regressor names between ggt.sp.obj$X and pred.covars; otherwise, column order must match between ggt.sp.obj$X and pred.covars.

start

specifies the first sample included in the prediction calculation. This is useful for those who choose to acknowledge chain burn-in.

end

specifies the last sample included in the prediction calculation. 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

obs.coordsthe $n \times 2$ matrix of the observation coordinates from ggt.sp.obj$coords.
pred.coordsthe $m \times 2$ matrix of prediction point coordinates specified by pred.coords.
pp.samplesan $n \times m$ matrix that holds $m$ samples from the posterior predictive distribution of the response variable specified in the ggt.sp.obj.

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

###############################################
##prediction
###############################################

##resample
n.pred <- 100
BEF.pred <- BEF[sample(1:nrow(BEF),n.pred),]

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

pred <- sp.predict(ggt.sp.out, pred.joint=TRUE, pred.coords=pred.coords, pred.covars=BEF.pred, start=500, thin=2)

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")
contour(int.obs, add=TRUE)

int.pred <- interp(BEF.pred$XUTM, BEF.pred$YUTM, rowMeans(pred$pp.samples))
image(int.pred, xlab="Longitude (meters)", ylab="Latitude (meters)", main="Predicted density
of American beech")
contour(int.pred, add=TRUE)