Compute Widely Applicable Information Criterion for spOccupancy Model Objects
waicOcc.RdFunction for computing the Widely Applicable Information Criterion
(WAIC; Watanabe 2010) for spOccupancy model objects.
Arguments
- object
an object of class
PGOcc,spPGOcc,msPGOcc,spMsPGOcc,intPGOcc,spIntPGOcc,lfJSDM,sfJSDM,lfMsPGOcc,sfMsPGOcc,tPGOcc,stPGOcc,svcPGBinom,svcPGOcc,svcTPGBinom,svcTPGOcc, orintMsPGOcc,svcMsPGOcc,tMsPGOcc,stMsPGOcc,svcTMsPGOcc.- by.sp
a logical value indicating whether to return a separate WAIC value for each species in a multi-species occupancy model or a single value for all species.
- ...
currently no additional arguments
References
Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11:3571-3594.
Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. (2013). Bayesian Data Analysis. 3rd edition. CRC Press, Taylor and Francis Group
Gelman, A., J. Hwang, and A. Vehtari (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24:997-1016.
Author
Jeffrey W. Doser doserjef@msu.edu,
Andrew O. Finley finleya@msu.edu
Value
When object is of class PGOcc, spPGOcc, msPGOcc, spMsPGOcc, lfJSDM, sfJSDM, lfMsPGOcc, sfMsPGOcc, tPGOcc, stPGOcc, svcPGBinom, svcPGOcc, svcTPGOcc, svcTPGBinom, svcMsPGOcc, tMsPGOcc, stMsPGOcc, svcTMsPGOcc
returns a vector with three elements corresponding to
estimates of the expected log pointwise predictive density (elpd), the
effective number of parameters (pD), and the WAIC. When by.sp = TRUE for multi-species models, object is a data frame with each row corresponding to a different species. When object is
of class intPGOcc or spIntPGOcc, returns a data frame with
columns elpd, pD, and WAIC, with each row corresponding to the estimated
values for each data source in the integrated model.
Details
The effective number of parameters is calculated following the recommendations of Gelman et al. (2014).
Examples
set.seed(400)
# Simulate Data -----------------------------------------------------------
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(2:4, J, replace = TRUE)
beta <- c(0.5, -0.15)
p.occ <- length(beta)
alpha <- c(0.7, 0.4)
p.det <- length(alpha)
dat <- simOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha,
sp = FALSE)
occ.covs <- dat$X[, 2, drop = FALSE]
colnames(occ.covs) <- c('occ.cov')
det.covs <- list(det.cov = dat$X.p[, , 2])
# Data bundle
data.list <- list(y = dat$y,
occ.covs = occ.covs,
det.covs = det.covs)
# Priors
prior.list <- list(beta.normal = list(mean = rep(0, p.occ),
var = rep(2.72, p.occ)),
alpha.normal = list(mean = rep(0, p.det),
var = rep(2.72, p.det)))
# Initial values
inits.list <- list(alpha = rep(0, p.det),
beta = rep(0, p.occ),
z = apply(data.list$y, 1, max, na.rm = TRUE))
n.samples <- 5000
n.report <- 1000
out <- PGOcc(occ.formula = ~ occ.cov,
det.formula = ~ det.cov,
data = data.list,
inits = inits.list,
n.samples = n.samples,
priors = prior.list,
n.omp.threads = 1,
verbose = TRUE,
n.report = n.report,
n.burn = 4000,
n.thin = 1)
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> Occupancy model with Polya-Gamma latent
#> variable fit with 64 sites.
#>
#> Samples per Chain: 5000
#> Burn-in: 4000
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 1000
#>
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#>
#> ----------------------------------------
#> Chain 1
#> ----------------------------------------
#> Sampling ...
#> Sampled: 1000 of 5000, 20.00%
#> -------------------------------------------------
#> Sampled: 2000 of 5000, 40.00%
#> -------------------------------------------------
#> Sampled: 3000 of 5000, 60.00%
#> -------------------------------------------------
#> Sampled: 4000 of 5000, 80.00%
#> -------------------------------------------------
#> Sampled: 5000 of 5000, 100.00%
# Calculate WAIC
waicOcc(out)
#> elpd pD WAIC
#> -110.839502 3.482209 228.643422