tests/testthat/test-occumb.R
In fukayak/occumb: Site Occupancy Modeling for Environmental DNA Metabarcoding
### Tests for quality controls -------------------------------------------------
test_that("Checks for data work", {
expect_error(occumb(data = array(1, dim = rep(2, 3))),
"An occumbData class object is expected for data")
})
test_that("Checks for formula work", {
formulas <- c("formula_phi",
"formula_theta",
"formula_psi",
"formula_phi_shared",
"formula_theta_shared",
"formula_psi_shared")
expect_error(occumb(1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_phi"))
expect_error(occumb(~ 1, 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_theta"))
expect_error(occumb(~ 1, ~ 1, 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_psi"))
expect_error(occumb(~ 1, ~ 1, ~ 1, 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_phi_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_theta_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_psi_shared"))
expect_error(occumb(1, 1, 1, 1, 1, 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s",
paste(formulas, collapse = ", ")))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 0, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_phi_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ -1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_phi_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 0, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_theta_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ -1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_theta_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 0,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ -1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi_shared"))
expect_error(occumb(~ 0, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "phi"))
expect_error(occumb(~ -1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "phi"))
expect_error(occumb(~ 1, ~ 0, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "theta"))
expect_error(occumb(~ 1, ~ -1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "theta"))
expect_error(occumb(~ 1, ~ 1, ~ 0, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi"))
expect_error(occumb(~ 1, ~ 1, ~ -1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi"))
})
test_that("Checks for priors work", {
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = Inf),
"'prior_prec' should have a finite value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = 0),
"'prior_prec' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = -0.01),
"'prior_prec' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_ulim = 0),
"'prior_ulim' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_ulim = -0.01),
"'prior_ulim' should have a positive value")
})
### Tests for set_const() ------------------------------------------------------
test_that("Replacement of missing y to NA works", {
y <- array(1:8, dim = rep(2, 3))
y[, 1, 1] <- 0
result <- set_const(occumbData(y = y))
expect_identical(result$y[, 1, 1], as.numeric(rep(NA, 2)))
expect_identical(result$y[, 2, 1], c(3, 4))
expect_identical(result$y[, 1, 2], c(5, 6))
expect_identical(result$y[, 2, 2], c(7, 8))
})
test_that("Replacement of sequence depth 0 to 1 works", {
y <- array(1:8, dim = rep(2, 3))
y[, 1, 1] <- 0
result <- set_const(occumbData(y = y))
expect_equal(result$N[1, 1], 1)
expect_equal(result$N[2, 1], sum(3:4))
expect_equal(result$N[1, 2], sum(5:6))
expect_equal(result$N[2, 2], sum(7:8))
})
### Test cases for set_data/inits_code/set_params_monitored/write_jags_model ---
phi <- theta <- c("i", "ij", "ijk")
psi <- c("i", "ij")
phi_shared <- theta_shared <- psi_shared <- c(FALSE, TRUE)
cases <- expand.grid(phi, theta, psi, phi_shared, theta_shared, psi_shared)
colnames(cases) <- c("phi", "theta", "psi",
"phi_shared", "theta_shared", "psi_shared")
### Tests for set_data() -------------------------------------------------------
test_that("Data list is correct for 144 available models", {
const <- list(I = sample.int(1E3, 1),
J = sample.int(1E3, 1),
K = sample.int(1E3, 1),
N = sample.int(1E3, 1),
y = sample.int(1E3, 1))
prior_prec <- rnorm(1)
prior_ulim <- rnorm(1)
for (i in 1:nrow(cases)) {
margs <- list(cov_phi = rnorm(1),
cov_theta = rnorm(1),
cov_psi = rnorm(1),
cov_phi_shared = rnorm(1),
cov_theta_shared = rnorm(1),
cov_psi_shared = rnorm(1),
M = sample.int(1E3, 1),
M_phi_shared = sample.int(1E3, 1),
M_theta_shared = sample.int(1E3, 1),
M_psi_shared = sample.int(1E3, 1),
m_phi = sample.int(1E3, 1),
m_theta = sample.int(1E3, 1),
m_psi = sample.int(1E3, 1),
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
ans <- list(I = const$I,
J = const$J,
K = const$K,
N = const$N,
y = const$y,
cov_phi = margs$cov_phi,
cov_theta = margs$cov_theta,
cov_psi = margs$cov_psi,
M = margs$M,
m_phi = margs$m_phi,
m_theta = margs$m_theta,
m_psi = margs$m_psi,
prior_prec = prior_prec,
prior_ulim = prior_ulim)
if (cases$phi_shared[i]) {
ans$cov_phi_shared <- margs$cov_phi_shared
ans$M_phi_shared <- margs$M_phi_shared
}
if (cases$theta_shared[i]) {
ans$cov_theta_shared <- margs$cov_theta_shared
ans$M_theta_shared <- margs$M_theta_shared
}
if (cases$psi_shared[i]) {
ans$cov_psi_shared <- margs$cov_psi_shared
ans$M_psi_shared <- margs$M_psi_shared
}
res <- set_data(const, margs, prior_prec, prior_ulim)
expect_equal(res, ans)
}
})
### Tests for inits_code() -----------------------------------------------------
test_that("Code for initial value function is correct for 144 available models", {
const <- list(I = sample.int(1E3, 1),
J = sample.int(1E3, 1),
K = sample.int(1E3, 1),
N = sample.int(1E3, 1),
y = sample.int(1E3, 1))
dat <- list()
for (i in 1:nrow(cases)) {
margs <- list(cov_phi = rnorm(1),
cov_theta = rnorm(1),
cov_psi = rnorm(1),
cov_phi_shared = rnorm(1),
cov_theta_shared = rnorm(1),
cov_psi_shared = rnorm(1),
M = sample.int(1E3, 1),
M_phi_shared = sample.int(1E3, 1),
M_theta_shared = sample.int(1E3, 1),
M_psi_shared = sample.int(1E3, 1),
m_phi = sample.int(1E3, 1),
m_theta = sample.int(1E3, 1),
m_psi = sample.int(1E3, 1),
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
# Tests for prior_ulim > 10 cases
dat$prior_ulim <- runif(1, min = 10, max = 1000)
ans <- c("function() {",
" list(z = matrix(1, const$I, const$J),",
" u = array(1, dim = c(const$I, const$J, const$K)),",
" r = array(stats::rnorm(const$I * const$J * const$K,",
" mean = 1, sd = 0.1),",
" dim = c(const$I, const$J, const$K)),")
if (cases$phi_shared[i]) {
ans <- c(ans, " alpha_shared = stats::rnorm(margs$M_phi_shared, sd = 0.1),")
}
if (cases$theta_shared[i]) {
ans <- c(ans, " beta_shared = stats::rnorm(margs$M_theta_shared, sd = 0.1),")
}
if (cases$psi_shared[i]) {
ans <- c(ans, " gamma_shared = stats::rnorm(margs$M_psi_shared, sd = 0.1),")
}
ans <- c(ans,
" spec_eff = matrix(stats::rnorm(const$I * margs$M, sd = 0.1),",
" const$I, margs$M),",
" Mu = stats::rnorm(margs$M, sd = 0.1),",
" sigma = stats::rnorm(margs$M, mean = 1, sd = 0.1))",
"}")
res <- inits_code(const, margs, dat)
expect_equal(res, ans)
# Tests for prior_ulim < 10 cases
dat$prior_ulim <- runif(1, min = 0, max = 10)
ans <- c("function() {",
" list(z = matrix(1, const$I, const$J),",
" u = array(1, dim = c(const$I, const$J, const$K)),",
" r = array(stats::rnorm(const$I * const$J * const$K,",
" mean = 1, sd = 0.1),",
" dim = c(const$I, const$J, const$K)),")
if (cases$phi_shared[i]) {
ans <- c(ans, " alpha_shared = stats::rnorm(margs$M_phi_shared, sd = 0.1),")
}
if (cases$theta_shared[i]) {
ans <- c(ans, " beta_shared = stats::rnorm(margs$M_theta_shared, sd = 0.1),")
}
if (cases$psi_shared[i]) {
ans <- c(ans, " gamma_shared = stats::rnorm(margs$M_psi_shared, sd = 0.1),")
}
ans <- c(ans,
" spec_eff = matrix(stats::rnorm(const$I * margs$M, sd = 0.1),",
" const$I, margs$M),",
" Mu = stats::rnorm(margs$M, sd = 0.1),")
ans <- c(ans,
sprintf(" sigma = stats::runif(margs$M, min = 0, max = %s))",
dat$prior_ulim))
ans <- c(ans, "}")
res <- inits_code(const, margs, dat)
expect_equal(res, ans)
}
})
### Tests for set_params_monitored() -------------------------------------------
test_that("Parameter list is correct for 144 available models", {
for (i in 1:nrow(cases)) {
ans <- c("Mu", "sigma", "rho", "alpha", "beta", "gamma")
if (cases$phi_shared[i]) {
ans <- c(ans, "alpha_shared")
}
if (cases$theta_shared[i]) {
ans <- c(ans, "beta_shared")
}
if (cases$psi_shared[i]) {
ans <- c(ans, "gamma_shared")
}
ans <- c(ans, c("phi", "theta", "psi", "z", "pi"))
res <- set_params_monitored(phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
expect_equal(res, ans)
}
})
### Tests for write_jags_model() -----------------------------------------------
test_that("JAGS code is correct for 144 available models", {
for (i in 1:nrow(cases)) {
ans <- readLines(system.file("jags",
"occumb_template1.jags",
package = "occumb"))
if (cases$phi[i] == "i") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i], 1)")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i, j], 1)")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i, j, k], 1)")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template2.jags",
package = "occumb")))
if (cases$theta[i] == "i") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i, j])")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i, j, k])")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template3.jags",
package = "occumb")))
if (cases$psi[i] == "i") {
ans <- c(ans,
" z[i, j] ~ dbern(psi[i])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" z[i, j] ~ dbern(psi[i, j])")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template4.jags",
package = "occumb")))
if (cases$phi_shared[i]) {
if (cases$phi[i] == "i") {
ans <- c(ans,
" log(phi[i]) <- inprod(alpha[i, ], cov_phi[]) + inprod(alpha_shared[], cov_phi_shared[i, ])")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" log(phi[i, j]) <- inprod(alpha[i, ], cov_phi[j, ]) + inprod(alpha_shared[], cov_phi_shared[i, j, ])",
" }")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" log(phi[i, j, k]) <- inprod(alpha[i, ], cov_phi[j, k, ]) + inprod(alpha_shared[], cov_phi_shared[i, j, k, ])",
" }",
" }")
}
} else {
if (cases$phi[i] == "i") {
ans <- c(ans,
" log(phi[i]) <- inprod(alpha[i, ], cov_phi[])")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" log(phi[i, j]) <- inprod(alpha[i, ], cov_phi[j, ])",
" }")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" log(phi[i, j, k]) <- inprod(alpha[i, ], cov_phi[j, k, ])",
" }",
" }")
}
}
if (cases$theta_shared[i]) {
if (cases$theta[i] == "i") {
ans <- c(ans,
" logit(theta[i]) <- inprod(beta[i, ], cov_theta[]) + inprod(beta_shared[], cov_theta_shared[i, ])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(theta[i, j]) <- inprod(beta[i, ], cov_theta[j, ]) + inprod(beta_shared[], cov_theta_shared[i, j, ])",
" }")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" logit(theta[i, j, k]) <- inprod(beta[i, ], cov_theta[j, k, ]) + inprod(beta_shared[], cov_theta_shared[i, j, k, ])",
" }",
" }")
}
} else {
if (cases$theta[i] == "i") {
ans <- c(ans,
" logit(theta[i]) <- inprod(beta[i, ], cov_theta[])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(theta[i, j]) <- inprod(beta[i, ], cov_theta[j, ])",
" }")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" logit(theta[i, j, k]) <- inprod(beta[i, ], cov_theta[j, k, ])",
" }",
" }")
}
}
if (cases$psi_shared[i]) {
if (cases$psi[i] == "i") {
ans <- c(ans,
" logit(psi[i]) <- inprod(gamma[i, ], cov_psi[]) + inprod(gamma_shared[], cov_psi_shared[i, ])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(psi[i, j]) <- inprod(gamma[i, ], cov_psi[j, ]) + inprod(gamma_shared[], cov_psi_shared[i, j, ])",
" }")
}
} else {
if (cases$psi[i] == "i") {
ans <- c(ans,
" logit(psi[i]) <- inprod(gamma[i, ], cov_psi[])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(psi[i, j]) <- inprod(gamma[i, ], cov_psi[j, ])",
" }")
}
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template5.jags",
package = "occumb")))
if (cases$phi_shared[i]) {
ans <- c(ans,
" for (m in 1:M_phi_shared) {",
" alpha_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
if (cases$theta_shared[i]) {
ans <- c(ans,
" for (m in 1:M_theta_shared) {",
" beta_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
if (cases$psi_shared[i]) {
ans <- c(ans,
" for (m in 1:M_psi_shared) {",
" gamma_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
ans <- c(ans, "}", "")
res <- write_jags_model(phi = cases$phi[i],
theta = cases$theta[i],
psi = cases$psi[i],
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
expect_equal(res, ans)
}
})
### Test for get_data() --------------------------------------------------------
test_that("Outputs are correct when proper variable names are given", {
I <- 2
J <- 2
K <- 2
y <- array(sample.int(I * J * K), dim = c(I, J, K))
spec_cov <- list(cov1 = rnorm(I))
site_cov <- list(cov2 = rnorm(J), cov3 = factor(1:J))
repl_cov <- list(cov4 = matrix(rnorm(J * K), J, K))
data <- occumbData(y = y,
spec_cov = spec_cov,
site_cov = site_cov,
repl_cov = repl_cov)
fit <- occumb(data = data,
n.chains = 1, n.burnin = 10, n.thin = 1, n.iter = 20,
verbose = FALSE)
expect_identical(get_data(fit, "y"), y)
expect_identical(get_data(fit, "spec_cov"), spec_cov)
expect_identical(get_data(fit, "site_cov"), site_cov)
expect_identical(get_data(fit, "repl_cov"), repl_cov)
})
fukayak/occumb documentation built on April 17, 2025, 11:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### Tests for quality controls -------------------------------------------------
test_that("Checks for data work", {
expect_error(occumb(data = array(1, dim = rep(2, 3))),
"An occumbData class object is expected for data")
})
test_that("Checks for formula work", {
formulas <- c("formula_phi",
"formula_theta",
"formula_psi",
"formula_phi_shared",
"formula_theta_shared",
"formula_psi_shared")
expect_error(occumb(1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_phi"))
expect_error(occumb(~ 1, 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_theta"))
expect_error(occumb(~ 1, ~ 1, 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_psi"))
expect_error(occumb(~ 1, ~ 1, ~ 1, 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_phi_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_theta_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s", "formula_psi_shared"))
expect_error(occumb(1, 1, 1, 1, 1, 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("Formula is expected for: %s",
paste(formulas, collapse = ", ")))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 0, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_phi_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ -1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_phi_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 0, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_theta_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ -1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
"No intercept in formula_theta_shared: remove 0 or -1 from the formula")
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 0,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi_shared"))
expect_error(occumb(~ 1, ~ 1, ~ 1, ~ 1, ~ 1, ~ -1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi_shared"))
expect_error(occumb(~ 0, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "phi"))
expect_error(occumb(~ -1, ~ 1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "phi"))
expect_error(occumb(~ 1, ~ 0, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "theta"))
expect_error(occumb(~ 1, ~ -1, ~ 1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "theta"))
expect_error(occumb(~ 1, ~ 1, ~ 0, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi"))
expect_error(occumb(~ 1, ~ 1, ~ -1, ~ 1, ~ 1, ~ 1,
data = occumbData(y = array(1, dim = rep(2, 3)))),
sprintf("No intercept in formula_%s: remove 0 or -1 from the formula", "psi"))
})
test_that("Checks for priors work", {
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = Inf),
"'prior_prec' should have a finite value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = 0),
"'prior_prec' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_prec = -0.01),
"'prior_prec' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_ulim = 0),
"'prior_ulim' should have a positive value")
expect_error(occumb(data = occumbData(y = array(1, dim = rep(2, 3))),
prior_ulim = -0.01),
"'prior_ulim' should have a positive value")
})
### Tests for set_const() ------------------------------------------------------
test_that("Replacement of missing y to NA works", {
y <- array(1:8, dim = rep(2, 3))
y[, 1, 1] <- 0
result <- set_const(occumbData(y = y))
expect_identical(result$y[, 1, 1], as.numeric(rep(NA, 2)))
expect_identical(result$y[, 2, 1], c(3, 4))
expect_identical(result$y[, 1, 2], c(5, 6))
expect_identical(result$y[, 2, 2], c(7, 8))
})
test_that("Replacement of sequence depth 0 to 1 works", {
y <- array(1:8, dim = rep(2, 3))
y[, 1, 1] <- 0
result <- set_const(occumbData(y = y))
expect_equal(result$N[1, 1], 1)
expect_equal(result$N[2, 1], sum(3:4))
expect_equal(result$N[1, 2], sum(5:6))
expect_equal(result$N[2, 2], sum(7:8))
})
### Test cases for set_data/inits_code/set_params_monitored/write_jags_model ---
phi <- theta <- c("i", "ij", "ijk")
psi <- c("i", "ij")
phi_shared <- theta_shared <- psi_shared <- c(FALSE, TRUE)
cases <- expand.grid(phi, theta, psi, phi_shared, theta_shared, psi_shared)
colnames(cases) <- c("phi", "theta", "psi",
"phi_shared", "theta_shared", "psi_shared")
### Tests for set_data() -------------------------------------------------------
test_that("Data list is correct for 144 available models", {
const <- list(I = sample.int(1E3, 1),
J = sample.int(1E3, 1),
K = sample.int(1E3, 1),
N = sample.int(1E3, 1),
y = sample.int(1E3, 1))
prior_prec <- rnorm(1)
prior_ulim <- rnorm(1)
for (i in 1:nrow(cases)) {
margs <- list(cov_phi = rnorm(1),
cov_theta = rnorm(1),
cov_psi = rnorm(1),
cov_phi_shared = rnorm(1),
cov_theta_shared = rnorm(1),
cov_psi_shared = rnorm(1),
M = sample.int(1E3, 1),
M_phi_shared = sample.int(1E3, 1),
M_theta_shared = sample.int(1E3, 1),
M_psi_shared = sample.int(1E3, 1),
m_phi = sample.int(1E3, 1),
m_theta = sample.int(1E3, 1),
m_psi = sample.int(1E3, 1),
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
ans <- list(I = const$I,
J = const$J,
K = const$K,
N = const$N,
y = const$y,
cov_phi = margs$cov_phi,
cov_theta = margs$cov_theta,
cov_psi = margs$cov_psi,
M = margs$M,
m_phi = margs$m_phi,
m_theta = margs$m_theta,
m_psi = margs$m_psi,
prior_prec = prior_prec,
prior_ulim = prior_ulim)
if (cases$phi_shared[i]) {
ans$cov_phi_shared <- margs$cov_phi_shared
ans$M_phi_shared <- margs$M_phi_shared
}
if (cases$theta_shared[i]) {
ans$cov_theta_shared <- margs$cov_theta_shared
ans$M_theta_shared <- margs$M_theta_shared
}
if (cases$psi_shared[i]) {
ans$cov_psi_shared <- margs$cov_psi_shared
ans$M_psi_shared <- margs$M_psi_shared
}
res <- set_data(const, margs, prior_prec, prior_ulim)
expect_equal(res, ans)
}
})
### Tests for inits_code() -----------------------------------------------------
test_that("Code for initial value function is correct for 144 available models", {
const <- list(I = sample.int(1E3, 1),
J = sample.int(1E3, 1),
K = sample.int(1E3, 1),
N = sample.int(1E3, 1),
y = sample.int(1E3, 1))
dat <- list()
for (i in 1:nrow(cases)) {
margs <- list(cov_phi = rnorm(1),
cov_theta = rnorm(1),
cov_psi = rnorm(1),
cov_phi_shared = rnorm(1),
cov_theta_shared = rnorm(1),
cov_psi_shared = rnorm(1),
M = sample.int(1E3, 1),
M_phi_shared = sample.int(1E3, 1),
M_theta_shared = sample.int(1E3, 1),
M_psi_shared = sample.int(1E3, 1),
m_phi = sample.int(1E3, 1),
m_theta = sample.int(1E3, 1),
m_psi = sample.int(1E3, 1),
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
# Tests for prior_ulim > 10 cases
dat$prior_ulim <- runif(1, min = 10, max = 1000)
ans <- c("function() {",
" list(z = matrix(1, const$I, const$J),",
" u = array(1, dim = c(const$I, const$J, const$K)),",
" r = array(stats::rnorm(const$I * const$J * const$K,",
" mean = 1, sd = 0.1),",
" dim = c(const$I, const$J, const$K)),")
if (cases$phi_shared[i]) {
ans <- c(ans, " alpha_shared = stats::rnorm(margs$M_phi_shared, sd = 0.1),")
}
if (cases$theta_shared[i]) {
ans <- c(ans, " beta_shared = stats::rnorm(margs$M_theta_shared, sd = 0.1),")
}
if (cases$psi_shared[i]) {
ans <- c(ans, " gamma_shared = stats::rnorm(margs$M_psi_shared, sd = 0.1),")
}
ans <- c(ans,
" spec_eff = matrix(stats::rnorm(const$I * margs$M, sd = 0.1),",
" const$I, margs$M),",
" Mu = stats::rnorm(margs$M, sd = 0.1),",
" sigma = stats::rnorm(margs$M, mean = 1, sd = 0.1))",
"}")
res <- inits_code(const, margs, dat)
expect_equal(res, ans)
# Tests for prior_ulim < 10 cases
dat$prior_ulim <- runif(1, min = 0, max = 10)
ans <- c("function() {",
" list(z = matrix(1, const$I, const$J),",
" u = array(1, dim = c(const$I, const$J, const$K)),",
" r = array(stats::rnorm(const$I * const$J * const$K,",
" mean = 1, sd = 0.1),",
" dim = c(const$I, const$J, const$K)),")
if (cases$phi_shared[i]) {
ans <- c(ans, " alpha_shared = stats::rnorm(margs$M_phi_shared, sd = 0.1),")
}
if (cases$theta_shared[i]) {
ans <- c(ans, " beta_shared = stats::rnorm(margs$M_theta_shared, sd = 0.1),")
}
if (cases$psi_shared[i]) {
ans <- c(ans, " gamma_shared = stats::rnorm(margs$M_psi_shared, sd = 0.1),")
}
ans <- c(ans,
" spec_eff = matrix(stats::rnorm(const$I * margs$M, sd = 0.1),",
" const$I, margs$M),",
" Mu = stats::rnorm(margs$M, sd = 0.1),")
ans <- c(ans,
sprintf(" sigma = stats::runif(margs$M, min = 0, max = %s))",
dat$prior_ulim))
ans <- c(ans, "}")
res <- inits_code(const, margs, dat)
expect_equal(res, ans)
}
})
### Tests for set_params_monitored() -------------------------------------------
test_that("Parameter list is correct for 144 available models", {
for (i in 1:nrow(cases)) {
ans <- c("Mu", "sigma", "rho", "alpha", "beta", "gamma")
if (cases$phi_shared[i]) {
ans <- c(ans, "alpha_shared")
}
if (cases$theta_shared[i]) {
ans <- c(ans, "beta_shared")
}
if (cases$psi_shared[i]) {
ans <- c(ans, "gamma_shared")
}
ans <- c(ans, c("phi", "theta", "psi", "z", "pi"))
res <- set_params_monitored(phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
expect_equal(res, ans)
}
})
### Tests for write_jags_model() -----------------------------------------------
test_that("JAGS code is correct for 144 available models", {
for (i in 1:nrow(cases)) {
ans <- readLines(system.file("jags",
"occumb_template1.jags",
package = "occumb"))
if (cases$phi[i] == "i") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i], 1)")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i, j], 1)")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" r[i, j, k] ~ dgamma(phi[i, j, k], 1)")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template2.jags",
package = "occumb")))
if (cases$theta[i] == "i") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i, j])")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" u[i, j, k] ~ dbern(z[i, j] * theta[i, j, k])")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template3.jags",
package = "occumb")))
if (cases$psi[i] == "i") {
ans <- c(ans,
" z[i, j] ~ dbern(psi[i])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" z[i, j] ~ dbern(psi[i, j])")
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template4.jags",
package = "occumb")))
if (cases$phi_shared[i]) {
if (cases$phi[i] == "i") {
ans <- c(ans,
" log(phi[i]) <- inprod(alpha[i, ], cov_phi[]) + inprod(alpha_shared[], cov_phi_shared[i, ])")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" log(phi[i, j]) <- inprod(alpha[i, ], cov_phi[j, ]) + inprod(alpha_shared[], cov_phi_shared[i, j, ])",
" }")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" log(phi[i, j, k]) <- inprod(alpha[i, ], cov_phi[j, k, ]) + inprod(alpha_shared[], cov_phi_shared[i, j, k, ])",
" }",
" }")
}
} else {
if (cases$phi[i] == "i") {
ans <- c(ans,
" log(phi[i]) <- inprod(alpha[i, ], cov_phi[])")
}
if (cases$phi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" log(phi[i, j]) <- inprod(alpha[i, ], cov_phi[j, ])",
" }")
}
if (cases$phi[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" log(phi[i, j, k]) <- inprod(alpha[i, ], cov_phi[j, k, ])",
" }",
" }")
}
}
if (cases$theta_shared[i]) {
if (cases$theta[i] == "i") {
ans <- c(ans,
" logit(theta[i]) <- inprod(beta[i, ], cov_theta[]) + inprod(beta_shared[], cov_theta_shared[i, ])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(theta[i, j]) <- inprod(beta[i, ], cov_theta[j, ]) + inprod(beta_shared[], cov_theta_shared[i, j, ])",
" }")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" logit(theta[i, j, k]) <- inprod(beta[i, ], cov_theta[j, k, ]) + inprod(beta_shared[], cov_theta_shared[i, j, k, ])",
" }",
" }")
}
} else {
if (cases$theta[i] == "i") {
ans <- c(ans,
" logit(theta[i]) <- inprod(beta[i, ], cov_theta[])")
}
if (cases$theta[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(theta[i, j]) <- inprod(beta[i, ], cov_theta[j, ])",
" }")
}
if (cases$theta[i] == "ijk") {
ans <- c(ans,
" for (j in 1:J) {",
" for (k in 1:K) {",
" logit(theta[i, j, k]) <- inprod(beta[i, ], cov_theta[j, k, ])",
" }",
" }")
}
}
if (cases$psi_shared[i]) {
if (cases$psi[i] == "i") {
ans <- c(ans,
" logit(psi[i]) <- inprod(gamma[i, ], cov_psi[]) + inprod(gamma_shared[], cov_psi_shared[i, ])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(psi[i, j]) <- inprod(gamma[i, ], cov_psi[j, ]) + inprod(gamma_shared[], cov_psi_shared[i, j, ])",
" }")
}
} else {
if (cases$psi[i] == "i") {
ans <- c(ans,
" logit(psi[i]) <- inprod(gamma[i, ], cov_psi[])")
}
if (cases$psi[i] == "ij") {
ans <- c(ans,
" for (j in 1:J) {",
" logit(psi[i, j]) <- inprod(gamma[i, ], cov_psi[j, ])",
" }")
}
}
ans <- c(ans, readLines(system.file("jags",
"occumb_template5.jags",
package = "occumb")))
if (cases$phi_shared[i]) {
ans <- c(ans,
" for (m in 1:M_phi_shared) {",
" alpha_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
if (cases$theta_shared[i]) {
ans <- c(ans,
" for (m in 1:M_theta_shared) {",
" beta_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
if (cases$psi_shared[i]) {
ans <- c(ans,
" for (m in 1:M_psi_shared) {",
" gamma_shared[m] ~ dnorm(0, prior_prec)",
" }")
}
ans <- c(ans, "}", "")
res <- write_jags_model(phi = cases$phi[i],
theta = cases$theta[i],
psi = cases$psi[i],
phi_shared = cases$phi_shared[i],
theta_shared = cases$theta_shared[i],
psi_shared = cases$psi_shared[i])
expect_equal(res, ans)
}
})
### Test for get_data() --------------------------------------------------------
test_that("Outputs are correct when proper variable names are given", {
I <- 2
J <- 2
K <- 2
y <- array(sample.int(I * J * K), dim = c(I, J, K))
spec_cov <- list(cov1 = rnorm(I))
site_cov <- list(cov2 = rnorm(J), cov3 = factor(1:J))
repl_cov <- list(cov4 = matrix(rnorm(J * K), J, K))
data <- occumbData(y = y,
spec_cov = spec_cov,
site_cov = site_cov,
repl_cov = repl_cov)
fit <- occumb(data = data,
n.chains = 1, n.burnin = 10, n.thin = 1, n.iter = 20,
verbose = FALSE)
expect_identical(get_data(fit, "y"), y)
expect_identical(get_data(fit, "spec_cov"), spec_cov)
expect_identical(get_data(fit, "site_cov"), site_cov)
expect_identical(get_data(fit, "repl_cov"), repl_cov)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.