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 X variables, two latent variables, and no row effects
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
lv.control = list(num.lv = 2),
save.model = TRUE, mcmc.control = example_mcmc_control,
model.name = testpath)
## Partition variance for each species into that explained by covariates
## and by the latent variables
dovar <- calc.varpart(spiderfit_nb)
## Consider the intercept and first two covariates in X as one group,
## and remaining four covariates in X as another group,
## then partition variance for each species based on these groups.
dovar <- calc.varpart(spiderfit_nb, groupX = c(1,1,1,2,2,2,2))
## Example 1b - model with X variables, two latent variables, and
## species-specific random intercepts at a so-called region level
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
lv.control = list(num.lv = 2),
ranef.ids = data.frame(subregion = rep(1:7,each=4)),
save.model = TRUE, mcmc.control = example_mcmc_control,
model.name = testpath)
## Partition variance for each species into that explained by covariates
## and by the latent variables
dovar <- calc.varpart(spiderfit_nb)
## Consider the intercept and first two covariates in X as one group,
## and remaining four covariates in X as another group,
## then partition variance for each species based on these groups.
dovar <- calc.varpart(spiderfit_nb, groupX = c(1,1,1,2,2,2,2))
## Example 2 - model fitted to count data, no site effects, and
## two latent variables, plus traits included to explain environmental responses
data(antTraits)
y <- antTraits$abun
X <- as.matrix(scale(antTraits$env))
## Include only traits 1, 2, and 5
traits <- as.matrix(antTraits$traits[,c(1,2,5)])
example_which_traits <- vector("list",ncol(X)+1)
for(i in 1:length(example_which_traits))
example_which_traits[[i]] <- 1:ncol(traits)
## Just for fun, the regression coefficients for the second column of X,
## corresponding to the third element in the list example_which_traits,
## will be estimated separately and not regressed against traits.
example_which_traits[[3]] <- 0
fit_traits <- boral(y, X = X, traits = traits, which.traits = example_which_traits,
family = "negative.binomial", mcmc.control = example_mcmc_control,
save.model = TRUE, model.name = testpath)
## Partition variance for each species due to covariates in X
## and latent variables. Also calculate proportion of variance
## due to the covariates which can be explained by traits
dovar <- calc.varpart(fit_traits)
}
Run the code above in your browser using DataLab