# NOT RUN {
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)
## 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)
summary(spiderfit_nb)
par(mfrow = c(2,2))
plot(spiderfit_nb) ## Plots used in residual analysis,
## Used to check if assumptions such an mean-variance relationship
## are adequately satisfied.
lvsplot(spiderfit_nb) ## Biplot of the latent variables,
## which can be interpreted in the same manner as an ordination plot.
# }
# NOT RUN {
## 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)
summary(spiderfit_nb)
## The results can be compared with the default example from
## the manyglm() function in mvabund.
## Caterpillar plots for the coefficients
par(mfrow=c(2,3), mar = c(5,6,1,1))
sapply(colnames(spiderfit_nb$X), coefsplot, x = spiderfit_nb)
## 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)
spiderfit_nb$row.coefs
## 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)
spiderfit_nb$row.coefs
## 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))
summary(spiderfit_lvstruc)
## Example 3a - Extend example 2 to demonstrate grouped covariate selection
## on the last three covariates.
example_prior_control <- list(type = c("normal","normal","normal","uniform"),
ssvs.index = c(-1,-1,-1,1,2,3))
spiderfit_nb2 <- boral(y, X = X, family = "negative.binomial",
mcmc.control = example_mcmc_control,
prior.control = example_prior_control)
summary(spiderfit_nb2)
## Example 3b - Extend example 2 to demonstrate individual covariate selection
## on the last three covariates.
example_prior_control <- list(type = c("normal","normal","normal","uniform"),
ssvs.index = c(-1,-1,-1,0,0,0))
spiderfit_nb3 <- boral(y, X = X, family = "negative.binomial",
mcmc.control = example_mcmc_control,
prior.control = example_prior_control)
summary(spiderfit_nb3)
## Example 4 - model fitted to presence-absence data, no site effects, and
## two latent variables
data(tikus)
y <- tikus$abun
y[y > 0] <- 1
y <- y[1:20,] ## Consider only years 1981 and 1983
y <- y[,apply(y > 0,2,sum) > 2] ## Consider only spp with more than 2 presences
tikus.fit <- boral(y, family = "binomial",
lv.control = list(num.lv = 2),
mcmc.control = example_mcmc_control)
lvsplot(tikus.fit, biplot = FALSE)
## A strong location between the two sampling years
## Example 5a - 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,
lv.control = list(num.lv = 2),
which.traits = example_which_traits, family = "negative.binomial",
mcmc.control = example_mcmc_control, save.model = TRUE)
summary(fit_traits)
## Example 5b - perform selection on trait coefficients
ssvs_traitsindex <- vector("list",ncol(X)+1)
for(i in 1:length(ssvs_traitsindex))
ssvs_traitsindex[[i]] <- rep(0,ncol(traits))
ssvs_traitsindex[[3]] <- -1
fit_traits <- boral(y, X = X, traits = traits, which.traits = example_which_traits,
family = "negative.binomial",
lv.control = list(num.lv = 2), mcmc.control = example_mcmc_control,
save.model = TRUE, prior.control = list(ssvs.traitsindex = ssvs_traitsindex))
summary(fit_traits)
## Example 6 - simulate Bernoulli data, based on a model with two latent variables,
## no site variables, with two traits and one environmental covariates
## This example is a proof of concept that traits can used to
## explain environmental responses
library(mvtnorm)
n <- 100; s <- 50
X <- as.matrix(scale(1:n))
colnames(X) <- c("elevation")
traits <- cbind(rbinom(s,1,0.5), rnorm(s))
## one categorical and one continuous variable
colnames(traits) <- c("thorns-dummy","SLA")
simfit <- list(true.lv = rmvnorm(n, mean = rep(0,2)),
lv.coefs = cbind(rnorm(s), rmvnorm(s, mean = rep(0,2))),
traits.coefs = matrix(c(0.1,1,-0.5,1,0.5,0,-1,1), 2, byrow = TRUE))
rownames(simfit$traits.coefs) <- c("beta0","elevation")
colnames(simfit$traits.coefs) <- c("kappa0","thorns-dummy","SLA","sigma")
simy = create.life(true.lv = simfit$true.lv, lv.coefs = simfit$lv.coefs, X = X,
traits = traits, traits.coefs = simfit$traits.coefs, family = "binomial")
example_which_traits <- vector("list",ncol(X)+1)
for(i in 1:length(example_which_traits))
example_which_traits[[i]] <- 1:ncol(traits)
fit_traits <- boral(y = simy, X = X, traits = traits,
which.traits = example_which_traits, family = "binomial",
lv.control = list(num.lv = 2), save.model = TRUE,
mcmc.control = example_mcmc_control)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab