spBayes: spMvLM – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

spMvLM

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 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 `environment(formula)`, typically the environment from which `spMvLM` is called.
`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 `coords`.
`starting`	a list with tags corresponding to `beta`, `A`, `phi`, and `nu`. Depending on the specification of the non-spatial residual, tags are `L` or `Psi` for a block diagonal or diagonal covariance matrix, respectively. The value portion of each tag is a vector that holds the parameter's starting values and are of length p for `beta` (where p is the total number of regression coefficients in the multivariate model), q(q+1)/2 for `A` and `L`, and q for `Psi`, `phi`, and `nu`. Here, `A` and `L` hold the lower-triangle elements in column major ordering of the Cholesky square root of the spatial and non-spatial cross-covariance matrices, respectively.
`tuning`	a list with tags `A`, `phi`, and `nu`. Depending on the specification of the non-spatial residual, tags are `L` or `Psi` for a block diagonal or diagonal covariance matrix, respectively. The value portion of each tag defines the variance of the Metropolis sampler Normal proposal distribution. For `A` and `L` the vectors are of length q(q+1)/2 and q for `Psi`, `phi`, and `nu`.
`priors`	a list with tags `beta.flat`, `beta.norm`, `K.iw`, `Psi.iw`, `Psi.ig`, `phi.unif`, and `nu.unif`. If the regression coefficients, i.e., `beta` vector, are assumed to follow a multivariate Normal distribution then pass the hyperparameters as a list of length two with the first and second elements corresponding to the mean vector and positive definite covariance matrix, respectively. If `beta` is assumed flat then no arguments are passed. The default is a flat prior. Use `Psi.iw` if the non-spatial residual covariance matrix is assumed block diagonal. Otherwise if the non-spatial residual covariance matrix is assumed diagonal then each of the q diagonal element are assumed to follow an inverse-Gamma in which case use `Psi.ig`. The hyperparameters of the inverse-Wishart, i.e., for cross-covariance matrices AA' `K.iw` and LL' `Psi.iw`, are passed as a list of length two, with the first and second elements corresponding to the df and qxq scale matrix, respectively. If `Psi.ig` is specified, the inverse-Gamma hyperparameters of the diagonal variance elements are pass using a list of length two with the first and second list elements consisting of vectors of the q shape and scale hyperparameters, respectively. The hyperparameters of the Uniform `phi.unif`, and `nu.unif` are also passed as a list of vectors with the first and second list elements corresponding to the lower and upper support, respectively.
`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: `"exponential"`, `"matern"`, `"spherical"`, and `"gaussian"`. See below for details.
`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., `knots` is not specified) then this argument is ignored.
`amcmc`	a list with tags `n.batch`, `batch.length`, and `accept.rate`. Specifying this argument invokes an adaptive MCMC sampler see Roberts and Rosenthal (2007) for an explanation.
`n.samples`	the number of MCMC iterations. This argument is ignored if `amcmc` is specified.
`verbose`	if `TRUE`, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
`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 `coords`.
`knot.coords`	the m x 2 matrix as specified by `knots`.
`p.theta.samples`	a `coda` object of posterior samples for the defined parameters.
`acceptance`	the Metropolis sampling acceptance percent. Reported at `batch.length` or `n.report` intervals for `amcmc` specified and non-specified, respectively

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

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)

spBayes

Univariate and Multivariate Spatial-Temporal Modeling

v0.4-3

GPL (>= 2)

Authors

Andrew Finley [aut, cre], Sudipto Banerjee [aut]

Initial release

2020-02-06