set.seed(400)
# Simulate Data -----------------------------------------------------------
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(2:4, J, replace = TRUE)
beta <- c(0.5, 2)
p.occ <- length(beta)
alpha <- c(0, 1)
p.det <- length(alpha)
phi <- c(3 / .6, 3 / .8)
sigma.sq <- c(0.5, 0.9)
svc.cols <- c(1, 2)
dat <- simOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha,
sigma.sq = sigma.sq, phi = phi, sp = TRUE, cov.model = 'exponential',
svc.cols = svc.cols)
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .5), replace = FALSE)
y <- dat$y[-pred.indx, ]
# Occupancy covariates
X <- dat$X[-pred.indx, ]
# Prediction covariates
X.0 <- dat$X[pred.indx, ]
# Detection covariates
X.p <- dat$X.p[-pred.indx, , ]
coords <- as.matrix(dat$coords[-pred.indx, ])
coords.0 <- as.matrix(dat$coords[pred.indx, ])
psi.0 <- dat$psi[pred.indx]
w.0 <- dat$w[pred.indx, , drop = FALSE]
# Package all data into a list
occ.covs <- X[, -1, drop = FALSE]
colnames(occ.covs) <- c('occ.cov')
det.covs <- list(det.cov.1 = X.p[, , 2])
data.list <- list(y = y,
occ.covs = occ.covs,
det.covs = det.covs,
coords = coords)
# Number of batches
n.batch <- 10
# Batch length
batch.length <- 25
n.iter <- n.batch * batch.length
# Priors
prior.list <- list(beta.normal = list(mean = 0, var = 2.72),
alpha.normal = list(mean = 0, var = 2.72),
sigma.sq.ig = list(a = 2, b = 0.5),
phi.unif = list(a = 3/1, b = 3/.1))
# Initial values
inits.list <- list(alpha = 0, beta = 0,
phi = 3 / .5,
sigma.sq = 0.5,
z = apply(y, 1, max, na.rm = TRUE))
# Tuning
tuning.list <- list(phi = 1)
# Note that this is just a test case and more iterations/chains may need to
# be run to ensure convergence.
out <- svcPGOcc(occ.formula = ~ occ.cov,
det.formula = ~ det.cov.1,
data = data.list,
inits = inits.list,
n.batch = n.batch,
batch.length = batch.length,
accept.rate = 0.43,
priors = prior.list,
cov.model = 'exponential',
tuning = tuning.list,
n.omp.threads = 1,
verbose = TRUE,
NNGP = TRUE,
svc.cols = c(1, 2),
n.neighbors = 15,
search.type = 'cb',
n.report = 10,
n.burn = 50,
n.thin = 1)
summary(out)
# Predict at new locations ------------------------------------------------
out.pred <- predict(out, X.0, coords.0, verbose = FALSE)
Run the code above in your browser using DataLab