if (FALSE) {
library(mvabund) ## Load a dataset from the mvabund package
library(mvtnorm)
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, random site effects,
## and environmental covariates
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
row.eff = "random", lv.control = list(num.lv = 2),
mcmc.control = example_mcmc_control, save.model = TRUE,
model.name = testpath)
## Predictions conditional on predicted latent variables
getcondpreds <- predict(spiderfit_nb)
## Predictions marginal on latent variables, random row effects
## The intervals for these will generally be wider than the
## conditional intervals.
getmargpreds <- predict(spiderfit_nb, predict.type = "marginal")
## Now suppose you extrpolate to new sites
newX <- rmvnorm(100, mean = rep(0,ncol(X)))
## Below won't work since conditional predictions are made to the same sites
#getcondpreds <- predict(spiderfit_nb, newX = newX)
## Marginal predictions will work though, provided newrow.ids is set up
## properly. For example,
new_row_ids <- matrix(sample(1:28,100,replace=TRUE), 100, 1)
while(length(table(new_row_ids)) != 28) {
new_row_ids <- matrix(sample(1:28,100,replace=TRUE), 100, 1)
}
getmargpreds <- predict(spiderfit_nb, newX = newX, predict.type = "marginal",
newrow.ids = new_row_ids)
## Example 1b - Similar to 1 except with no random site effects,
## and a non-independence correlation structure for the latent variables
## based on a fake distance matrix
fakedistmat <- as.matrix(dist(1:n))
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
lv.control = list(type = "squared.exponential", num.lv = 2,
distmat = fakedistmat), model.name = testpath,
mcmc.control = example_mcmc_control, save.model = TRUE)
getmargpreds <- predict(spiderfit_nb, predict.type = "marginal",
distmat = fakedistmat)
## Now suppose you extrpolate to new sites
newfakedistmat <- as.matrix(dist(1:100))
getmargpreds <- predict(spiderfit_nb, newX = newX,
predict.type = "marginal", distmat = newfakedistmat)
## Example 1c - similar to 1 except instead of random site effects,
## there are species-specific random intercepts at a so-called
## "region" level
y <- spider$abun
X <- scale(spider$x)
spiderfit_nb <- boral(y, X = X, family = "negative.binomial",
lv.control = list(num.lv = 2),
ranef.ids = data.frame(region = rep(1:7,each=4)),
mcmc.control = example_mcmc_control, model.name = testpath,
save.model = TRUE)
## Predictions conditional on predicted latent variables and
## random intercepts
getcondpreds <- predict(spiderfit_nb)
## Predictions marginal on latent variables, random intercepts
## The intervals for these will generally be wider than the
## conditional intervals.
getmargpreds <- predict(spiderfit_nb, predict.type = "marginal")
## Now suppose you extrpolate to new sites
newX <- rmvnorm(100, mean = rep(0,ncol(X)))
## Marginal predictions will work though, provided newranef.ids is set up
## properly. For example,
new_ranef_ids <- matrix(sample(1:7,100,replace=TRUE), 100, 1)
getmargpreds <- predict(spiderfit_nb, newX = newX, predict.type = "marginal",
newranef.ids = new_ranef_ids)
## Example 2 - simulate count data, based on a model with two latent variables,
## no site variables, with two traits and one environmental covariates
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)), 1),
traits.coefs = matrix(c(0.1,1,-0.5,0.1,0.5,0,-1,0.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 = "normal")
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 = "normal",
lv.control = list(num.lv = 2), save.model = TRUE,
mcmc.control = example_mcmc_control,
model.name = testpath)
## Predictions conditional on predicted latent variables
getcondpreds <- predict(fit_traits)
## Predictions marginal on latent variables
## The intervals for these will generally be wider than the
## conditional intervals.
getmargpreds <- predict(fit_traits, predict.type = "marginal")
}
Run the code above in your browser using DataLab