set.seed(408)
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(3, size = J, replace = TRUE)
n.sp <- 6
# Community-level covariate effects
beta.mean <- c(-2, 0.5)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.2, 1.2)
mu.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = n.sp, ncol = p.abund)
for (i in 1:p.abund) {
beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
sp <- TRUE
kappa <- runif(n.sp, 0.1, 1)
factor.model <- TRUE
n.factors <- 3
cov.model <- 'spherical'
phi <- runif(n.factors, 3 / 1, 3 / .1)
dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta,
mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB',
factor.model = factor.model, n.factors = n.factors,
phi = phi, cov.model = cov.model)
# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y <- dat$y[, -pred.indx, , drop = FALSE]
# Occupancy covariates
X <- dat$X[-pred.indx, , , drop = FALSE]
# Coordinates
coords <- dat$coords[-pred.indx, ]
# Prediction values
y.0 <- dat$y[, pred.indx, , drop = FALSE]
X.0 <- dat$X[pred.indx, , , drop = FALSE]
coords.0 <- dat$coords[pred.indx, ]
# Package all data into a list
covs <- list(int = X[, , 1], abund.cov.1 = X[, , 2])
data.list <- list(y = y, covs = covs, coords = coords)
prior.list <- list(beta.comm.normal = list(mean = 0, var = 100),
kappa.unif = list(a = 0, b = 10),
phi.unif = list(a = 3 / 1, b = 3 / .1),
tau.sq.beta.ig = list(a = .1, b = .1))
inits.list <- list(beta.comm = 0,
beta = 0,
kappa = 0.5,
phi = 3 / .5,
tau.sq.beta = 1)
tuning.list <- list(kappa = 0.3, beta = 0.1, lambda = 0.5, w = 0.5,
phi = 1)
# Small
n.batch <- 2
batch.length <- 25
n.burn <- 20
n.thin <- 1
n.chains <- 1
out <- sfMsAbund(formula = ~ abund.cov.1,
data = data.list,
n.batch = n.batch,
inits = inits.list,
priors = prior.list,
tuning = tuning.list,
batch.length = batch.length,
n.omp.threads = 1,
n.factors = 1,
verbose = TRUE,
n.neighbors = 5,
n.report = 1,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
# Predict at new locations
out.pred <- predict(out, X.0, coords.0)
str(out.pred)
Run the code above in your browser using DataLab