Function for performing posterior predictive checks
ppcAbund.RdFunction for performing posterior predictive checks on spAbundance model objects.
Arguments
- object
an object of class
NMix,spNMix,msNMix,lfMsNMix,sfMsNMix,abund,spAbund,msAbund,lfMsAbund,sfMsAbund,DS,spDS,msDS,lfMsDS,sfMsDS.- fit.stat
a quoted keyword that specifies the fit statistic to use in the posterior predictive check. Supported fit statistics are
"freeman-tukey"and"chi-squared".- group
a positive integer indicating the way to group the abundance data for the posterior predictive check. Value 0 will not group the data and use the raw counts, 1 will group values by row (site), and value 2 will group values by column (replicate).
- type
a character string indicating whether fitted values should be generated conditional on the estimated latent abundance values (
type = 'conditional') estimated during the model or based on the marginal expected abundance values (type = 'marginal'). This is only relevant for N-mixture models.- ...
currently no additional arguments
Author
Jeffrey W. Doser doserjef@msu.edu,
Note
ppcAbund will return an error for Gaussian or zero-inflated Gaussian models.
For Gaussian models, standard residual diagnostics can be used to assess model fit.
Value
An object of class ppcAbund that is a list comprised of:
- fit.y
a numeric vector of posterior samples for the fit statistic calculated on the observed data when
objectis of classNMix,spNMix,abund,spAbund,DS, orspDS. Whenobjectis of classmsNMix,lfMsNMix,sfMsNMix,msAbund,lfMsAbund,sfMsAbund,msDS,lfMsDS,sfMsDS, this is a numeric matrix with rows corresponding to posterior samples and columns corresponding to species.- fit.y.rep
a numeric vector of posterior samples for the fit statistic calculated on a replicate data set generated from the model when
objectis of classNMix,spNMix,abund,spAbund,DS,spDS. Whenobjectis of classmsNMix,lfMsNMix,sfMsNMix,msAbund,lfMsAbund,sfMsAbund,msDS,lfMsDS,sfMsDS, this is a numeric matrix with rows corresponding to posterior samples and columns corresponding to species.- fit.y.group.quants
a matrix consisting of posterior quantiles for the fit statistic using the observed data for each unique element the fit statistic is calculated for (i.e., observations when group = 0, sites when group = 1, replicates when group = 2) when
objectis of classNMix,spNMix,abund,spAbund,DS, orspDS. Whenobjectis of classmsNMix,lfMsNMix,sfMsNMix,msAbund,lfMsAbund,sfMsAbund,msDS,lfMsDS,sfMsDS, this is a three-dimensional array with the additional dimension corresponding to species.- fit.y.rep.group.quants
a matrix consisting of posterior quantiles for the fit statistic using the model replicated data for each unique element the fit statistic is calculated for (i.e., observations when group = 0, sites when group = 1, replicates when group = 2) when
objectis of classNMix,spNMix,abund,spAbund,DS,spDS. Whenobjectis of classmsNMix,sfMsNMix,msAbund,lfMsAbund,sfMsAbund,msDS,lfMsDS,sfMsDS, this is a three-dimensional array with the additional dimension corresponding to species.
The return object will include additional objects used for standard extractor functions.
Examples
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)
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> Poisson N-mixture model with 100 sites.
#>
#> Samples per Chain: 100 (4 batches of length 25)
#> Burn-in: 50
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 50
#>
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#>
#> Adaptive Metropolis with target acceptance rate: 43.0
#> ----------------------------------------
#> Chain 1
#> ----------------------------------------
#> Sampling ...
#> Batch: 1 of 4, 25.00%
#> Parameter Acceptance Tuning
#> beta[1] 12.0 0.98020
#> beta[2] 4.0 0.98020
#> alpha[1] 24.0 0.98020
#> alpha[2] 20.0 0.98020
#> alpha[3] 24.0 0.98020
#> -------------------------------------------------
#> Batch: 2 of 4, 50.00%
#> Parameter Acceptance Tuning
#> beta[1] 12.0 0.97045
#> beta[2] 12.0 0.97045
#> alpha[1] 12.0 0.97045
#> alpha[2] 28.0 0.97045
#> alpha[3] 28.0 0.97045
#> -------------------------------------------------
#> Batch: 3 of 4, 75.00%
#> Parameter Acceptance Tuning
#> beta[1] 12.0 0.96079
#> beta[2] 4.0 0.96079
#> alpha[1] 16.0 0.96079
#> alpha[2] 16.0 0.96079
#> alpha[3] 24.0 0.96079
#> -------------------------------------------------
#> Batch: 4 of 4, 100.00%
# Posterior predictive check
ppc.out <- ppcAbund(out, fit.stat = 'chi-squared', group = 0)
summary(ppc.out)
#>
#> Call:
#> ppcAbund(object = out, fit.stat = "chi-squared", group = 0)
#>
#> Samples per Chain: 100
#> Burn-in: 50
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 50
#>
#> Bayesian p-value: 0
#> Fit statistic: chi-squared