Function for fitting multivariate Bayesian spatial regression models
Description
The function spMvLM
fits Gaussian multivariate Bayesian
spatial regression models. Given a set of knots, spMvLM
will
also fit a predictive process model (see references below).
Usage
1 2 3 4 
Arguments
formula 
a list of q symbolic regression model descriptions to be fit. See example below. 
data 
an optional data frame containing the variables in the
model. If not found in data, the variables are taken from 
coords 
an n x 2 matrix of the observation coordinates in R^2 (e.g., easting and northing). 
knots 
either a m x 2 matrix of the predictive process knot coordinates in R^2 (e.g., easting and northing) or a vector of length two or three with the first and second elements recording the number of columns and rows in the desired knot grid. The third, optional, element sets the offset of the outermost knots from the extent of the 
starting 
a list with tags corresponding to The value portion of each tag is a vector
that holds the parameter's starting values and are of length
p for 
tuning 
a list with tags 
priors 
a list with tags 
cov.model 
a quoted keyword that specifies the covariance
function used to model the spatial dependence structure among the
observations. Supported covariance model key words are:

modified.pp 
a logical value indicating if the modified
predictive process should be used (see references below for
details). Note, if a predictive process model is not used (i.e., 
amcmc 
a list with tags 
n.samples 
the number of MCMC iterations. This argument is
ignored if 
verbose 
if 
n.report 
the interval to report Metropolis acceptance and MCMC progress. 
... 
currently no additional arguments. 
Details
Model parameters can be fixed at their starting
values by setting their
tuning
values to zero.
The no nugget model is specified by removing Psi
and L
from the starting
list.
Value
An object of class spMvLM
, which is a list with the following
tags:
coords 
the n x 2 matrix specified by

knot.coords 
the m x 2 matrix as specified by 
p.theta.samples 
a 
acceptance 
the Metropolis sampling
acceptance percent. Reported at 
The return object might include additional data used for subsequent prediction and/or model fit evaluation.
Author(s)
Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee baner009@umn.edu
References
Banerjee, S., A.E. Gelfand, A.O. Finley, and H. Sang. (2008) Gaussian Predictive Process Models for Large Spatial Datasets. Journal of the Royal Statistical Society Series B, 70:825–848.
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.
Finley, A.O., S. Banerjee, and A.E. Gelfand. (2015) spBayes for large univariate and multivariate pointreferenced spatiotemporal data models. Journal of Statistical Software, 63:1–28. http://www.jstatsoft.org/v63/i13.
Finley, A.O., H. Sang, S. Banerjee, and A.E. Gelfand. (2009) Improving the performance of predictive process modeling for large datasets. Computational Statistics and Data Analysis, 53:2873–2884.
Finley, A.O., S. Banerjee, A.R. Ek, and R.E. McRoberts. (2008) Bayesian multivariate process modeling for prediction of forest attributes. Journal of Agricultural, Biological, and Environmental Statistics, 13:60–83.
See Also
spLM
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87  ## Not run:
rmvn < function(n, mu=0, V = matrix(1)){
p < length(mu)
if(any(is.na(match(dim(V),p)))){stop("Dimension problem!")}
D < chol(V)
t(matrix(rnorm(n*p), ncol=p)%*%D + rep(mu,rep(n,p)))
}
set.seed(1)
##Generate some data
n < 25 ##number of locations
q < 2 ##number of outcomes at each location
nltr < q*(q+1)/2 ##number of triangular elements in the crosscovariance matrix
coords < cbind(runif(n,0,1), runif(n,0,1))
##Parameters for the bivariate spatial random effects
theta < rep(3/0.5,q)
A < matrix(0,q,q)
A[lower.tri(A,TRUE)] < c(1,1,0.25)
K < A%*%t(A)
Psi < diag(0,q)
C < mkSpCov(coords, K, Psi, theta, cov.model="exponential")
w < rmvn(1, rep(0,nrow(C)), C)
w.1 < w[seq(1,length(w),q)]
w.2 < w[seq(2,length(w),q)]
##Covariate portion of the mean
x.1 < cbind(1, rnorm(n))
x.2 < cbind(1, rnorm(n))
x < mkMvX(list(x.1, x.2))
B.1 < c(1,1)
B.2 < c(1,1)
B < c(B.1, B.2)
Psi < diag(c(0.1, 0.5))
y < rnorm(n*q, x%*%B+w, diag(n)%x%Psi)
y.1 < y[seq(1,length(y),q)]
y.2 < y[seq(2,length(y),q)]
##Call spMvLM
A.starting < diag(1,q)[lower.tri(diag(1,q), TRUE)]
n.samples < 1000
starting < list("phi"=rep(3/0.5,q), "A"=A.starting, "Psi"=rep(1,q))
tuning < list("phi"=rep(1,q), "A"=rep(0.01,length(A.starting)), "Psi"=rep(0.01,q))
priors < list("beta.Flat", "phi.Unif"=list(rep(3/0.75,q), rep(3/0.25,q)),
"K.IW"=list(q+1, diag(0.1,q)), "Psi.ig"=list(c(2,2), c(0.1,0.1)))
m.1 < spMvLM(list(y.1~x.11, y.2~x.21),
coords=coords, starting=starting, tuning=tuning, priors=priors,
n.samples=n.samples, cov.model="exponential", n.report=100)
burn.in < 0.75*n.samples
m.1 < spRecover(m.1, start=burn.in)
round(summary(m.1$p.theta.recover.samples)$quantiles[,c(3,1,5)],2)
round(summary(m.1$p.beta.recover.samples)$quantiles[,c(3,1,5)],2)
m.1.w.hat < summary(mcmc(t(m.1$p.w.recover.samples)))$quantiles[,c(3,1,5)]
m.1.w.1.hat < m.1.w.hat[seq(1, nrow(m.1.w.hat), q),]
m.1.w.2.hat < m.1.w.hat[seq(2, nrow(m.1.w.hat), q),]
par(mfrow=c(1,2))
plot(w.1, m.1.w.1.hat[,1], xlab="Observed w.1", ylab="Fitted w.1",
xlim=range(w), ylim=range(m.1.w.hat), main="Spatial random effects w.1")
arrows(w.1, m.1.w.1.hat[,1], w.1, m.1.w.1.hat[,2], length=0.02, angle=90)
arrows(w.1, m.1.w.1.hat[,1], w.1, m.1.w.1.hat[,3], length=0.02, angle=90)
lines(range(w), range(w))
plot(w.2, m.1.w.2.hat[,1], xlab="Observed w.2", ylab="Fitted w.2",
xlim=range(w), ylim=range(m.1.w.hat), main="Spatial random effects w.2")
arrows(w.2, m.1.w.2.hat[,1], w.2, m.1.w.2.hat[,2], length=0.02, angle=90)
arrows(w.2, m.1.w.2.hat[,1], w.2, m.1.w.2.hat[,3], length=0.02, angle=90)
lines(range(w), range(w))
## End(Not run)
