spBayes (version 0.1-2)

spGLM: Function for fitting univariate Bayesian generalized linear spatial regression models

Description

The function spGLM fits univariate stationary Bayesian generalized linear spatial regression models. Given a set of knots, spGLM fits a predictive process model (see references below).

Usage

spGLM(formula, family="binomial", data = parent.frame(), coords, knots,
      starting, tuning, priors, cov.model,
      n.samples, verbose=TRUE, n.report=100, ...)

Arguments

formula

a symbolic description of the regression model to be fit. See example below.

family

currently only supports binomial and poisson data using the logit and log link functions, respectively.

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 spGLM is called.

coords

an $n \times 2$ matrix of the observation coordinates in $R^2$ (e.g., easting and northing).

knots

either a $m \times 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

starting

a list with each tag corresponding to a parameter name. Valid list tags are beta, sigma.sq, phi, nu, and w. The value portion of each tag is the parameter's starting value. If the p

tuning

a list with each tag corresponding to a parameter name. Valid list tags are beta, sigma.sq, phi, nu, and w. The value portion of each tag defines the variance of the Metropolis norma

priors

a list with each tag corresponding to a parameter name. Valid list tags are beta.flat, beta.normal, sigma.sq.ig, phi.unif, and nu.unif. simga.sq is assumed to follow an

cov.model

a quoted key word that specifies the covariance function used to model the spatial dependence structure among the observations. Supported covariance model key words are: "exponential", "matern", "spherical"

n.samples

the number of MCMC iterations.

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.

Value

• An object of class spGLM, which is a list with the following tags:
• coords
  the $n \times 2$ matrix specified by coords.
• knot.coords
  the $m \times 2$ matrix as specified by knots.
• p.samples
  a coda object of posterior samples for the defined parameters.
• acceptance
  the Metropolis sampling acceptance rate.
• sp.effects
  a matrix that holds samples from the posterior distribution of the spatial random effects. The rows of this matrix correspond to the $n$ point observations and the columns are the posterior samples.
• The return object might include additional data used for subsequent prediction and/or model fit evaluation.

References

Finley, A.O., S. Banerjee, and R.E. McRoberts. (2008) A Bayesian approach to quantifying uncertainty in multi-source forest area estimates. Environmental and Ecological Statistics, 15:241--258. 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.

Finley, A.O,. H. Sang, S. Banerjee, and A.E. Gelfand. (2008). Improving the performance of predictive process modeling for large datasets. Computational Statistics and Data Analysis, DOI: 10.1016/j.csda.2008.09.008 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.

See Also

spGGT, spMvLM, spMvGLM

Examples

###########################
##Spatial poisson
###########################

##Generate count data
n <- 100

coords <- cbind(runif(n,1,100),runif(n,1,100))

phi <- 3/75
sigma.sq <- 2

R <- sigma.sq*exp(-phi*as.matrix(dist(coords)))
w <- mvrnorm(1, rep(0,n), R)

x <- as.matrix(rep(1,n))
beta <- 0.1
y <- rpois(n, exp(x%*%beta+w))

##Collect samples
beta.starting <- coefficients(glm(y~x-1, family="poisson"))
beta.tuning <- t(chol(vcov(glm(y~x-1, family="poisson"))))
            
n.samples <- 25000

m.1 <- spGLM(y~1, family="poisson", coords=coords, 
             starting=
             list("beta"=beta.starting, "phi"=0.06,"sigma.sq"=1, "w"=0),
             tuning=
             list("beta"=0.1, "phi"=0.01, "sigma.sq"=0.01, "w"=0.005),
             priors=
             list("beta.Flat", "phi.Unif"=c(0.03, 0.3), "sigma.sq.IG"=c(2, 1)),
             cov.model="exponential",
             n.samples=n.samples, verbose=TRUE, n.report=500)

m.1$p.samples[,"phi"] <- 3/m.1$p.samples[,"phi"]


burn.in <- 0.75*n.samples
print(summary(mcmc(m.1$p.samples[burn.in:n.samples,])))

beta.hat <- mean(m.1$p.samples[burn.in:n.samples,1])
w.hat <- rowMeans(m.1$sp.effects[,burn.in:n.samples])

y.hat <-exp(x%*%beta.hat+w.hat)

##Take a look
par(mfrow=c(1,2))
surf <- mba.surf(cbind(coords,y),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="obs")
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=y, cex=1, col="blue")

surf <- mba.surf(cbind(coords,y.hat),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="Fitted counts")
contour(surf, drawlabels=FALSE, add=TRUE)
text(coords, labels=round(y.hat,0), cex=1, col="blue")

###########################
##Spatial logistic
###########################

##Generate binary data
n <- 50

coords <- cbind(runif(n,1,100),runif(n,1,100))

phi <- 3/50
sigma.sq <- 1

R <- sigma.sq*exp(-phi*as.matrix(dist(coords)))
w <- mvrnorm(1, rep(0,n), R)

x <- as.matrix(rep(1,n))
beta <- 0.1
p <- 1/(1+exp(-(x%*%beta+w)))
y <- rbinom(n, 1, prob=p)

##Collect samples
beta.starting <- coefficients(glm(y~x-1, family="binomial"))
beta.tuning <- t(chol(vcov(glm(y~x-1, family="binomial"))))
            
n.samples <- 50000

m.1 <- spGLM(y~1, family="binomial", coords=coords, 
             starting=
             list("beta"=beta.starting, "phi"=0.06,"sigma.sq"=1, "w"=0),
             tuning=
             list("beta"=beta.tuning, "phi"=0.1, "sigma.sq"=0.1, "w"=0.01),
             priors=
             list("beta.Normal"=list(0,10), "phi.Unif"=c(0.03, 0.3),
                  "sigma.sq.IG"=c(2, 1)),
             cov.model="exponential",
             n.samples=n.samples, verbose=TRUE, n.report=500)

m.1$p.samples[,"phi"] <- 3/m.1$p.samples[,"phi"]

burn.in <- 0.75*n.samples
print(summary(mcmc(m.1$p.samples[burn.in:n.samples,])))

beta.hat <- mean(m.1$p.samples[burn.in:n.samples,1])
w.hat <- rowMeans(m.1$sp.effects[,burn.in:n.samples])

y.hat <- 1/(1+exp(-(x%*%beta.hat+w.hat)))

##Take a look
par(mfrow=c(1,2))
surf <- mba.surf(cbind(coords,y),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="Observed")
contour(surf, add=TRUE)
points(coords[y==1,], pch=19, cex=1)
points(coords[y==0,], cex=1)

surf <- mba.surf(cbind(coords,y.hat),no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(surf, main="Fitted probabilities")
contour(surf, add=TRUE)
points(coords[y==1,], pch=19, cex=1)
points(coords[y==0,], cex=1)