spMvLM {spBayes} | R Documentation |
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
spMvLM(formula, data = parent.frame(), coords, knots,
starting, tuning, priors, cov.model,
modified.pp = TRUE, amcmc, n.samples,
verbose=TRUE, n.report=100, ...)
Arguments
formula |
a list of |
data |
an optional data frame containing the variables in the
model. If not found in |
coords |
an |
knots |
either a |
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
|
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 |
knot.coords |
the |
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 point-referenced spatio-temporal data models. Journal of Statistical Software, 63:1–28. https://www.jstatsoft.org/article/view/v063i13.
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
Examples
## 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 cross-covariance 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.1-1, y.2~x.2-1),
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)