spBayes (version 0.2-9)

spPredict: Function for new locations given a model object

Description

The function spPredict collects posterior predictive samples for a set of new locations given a spLM, spMvLM, spGLM, or spMvGLM object.

Usage

spPredict(sp.obj, pred.coords, pred.covars, start=1, end, thin=1,
          verbose=TRUE, n.report=100, ...)

Arguments

sp.obj

an object returned by spLM, spMvLM, spGLM, or spMvGLM.

pred.coords

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

pred.covars

an $n \times p$ design matrix associated with the new locations. If this is a multivariate prediction defined by $q$ models, i.e., from spMvLM or spMvGLM

start

specifies the first sample included in the composition sampling.

end

specifies the last sample included in the composition. The default is to use all posterior samples in sp.obj.

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

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

...

currently no additional arguments.

Value

• p.y.predictive.samples
  a matrix that holds the response variable(s) posterior predictive samples. For multivariate models the rows of this matrix correspond to the predicted locations and the columns are the posterior predictive samples. If prediction is for $q$ response variables the p.predictive.samples matrix has $qn$ rows. The predictions for locations are held in rows 1:q, (q+1):2q, ..., ((n-1)q+1):nq (i.e., the samples for the first location's $q$ response variables are in rows 1:q, second location in rows (q+1):2q, etc.). For spGLM the p.predictive.samples matrix holds posterior predictive samples for $\frac{1}{1+\exp(-x(s)'\beta-w(s))}$ and $\exp(x(s)'\beta+w(s))$ for family binomial and poisson, respectively. Here $s$ indexes the prediction location, $\beta$ are the regression coefficients, and $w$ is the associated spatial random spatial effect. These values can be fed directly into rbinom or rpois to generate the realization from the respective distribution.
• p.w.predictive.samples
  a matrix, organized the same as p.y.predictive.samples, that holds the spatial random effects posterior predictive samples. This is only returned for spGLM and spMvGLM objects.

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

Examples

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)

n <- 200
coords <- cbind(runif(n,0,1), runif(n,0,1))
X <- as.matrix(cbind(1, rnorm(n)))

B <- as.matrix(c(1,5))
p <- length(B)
sigma.sq <- 10
tau.sq <- 0.01
phi <- 3/0.5

D <- as.matrix(dist(coords))
R <- exp(-phi*D)
w <- rmvn(1, rep(0,n), sigma.sq*R)
y <- rnorm(n, X%*%B + w, sqrt(tau.sq))

##partition the data for out of sample prediction
mod <- 1:100
y.mod <- y[mod]
X.mod <- X[mod,]
coords.mod <- coords[mod,]

n.samples <- 1000

starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1)
tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)
priors <- list("beta.Flat", "phi.Unif"=c(3/1, 3/0.1),
               "sigma.sq.IG"=c(2, 5), "tau.sq.IG"=c(2, 0.01))
cov.model <- "exponential"

m.1 <- spLM(y.mod~X.mod-1, coords=coords.mod, starting=starting, tuning=tuning,
priors=priors, cov.model=cov.model, n.samples=n.samples)

m.1.pred <- spPredict(m.1, pred.covars=X, pred.coords=coords,
start=0.5*n.samples)

y.hat <- apply(m.1.pred$p.predictive.samples, 1, mean)

quant <- function(x){quantile(x, prob=c(0.025, 0.5, 0.975))}

y.hat <- apply(m.1.pred$p.predictive.samples, 1, quant)

plot(y, y.hat[2,], pch=19, cex=0.5, xlab="observed y", ylab="predicted y")
arrows(y[-mod], y.hat[2,-mod], y[-mod], y.hat[1,-mod], angle=90, length=0.05)
arrows(y[-mod], y.hat[2,-mod], y[-mod], y.hat[3,-mod], angle=90, length=0.05)