## 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 1a - Simulate a response matrix of normally distributed data
library(mvtnorm)
## 30 rows (sites) with two latent variables
true_lv <- rbind(rmvnorm(n=15,mean=c(1,2)),rmvnorm(n=15,mean=c(-3,-1)))
## 30 columns (species)
true_lv_coefs <- cbind(matrix(runif(30*3),30,3),1)
X <- matrix(rnorm(30*4),30,4)
## 4 explanatory variables
X.coefs <- matrix(rnorm(30*4),30,4)
simy <- create.life(true.lv = true_lv, lv.coefs = true_lv_coefs,
X = X, X.coefs = X.coefs, family = "normal")
if (FALSE) {
fit_simdata <- boral(simy, X = X, family = "normal", lv.control = list(num.lv = 2),
mcmc.control = example_mcmc_control, model.name = testpath)
summary(fit_simdata)
}
## Example 1b - Include a nested random row effect
## 30 subregions nested within six regions
example_row_ids <- cbind(1:30, rep(1:6,each=5))
## Subregion has a small std deviation; region has a larger one
true_ranef_sigma <- list(ID1 = 0.5, ID2 = 2)
simy <- create.life(true.lv = true_lv, lv.coefs = true_lv_coefs,
X = X, X.coefs = X.coefs, row.eff = "random",
row.params = true_ranef_sigma, row.ids = example_row_ids, family = "normal",
save.params = TRUE)
## Example 1c - Same as example 1b except new, true latent variables are generated
simy <- create.life(true.lv = NULL, lv.coefs = true_lv_coefs,
X = X, X.coefs = X.coefs, row.eff = "random",
row.params = true_ranef_sigma, row.ids = example_row_ids, family = "normal",
save.params = TRUE)
## Example 1d - Same as example 1a except new, true latent variables are generated
## with a non-independent correlation structure using a fake distance matrix
makedistmat <- as.matrix(dist(1:30))
simy <- create.life(true.lv = NULL, lv.coefs = true_lv_coefs,
lv.control = list(num.lv = 2, type = "exponential", lv.covparams = 5,
distmat = makedistmat),
X = X, X.coefs = X.coefs, row.eff = "random",
row.params = true_ranef_sigma, row.ids = example_row_ids, family = "normal",
save.params = TRUE)
## Example 1e - Similar to 1b, except instead of a nested random row effect,
## it includes a species-specific random interept at the region level
example_ranef_ids <- data.frame(region = rep(1:6,each=5))
## Subregion has a small std deviation; region has a larger one
true_ranef_sigma <- matrix(runif(nrow(true_lv_coefs)),
nrow = nrow(true_lv_coefs), ncol = 1)
simy <- create.life(true.lv = true_lv, lv.coefs = true_lv_coefs,
X = X, X.coefs = X.coefs, ranef.ids = example_ranef_ids,
ranef.params = true_ranef_sigma, family = "normal",
save.params = TRUE)
## Example 2 - Simulate a response matrix of ordinal data
## 30 rows (sites) with two latent variables
true_lv <- rbind(rmvnorm(15,mean=c(-2,-2)),rmvnorm(15,mean=c(2,2)))
## 10 columns (species)
true_lv_coefs <- rmvnorm(10,mean = rep(0,3));
## Cutoffs for proportional odds regression (must be in increasing order)
true.ordinal.cutoffs <- seq(-2,10,length=10-1)
simy <- create.life(true.lv = true_lv, lv.coefs = true_lv_coefs,
family = "ordinal", cutoffs = true.ordinal.cutoffs, save.params = TRUE)
if (FALSE) {
fit_simdata <- boral(y = simy$resp, family = "ordinal", lv.control = list(num.lv = 2),
mcmc.control = example_mcmc_control, model.name = testpath)
}
if (FALSE) {
## Example 3 - Simulate a response matrix of count data based off
## a fitted model involving traits (ants data from mvabund)
library(mvabund)
data(antTraits)
y <- antTraits$abun
X <- as.matrix(antTraits$env)
## Include only traits 1, 2, and 5, plus an intercept
traits <- as.matrix(antTraits$traits[,c(1,2,5)])
## Please see help file for boral regarding the use of which.traits
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, X = X, traits = traits, which.traits = example_which_traits,
family = "negative.binomial", lv.control = list(num.lv = 2),
mcmc.control = example_mcmc_control, model.name = testpath)
## The hard way
simy <- create.life(true.lv = fit_traits$lv.mean,
lv.coefs = fit_traits$lv.coefs.median, X = X,
X.coefs = fit_traits$X.coefs.median, traits = traits,
traits.coefs = fit_traits$traits.coefs.median, family = "negative.binomial")
## The easy way, using the same latent variables as the fitted model
simy <- simulate(object = fit_traits)
## The easy way, generating new latent variables
simy <- simulate(object = fit_traits, new.lvs = TRUE)
## Example 4 - 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(lv.coefs = cbind(rnorm(s), rmvnorm(s, mean = rep(0,2))),
lv.control = list(num.lv = 2, type = "independent"),
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(lv.control = simfit$lv.control, 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, model.name = testpath)
## Example 5 -- extend Example 2e in the main boral help file to simulate data from a
## fitted model
library(mvabund) ## Load a dataset from the mvabund package
data(spider)
y <- spider$abun
X <- scale(spider$x)
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)
simulate(object = spiderfit_nb)
simulate(object = spiderfit_nb, new.ranefs = TRUE)
}
Run the code above in your browser using DataLab