if (FALSE) {
## 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")
library(mvabund) ## Load a dataset from the mvabund package
data(spider)
y <- spider$abun
n <- nrow(y)
p <- ncol(y)
X <- scale(spider$x)
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
lv.control = list(num.lv = 2), mcmc.control = example_mcmc_control,
model.name = testpath)
## Do separate line plots for all the coefficients of X
par(mfrow=c(2,3), mar = c(5,6,1,1))
sapply(colnames(spiderfit_nb$X), coefsplot,
spiderfit_nb)
## Consider a model based on Example 5a in the main boral help file
## The model is fitted to count data, no site effects, 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, model.name = testpath,
save.model = TRUE)
summary(fit_traits)
par(mar = c(3,10,2,10))
coefsplot(object = fit_traits, fourthcorner = TRUE)
}
Run the code above in your browser using DataLab