From bdc40e853c25451249595876ef2386f962485c9c Mon Sep 17 00:00:00 2001 From: Joe Thorley Date: Thu, 14 Dec 2023 06:32:21 -0600 Subject: [PATCH] added test for jags engine --- tests/testthat/test-zzz-analyse.R | 132 ++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) diff --git a/tests/testthat/test-zzz-analyse.R b/tests/testthat/test-zzz-analyse.R index 549a8f6..712990b 100644 --- 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")) +}) + +