set.seed(1010)
J.x <- 10
J.y <- 10
J <- J.x * J.y
n.rep <- sample(3, J, replace = TRUE)
beta <- c(0, -1.5)
p.abund <- length(beta)
alpha <- c(0.5, 1.2, -0.5)
p.det <- length(alpha)
mu.RE <- list()
p.RE <- list()
phi <- 3/.6
sigma.sq <- 2
kappa <- 0.3
sp <- FALSE
cov.model <- 'exponential'
dist <- 'NB'
dat <- simNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha,
kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp,
phi = phi, sigma.sq = sigma.sq, cov.model = cov.model,
family = 'NB')
y <- dat$y
X <- dat$X
X.p <- dat$X.p
abund.covs <- X
colnames(abund.covs) <- c('int', 'abund.cov.1')
det.covs <- list(det.cov.1 = X.p[, , 2], det.cov.2 = X.p[, , 3])
data.list <- list(y = y, abund.covs = abund.covs,
det.covs = det.covs)
# Priors
prior.list <- list(beta.normal = list(mean = rep(0, p.abund),
var = rep(100, p.abund)),
alpha.normal = list(mean = rep(0, p.det),
var = rep(2.72, p.det)),
kappa.unif = c(0, 10))
# Starting values
inits.list <- list(alpha = 0, beta = 0, kappa = kappa,
N = apply(y, 1, max, na.rm = TRUE))
tuning <- 0.5
n.batch <- 4
batch.length <- 25
n.burn <- 50
n.thin <- 1
n.chains <- 1
out <- NMix(abund.formula = ~ abund.cov.1,
det.formula = ~ det.cov.1 + det.cov.2,
data = data.list,
n.batch = n.batch,
batch.length = batch.length,
inits = inits.list,
priors = prior.list,
accept.rate = 0.43,
n.omp.threads = 1,
verbose = TRUE,
n.report = 1,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
# Posterior predictive check
ppc.out <- ppcAbund(out, fit.stat = 'chi-squared', group = 0)
summary(ppc.out)
Run the code above in your browser using DataLab