if (FALSE) {
library(mvabund) ## Load a dataset from the mvabund package
data(spider)
y <- spider$abun
X <- scale(spider$x)
n <- nrow(y)
p <- ncol(y)
## NOTE: The values below MUST NOT be used in a real application;
## they are only used here to make the examples run quick!!!
example_mcmc_control <- list(n.burnin = 10, n.iteration = 100,
n.thin = 1)
testpath <- file.path(tempdir(), "jagsboralmodel.txt")
## Example 1 - model with two latent variables, site effects,
## and no environmental covariates
spiderfit_nb <- boral(y, family = "negative.binomial",
lv.control = list(num.lv = 2), row.eff = "fixed",
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
## Example 2a - model with no latent variables, no site effects,
## and environmental covariates
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
## Example 2b - suppose now, for some reason, the 28 rows were
## sampled such into four replications of seven sites
## Let us account for this as a fixed effect
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
row.eff = "fixed", row.ids = matrix(rep(1:7,each=4),ncol=1),
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
## Example 2c - suppose now, for some reason, the 28 rows reflected
## a nested design with seven regions, each with four sub-regions
## We can account for this nesting as a random effect
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
row.eff = "random",
row.ids = cbind(1:n, rep(1:7,each=4)),
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
## Example 2d - model with environmental covariates and
## two structured latent variables using fake distance matrix
fakedistmat <- as.matrix(dist(1:n))
spiderfit_lvstruc <- boral(y, X = X, family = "negative.binomial",
lv.control = list(num.lv = 2, type = "exponential", distmat = fakedistmat),
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
## Example 2e - Similar to 2d, but we will species-specific random intercepts
## for the seven regions (with row effects in the model)
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
ranef.ids = data.frame(region = rep(1:7,each=4)),
mcmc.control = example_mcmc_control, model.name = testpath)
spiderfit_nb$hpdintervals
}
Run the code above in your browser using DataLab