The function spMisalignGLM
fits Gaussian multivariate Bayesian
generalized linear spatial regression models to misaligned data.
1 2 3  spMisalignGLM(formula, family="binomial", weights, data = parent.frame(), coords,
starting, tuning, priors, cov.model,
amcmc, n.samples, verbose=TRUE, n.report=100, ...)

formula 
a list of q symbolic regression models to be fit. See example below. 
family 
currently only supports 
weights 
an optional list of weight vectors associated with each model
in the formula list. Weights correspond to number of trials and offset for
each location for the 
data 
an optional data frame containing the variables in the
model. If not found in 
coords 
a list of q n_i x 2 matrices of the observation coordinates in R^2 (e.g., easting and northing) where i=(1,2,…,q). 
starting 
a list with tags corresponding to 
tuning 
a list with tags 
priors 
a list with each tag corresponding to a
parameter name. Valid tags are 
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:

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. 
If a binomial
model is specified the response vector is the
number of successful trials at each location and weights
is the
total number of trials at each location.
For a poisson
specification, the weights
vector is the
count offset, e.g., population, at each location. This differs from
the glm
offset
argument which is passed as the
log of this value.
An object of class spMisalignGLM
, which is a list with the following
tags:
p.beta.theta.samples 
a 
acceptance 
the Metropolis sampler
acceptance rate. If 
acceptance.w 
if 
p.w.samples 
a matrix that holds samples from the posterior distribution of the locations' spatial random effects. Posterior samples are organized with the first response variable n_1 locations held in rows 1,…,n_1 rows, then the next response variable samples held in the (n_1+1),…,(n_1+n_2), etc. 
The return object might include additional data used for subsequent prediction and/or model fit evaluation.
Andrew O. Finley finleya@msu.edu,
Sudipto Banerjee baner009@umn.edu
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 B.D. Cook. (2014) Bayesian hierarchical models for spatially misaligned data. Methods in Ecology and Evolution, 5:514–523.
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116  ## 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 < 100 ##number of locations
q < 3 ##number of outcomes at each location
nltr < q*(q+1)/2 ##number of triangular elements in the crosscovariance matrix
coords < cbind(runif(n,0,1), runif(n,0,1))
##parameters for generating a multivariate spatial GP covariance matrix
theta < rep(3/0.5,q) ##spatial decay
A < matrix(0,q,q)
A[lower.tri(A,TRUE)] < c(1,1,1,1,0.5,0.25)
K < A%*%t(A)
K ##spatial crosscovariance
cov2cor(K) ##spatial crosscorrelation
C < mkSpCov(coords, K, diag(0,q), theta, cov.model="exponential")
w < rmvn(1, rep(0,nrow(C)), C) ##spatial random effects
w.a < w[seq(1,length(w),q)]
w.b < w[seq(2,length(w),q)]
w.c < w[seq(3,length(w),q)]
##covariate portion of the mean
x.a < cbind(1, rnorm(n))
x.b < cbind(1, rnorm(n))
x.c < cbind(1, rnorm(n))
x < mkMvX(list(x.a, x.b, x.c))
B.1 < c(1,1)
B.2 < c(1,1)
B.3 < c(1,1)
B < c(B.1, B.2, B.3)
y < rpois(nrow(C), exp(x%*%B+w))
y.a < y[seq(1,length(y),q)]
y.b < y[seq(2,length(y),q)]
y.c < y[seq(3,length(y),q)]
##subsample to make spatially misaligned data
sub.1 < 1:50
y.1 < y.a[sub.1]
w.1 < w.a[sub.1]
coords.1 < coords[sub.1,]
x.1 < x.a[sub.1,]
sub.2 < 25:75
y.2 < y.b[sub.2]
w.2 < w.b[sub.2]
coords.2 < coords[sub.2,]
x.2 < x.b[sub.2,]
sub.3 < 50:100
y.3 < y.c[sub.3]
w.3 < w.c[sub.3]
coords.3 < coords[sub.3,]
x.3 < x.c[sub.3,]
##call spMisalignGLM
q < 3
A.starting < diag(1,q)[lower.tri(diag(1,q), TRUE)]
n.batch < 200
batch.length < 25
n.samples < n.batch*batch.length
starting < list("beta"=rep(0,length(B)), "phi"=rep(3/0.5,q), "A"=A.starting, "w"=0)
tuning < list("beta"=rep(0.1,length(B)), "phi"=rep(1,q), "A"=rep(0.1,length(A.starting)), "w"=1)
priors < list("phi.Unif"=list(rep(3/0.75,q), rep(3/0.25,q)),
"K.IW"=list(q+1, diag(0.1,q)), rep(0.1,q))
m.1 < spMisalignGLM(list(y.1~x.11, y.2~x.21, y.3~x.31), family="poisson",
coords=list(coords.1, coords.2, coords.3),
starting=starting, tuning=tuning, priors=priors,
amcmc=list("n.batch"=n.batch, "batch.length"=batch.length, "accept.rate"=0.43),
cov.model="exponential", n.report=10)
burn.in < floor(0.75*n.samples)
plot(m.1$p.beta.theta.samples, density=FALSE)
##predict for all locations, i.e., observed and not observed
out < spPredict(m.1, start=burn.in, thin=10, pred.covars=list(x.a, x.b, x.c),
pred.coords=list(coords, coords, coords))
##summary and check
quants < function(x){quantile(x, prob=c(0.5,0.025,0.975))}
y.hat < apply(out$p.y.predictive.samples, 1, quants)
##unstack and plot
y.a.hat < y.hat[,1:n]
y.b.hat < y.hat[,(n+1):(2*n)]
y.c.hat < y.hat[,(2*n+1):(3*n)]
par(mfrow=c(1,3))
plot(y.a ,y.a.hat[1,], xlab="Observed y.a", ylab="Fitted & predicted y.a")
plot(y.b, y.b.hat[1,], xlab="Observed y.b", ylab="Fitted & predicted y.b")
plot(y.c, y.c.hat[1,], xlab="Observed y.c", ylab="Fitted & predicted y.c")
## End(Not run)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.