spBayes (version 0.3-8)

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, spMvGLM, spMisalignLM, or spMisalignGLM 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, spMvGLM,

pred.coords

for spLM, spMvLM, spGLM, and spMvGLM pred.coords i

pred.covars

for spLM, spMvLM, spGLM, and spMvGLM pred.covars i

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 spMvLM or spMvGLM 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.y.predictive.samples matrix has $qn$ rows. The predictions for locations are held in rows $1:q, (q+1):2q, \ldots, ((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 spMisalignLM and spMisalignGLM the posterior predictive samples are organized differently in p.y.predictive.samples with the first response variable $n_1$ locations held in rows $1\ldots,n_1$ rows, then the next response variable samples held in the $(n_1+1),\ldots,(n_1+n_2)$, etc.

  For spGLM and spMisalignGLM the p.y.predictive.samples matrix holds posterior predictive samples $\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$ is the vector of 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.

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)