added test for jags engine

poissonconsulting · Dec 14, 2023 · bdc40e8 · bdc40e8
1 parent a0106d5
commit bdc40e8
Showing 1 changed file with 132 additions and 0 deletions.
diff --git a/tests/testthat/test-zzz-analyse.R b/tests/testthat/test-zzz-analyse.R
@@ -532,3 +532,135 @@ random_effects = list(bSiteYear = c("Site", "YearFactor")))
   expect_identical(colnames(glance), c("n", "K", "nchains", "niters", "rhat_1", "rhat_2", "rhat_all", "converged"))
 })
 
+test_that("analyse full nimble notation vectorized embedded nested expression", {
+  set_analysis_mode("quick")
+
+  data <- embr::density99
+  data$YearFactor <- factor(data$Year)
+
+  model <- model(code = {
+
+  bIntercept ~ dnorm(0, sd = 5)
+  bYear ~ dnorm(0, sd = .5)
+
+  bHabitatQuality[1] <- 0
+  for(i in 2:nHabitatQuality) {
+    bHabitatQuality[i] ~ T(dnorm(0, sd = 5.), 0,)
+  }
+
+  log_sSiteYear ~ dlnorm(0, sd = 5)
+  log_sDensity ~ dt(0, sd = 5, 4.5)
+
+  log(sSiteYear) <- log_sSiteYear
+  log(sDensity) <- log_sDensity
+
+  for(i in 1:nSite) {
+    for(j in 1:nYearFactor) {
+      bSiteYear[i,j] ~ dnorm(0, sd = sSiteYear)
+    }
+  }
+
+  for(i in 1:length(Density)) {
+    eDensity[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]
+    Density[i] ~ dlnorm(eDensity[i], sd = sDensity)
+  }
+},
+new_expr = {
+  for(i in 1:length(Density)) {
+    fit[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]
+    log(prediction[i]) <- fit[i]
+    residual[i] <- res_lnorm(Density[i], fit[i], exp(log_sDensity))
+  }
+},
+new_expr_vec = TRUE,
+select_data = list("Year+" = numeric(), YearFactor = factor(),
+                   Site = factor(), Density = numeric(),
+                   HabitatQuality = factor()),
+fixed = "^(b|l)", derived = "eDensity",
+random_effects = list(bSiteYear = c("Site", "YearFactor")))
+
+  expect_output(expect_warning(analysis <- analyse(model, data = data, engine = "jags")))
+
+  expect_equal(as.data.frame(data_set(analysis)), data)
+  data2 <- data_set(analysis, marginalize_random_effects = TRUE)
+  expect_true(all(as.integer(data2$Site) == 1L))
+  expect_true(all(as.integer(data2$YearFactor) == 1L))
+
+  expect_identical(class(analysis), c("jmb_analysis", "mb_analysis"))
+
+  expect_identical(niters(analysis), 10L)
+  expect_identical(nchains(analysis), 2L)
+  expect_identical(nsims(analysis), 20L)
+  expect_identical(ngens(analysis), 40L)
+
+  expect_output(analysis <- reanalyse(analysis))
+
+  expect_identical(niters(analysis), 10L)
+  expect_identical(ngens(analysis), 40L)
+
+  expect_identical(pars(analysis, "fixed"), sort(c("bHabitatQuality", "bIntercept", "bYear", "log_sDensity", "log_sSiteYear")))
+  expect_identical(pars(analysis, "random"), "bSiteYear")
+  expect_identical(pars(analysis, "derived"), "eDensity")
+  expect_identical(pars(analysis, "primary"),
+                   c("bHabitatQuality", "bIntercept", "bSiteYear", "bYear", "log_sDensity", "log_sSiteYear"))
+  expect_identical(pars(analysis),
+                   c("bHabitatQuality", "bIntercept", "bSiteYear", "bYear", "eDensity", "log_sDensity", "log_sSiteYear"))
+
+  expect_is(as.mcmcr(analysis), "mcmcr")
+
+  glance <- glance(analysis)
+  expect_is(glance, "tbl")
+  expect_identical(colnames(glance), c("n", "K", "nchains", "niters", "nthin", "ess", "rhat", "converged"))
+  expect_identical(glance$n, 300L)
+  expect_identical(glance$nthin, 1L)
+  expect_identical(glance$K, 5L)
+
+  derived <- coef(analysis, param_type = "derived", simplify = TRUE)
+  expect_identical(colnames(derived), c("term", "estimate", "lower", "upper", "svalue"))
+  expect_identical(nrow(derived), 300L)
+
+  coef <- coef(analysis, simplify = TRUE)
+
+  expect_is(coef, "tbl")
+  expect_identical(colnames(coef), c("term", "estimate", "lower", "upper", "svalue"))
+
+  expect_identical(coef$term, as.term(c("bHabitatQuality[1]", "bHabitatQuality[2]",
+                                        "bIntercept", "bYear",
+                                        "log_sDensity", "log_sSiteYear")))
+
+  expect_identical(nrow(coef(analysis, "primary", simplify = TRUE)), 66L)
+  expect_identical(nrow(coef(analysis, "all", simplify = TRUE)), 366L)
+
+  tidy <- tidy(analysis)
+  expect_identical(colnames(tidy), c("term", "estimate", "lower", "upper", "esr", "rhat"))
+
+  year <- predict(analysis, new_data = "Year")
+
+  ppc <- posterior_predictive_check(analysis)
+
+  expect_is(ppc, "tbl_df")
+  expect_identical(colnames(ppc), c("moment", "observed", "median", "lower", "upper", "svalue"))
+  expect_identical(ppc$moment, structure(1:5, .Label = c("zeros", "mean", "variance", "skewness",
+                                                         "kurtosis"), class = "factor"))
+  expect_is(year, "tbl")
+  expect_identical(colnames(year), c("Site", "HabitatQuality", "Year", "Visit",
+                                     "Density", "YearFactor",
+                                     "estimate", "lower", "upper", "svalue"))
+  expect_true(all(year$lower < year$estimate))
+  expect_false(is.unsorted(year$estimate))
+
+  dd <- mcmc_derive_data(analysis, new_data = c("Site", "Year"), ref_data = TRUE)
+  expect_true(mcmcdata::is.mcmc_data(dd))
+
+  expect_warning(sensitivity <- sd_priors_by(analysis, by = 10, glance = FALSE), "incomplete adaptation")
+
+  analyses <- analyses(analysis, sensitivity)
+
+  glance <- glance(analyses, bound = TRUE)
+
+  expect_is(glance, "tbl")
+  expect_identical(nrow(glance), 1L)
+  expect_identical(colnames(glance), c("n", "K", "nchains", "niters", "rhat_1", "rhat_2", "rhat_all", "converged"))
+})
+
+