set.seed(1000)
# Sites
J.x <- 10
J.y <- 10
J <- J.x * J.y
# Binomial weights
weights <- sample(10, J, replace = TRUE)
beta <- c(0, 0.5, -0.2, 0.75)
p <- length(beta)
# No unstructured random effects
psi.RE <- list()
# Spatial parameters
sp <- TRUE
# Two spatially-varying covariates.
svc.cols <- c(1, 2)
p.svc <- length(svc.cols)
cov.model <- "exponential"
sigma.sq <- runif(p.svc, 0.4, 1.5)
phi <- runif(p.svc, 3/1, 3/0.2)
# Simulate the data
dat <- simBinom(J.x = J.x, J.y = J.y, weights = weights, beta = beta,
psi.RE = psi.RE, sp = sp, svc.cols = svc.cols,
cov.model = cov.model, sigma.sq = sigma.sq, phi = phi)
# Binomial data
y <- dat$y
# Covariates
X <- dat$X
# Spatial coordinates
coords <- dat$coords
# Subset data for prediction if desired
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y.0 <- y[pred.indx, drop = FALSE]
X.0 <- X[pred.indx, , drop = FALSE]
coords.0 <- coords[pred.indx, ]
y <- y[-pred.indx, drop = FALSE]
X <- X[-pred.indx, , drop = FALSE]
coords <- coords[-pred.indx, ]
weights.0 <- weights[pred.indx]
weights <- weights[-pred.indx]
# Package all data into a list
# Covariates
covs <- cbind(X)
colnames(covs) <- c('int', 'cov.1', 'cov.2', 'cov.3')
# Data list bundle
data.list <- list(y = y,
covs = covs,
coords = coords,
weights = weights)
# Priors
prior.list <- list(beta.normal = list(mean = 0, var = 2.72),
sigma.sq.ig = list(a = 2, b = 1),
phi.unif = list(a = 3 / 1, b = 3 / 0.1))
# Starting values
inits.list <- list(beta = 0, alpha = 0,
sigma.sq = 1, phi = phi)
# Tuning
tuning.list <- list(phi = 1)
n.batch <- 10
batch.length <- 25
n.burn <- 100
n.thin <- 1
# Note that this is just a test case and more iterations/chains may need to
# be run to ensure convergence.
out <- svcPGBinom(formula = ~ cov.1 + cov.2 + cov.3,
svc.cols = c(1, 2),
data = data.list,
n.batch = n.batch,
batch.length = batch.length,
inits = inits.list,
priors = prior.list,
accept.rate = 0.43,
cov.model = "exponential",
tuning = tuning.list,
n.omp.threads = 1,
verbose = TRUE,
NNGP = TRUE,
n.neighbors = 5,
n.report = 2,
n.burn = n.burn,
n.thin = n.thin,
n.chains = 1)
summary(out)
# Predict at new locations ------------------------------------------------
out.pred <- predict(out, X.0, coords.0, weights.0, verbose = FALSE)
str(out.pred)
Run the code above in your browser using DataLab