From cd4bdfe4f15540c2c72b739a54a3fb2c8efea390 Mon Sep 17 00:00:00 2001 From: Joe Thorley Date: Thu, 8 Aug 2024 17:02:26 -0700 Subject: [PATCH] style code --- DESCRIPTION | 2 +- R/analyse.R | 43 ++-- R/inits.R | 4 +- R/sd-priors-by.R | 40 ++-- README.Rmd | 7 +- tests/testthat/test-code.R | 13 +- tests/testthat/test-model.R | 18 +- tests/testthat/test-zzz-analyse.R | 353 ++++++++++++++++++------------ 8 files changed, 285 insertions(+), 195 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index eb4c8ea..91e775d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -33,5 +33,5 @@ Remotes: poissonconsulting/embr Encoding: UTF-8 LazyData: true -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.2 VignetteBuilder: knitr diff --git a/R/analyse.R b/R/analyse.R index 7ad5d31..fde6672 100644 --- a/R/analyse.R +++ b/R/analyse.R @@ -1,5 +1,4 @@ jmb_analyse_chain <- function(inits, loaded, data, monitor, niters, ngens, nthin, quiet) { - capture_output <- if (quiet) function(x) suppressWarnings(capture.output(x)) else eval capture_output( @@ -7,8 +6,10 @@ jmb_analyse_chain <- function(inits, loaded, data, monitor, niters, ngens, nthin ) niters <- niters * nthin - adapted <- rjags::adapt(jags_model, n.iter = floor(niters / 2), - progress.bar = "none", end.adaptation = TRUE) + adapted <- rjags::adapt(jags_model, + n.iter = floor(niters / 2), + progress.bar = "none", end.adaptation = TRUE + ) if (!adapted) warning("incomplete adaptation") @@ -16,15 +17,16 @@ jmb_analyse_chain <- function(inits, loaded, data, monitor, niters, ngens, nthin monitor <- monitor[monitor %in% stats::variable.names(jags_model)] - jags_samples <- rjags::jags.samples(model = jags_model, variable.names = monitor, - n.iter = niters, thin = nthin, progress.bar = "none") + jags_samples <- rjags::jags.samples( + model = jags_model, variable.names = monitor, + n.iter = niters, thin = nthin, progress.bar = "none" + ) list(jags_model = jags_model, jags_samples = jags_samples) } #' @export analyse1.jmb_model <- function(model, data, loaded, nchains, niters, nthin, quiet, glance, parallel, ...) { - timer <- timer::Timer$new() timer$start() @@ -37,23 +39,27 @@ analyse1.jmb_model <- function(model, data, loaded, nchains, niters, nthin, quie monitor <- embr::monitor(model) monitor <- monitor[!monitor %in% names(data)] - jags_chains <- llply(inits, .fun = jmb_analyse_chain, - .parallel = parallel, - loaded = loaded, - data = data, - monitor = monitor, - niters = niters, - nthin = nthin, - quiet = quiet) + jags_chains <- llply(inits, + .fun = jmb_analyse_chain, + .parallel = parallel, + loaded = loaded, + data = data, + monitor = monitor, + niters = niters, + nthin = nthin, + quiet = quiet + ) mcmc <- llply(jags_chains, function(x) x$jags_samples) mcmc <- lapply(mcmc, function(x) mcmcr::as.mcmcr(lapply(x, mcmcr::as.mcmcarray))) mcmc <- Reduce(mcmcr::bind_chains, mcmc) - obj <- c(obj, inits = list(inits), - jags_chains = list(jags_chains), - mcmcr = list(mcmc), - nthin = nthin) + obj <- c(obj, + inits = list(inits), + jags_chains = list(jags_chains), + mcmcr = list(mcmc), + nthin = nthin + ) obj$duration <- timer$elapsed() class(obj) <- c("jmb_analysis", "mb_analysis") @@ -62,4 +68,3 @@ analyse1.jmb_model <- function(model, data, loaded, nchains, niters, nthin, quie obj } - diff --git a/R/inits.R b/R/inits.R index 39be667..3b0e625 100644 --- a/R/inits.R +++ b/R/inits.R @@ -1,7 +1,9 @@ inits <- function(data, gen_inits, nchains) { rngs <- rjags::parallel.seeds("base::BaseRNG", nchains) - if (identical(gen_inits(data), list())) return(rngs) + if (identical(gen_inits(data), list())) { + return(rngs) + } inits <- list() for (i in 1:nchains) { diff --git a/R/sd-priors-by.R b/R/sd-priors-by.R index ae50c94..8e97524 100644 --- a/R/sd-priors-by.R +++ b/R/sd-priors-by.R @@ -1,16 +1,18 @@ #' @export sd_priors_by.jmb_code <- function( - x, by = 10, distributions = c("normal", "lognormal", "t"), ...) { + x, by = 10, distributions = c("normal", "lognormal", "t"), ...) { chk_number(by) chk_range(by, c(0.001, 1000)) chk_unused(...) chk_s3_class(distributions, "character") chk_unique(distributions) - chk_subset(distributions, c("laplace", "logistic", "lognormal", - "normal", "t", "nt")) + chk_subset(distributions, c( + "laplace", "logistic", "lognormal", + "normal", "t", "nt" + )) - if(!length(distributions)) { + if (!length(distributions)) { wrn("No prior distributions included.") return(x) } @@ -20,18 +22,24 @@ sd_priors_by.jmb_code <- function( pattern2 <- "\\s*[(][^,)]+,\\s*)((\\d+[.]{0,1}\\d*)|(\\d*[.]{0,1}\\d+))(\\s*\\^\\s*-\\s*2\\s*)([)])" pattern3 <- "\\s*[(][^,)]+,\\s*)((\\d+[.]{0,1}\\d*)|(\\d*[.]{0,1}\\d+))(\\s*\\^\\s*-\\s*2\\s*)(,[^,)]+[)])" replacement <- paste0("\\1(\\2 * ", by, ")^-2\\6") - if("laplace" %in% distributions) - x <- gsub(paste0("(~\\s*ddexp", pattern2), replacement, x) - if("logistic" %in% distributions) - x <- gsub(paste0("(~\\s*dlogis", pattern2), replacement, x) - if("lognormal" %in% distributions) - x <- gsub(paste0("(~\\s*dlnorm", pattern2), replacement, x) - if("normal" %in% distributions) - x <- gsub(paste0("(~\\s*dnorm", pattern2), replacement, x) - if("t" %in% distributions) - x <- gsub(paste0("(~\\s*dt", pattern3), replacement, x) - if("nt" %in% distributions) - x <- gsub(paste0("(~\\s*dnt", pattern3), replacement, x) + if ("laplace" %in% distributions) { + x <- gsub(paste0("(~\\s*ddexp", pattern2), replacement, x) + } + if ("logistic" %in% distributions) { + x <- gsub(paste0("(~\\s*dlogis", pattern2), replacement, x) + } + if ("lognormal" %in% distributions) { + x <- gsub(paste0("(~\\s*dlnorm", pattern2), replacement, x) + } + if ("normal" %in% distributions) { + x <- gsub(paste0("(~\\s*dnorm", pattern2), replacement, x) + } + if ("t" %in% distributions) { + x <- gsub(paste0("(~\\s*dt", pattern3), replacement, x) + } + if ("nt" %in% distributions) { + x <- gsub(paste0("(~\\s*dnt", pattern3), replacement, x) + } mb_code(x) } diff --git a/README.Rmd b/README.Rmd index af5b297..8d6bd96 100644 --- a/README.Rmd +++ b/README.Rmd @@ -62,10 +62,11 @@ for (i in 1:length(Pairs)) { }") # define data types and center year -model <- update_model(model, +model <- update_model(model, select_data = list("Pairs" = integer(), "Year*" = integer(), Annual = factor()), derived = "sAnnual", - random_effects = list(bAnnual = "Annual")) + random_effects = list(bAnnual = "Annual") +) data <- bauw::peregrine data$Annual <- factor(data$Year) @@ -107,7 +108,7 @@ devtools::install_github("poissonconsulting/jmbr") ## Citation ```{r, comment="", echo=FALSE} -citation(package = 'jmbr') +citation(package = "jmbr") ``` ## Contribution diff --git a/tests/testthat/test-code.R b/tests/testthat/test-code.R index 3350521..52e4e31 100644 --- a/tests/testthat/test-code.R +++ b/tests/testthat/test-code.R @@ -1,5 +1,4 @@ test_that("analyse", { - template <- "model{ bIntercept ~ dlogis(0, 5^-2) @@ -42,12 +41,16 @@ test_that("analyse", { code10 <- sd_priors_by(code, 10, distributions = c("logistic", "normal", "lognormal", "t")) expect_true(is.jmb_code(code10)) expect_identical(pars(code10), pars(code)) - expect_identical(as.character(code10), - "model{\n\n bIntercept ~ dlogis(0, (5 * 10)^-2)\n bYear ~ dnorm(0, (.5 * 10)^-2)\n\n bHabitatQuality[1] <- 0\n for(i in 2:nHabitatQuality) {\n bHabitatQuality[i] ~ dnorm(0, (5. * 10)^-2) T(0,)\n }\n\n log_sSiteYear ~ dlnorm(0, (5 * 10)^-2)\n log_sDensity ~ dt(0, (5 * 10)^-2, 4.5)\n\n log(sSiteYear) <- log_sSiteYear\n log(sDensity) <- log_sDensity\n\n for(i in 1:nSite) {\n for(j in 1:nYearFactor) {\n bSiteYear[i,j] ~ dnorm(0, sSiteYear^-2)\n }\n }\n\n for(i in 1:length(Density)) {\n eDensity[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]\n Density[i] ~ dlnorm(eDensity[i], sDensity^-2)\n }\n}") + expect_identical( + as.character(code10), + "model{\n\n bIntercept ~ dlogis(0, (5 * 10)^-2)\n bYear ~ dnorm(0, (.5 * 10)^-2)\n\n bHabitatQuality[1] <- 0\n for(i in 2:nHabitatQuality) {\n bHabitatQuality[i] ~ dnorm(0, (5. * 10)^-2) T(0,)\n }\n\n log_sSiteYear ~ dlnorm(0, (5 * 10)^-2)\n log_sDensity ~ dt(0, (5 * 10)^-2, 4.5)\n\n log(sSiteYear) <- log_sSiteYear\n log(sDensity) <- log_sDensity\n\n for(i in 1:nSite) {\n for(j in 1:nYearFactor) {\n bSiteYear[i,j] ~ dnorm(0, sSiteYear^-2)\n }\n }\n\n for(i in 1:length(Density)) {\n eDensity[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]\n Density[i] ~ dlnorm(eDensity[i], sDensity^-2)\n }\n}" + ) coder <- sd_priors_by(code, 10, distributions = "normal") expect_true(is.jmb_code(coder)) expect_identical(pars(coder), pars(code)) - expect_identical(as.character(coder), - "model{\n\n bIntercept ~ dlogis(0, 5^-2)\n bYear ~ dnorm(0, (.5 * 10)^-2)\n\n bHabitatQuality[1] <- 0\n for(i in 2:nHabitatQuality) {\n bHabitatQuality[i] ~ dnorm(0, (5. * 10)^-2) T(0,)\n }\n\n log_sSiteYear ~ dlnorm(0, 5^-2)\n log_sDensity ~ dt(0, 5^-2, 4.5)\n\n log(sSiteYear) <- log_sSiteYear\n log(sDensity) <- log_sDensity\n\n for(i in 1:nSite) {\n for(j in 1:nYearFactor) {\n bSiteYear[i,j] ~ dnorm(0, sSiteYear^-2)\n }\n }\n\n for(i in 1:length(Density)) {\n eDensity[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]\n Density[i] ~ dlnorm(eDensity[i], sDensity^-2)\n }\n}") + expect_identical( + as.character(coder), + "model{\n\n bIntercept ~ dlogis(0, 5^-2)\n bYear ~ dnorm(0, (.5 * 10)^-2)\n\n bHabitatQuality[1] <- 0\n for(i in 2:nHabitatQuality) {\n bHabitatQuality[i] ~ dnorm(0, (5. * 10)^-2) T(0,)\n }\n\n log_sSiteYear ~ dlnorm(0, 5^-2)\n log_sDensity ~ dt(0, 5^-2, 4.5)\n\n log(sSiteYear) <- log_sSiteYear\n log(sDensity) <- log_sDensity\n\n for(i in 1:nSite) {\n for(j in 1:nYearFactor) {\n bSiteYear[i,j] ~ dnorm(0, sSiteYear^-2)\n }\n }\n\n for(i in 1:length(Density)) {\n eDensity[i] <- bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]\n Density[i] ~ dlnorm(eDensity[i], sDensity^-2)\n }\n}" + ) }) diff --git a/tests/testthat/test-model.R b/tests/testthat/test-model.R index df507d7..df48cb3 100644 --- a/tests/testthat/test-model.R +++ b/tests/testthat/test-model.R @@ -32,13 +32,17 @@ test_that("analyse", { prediction[i] <- exp(bIntercept + bYear * Year[i] + bHabitatQuality[HabitatQuality[i]] + bSiteYear[Site[i], YearFactor[i]]) } " - model <- model(code = template, - 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")), - new_expr = new_expr) + model <- model( + code = template, + 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")), + new_expr = new_expr + ) expect_identical(class(model), c("jmb_model", "mb_model")) expect_true(is.jmb_model(model)) diff --git a/tests/testthat/test-zzz-analyse.R b/tests/testthat/test-zzz-analyse.R index 712990b..31724a3 100644 --- a/tests/testthat/test-zzz-analyse.R +++ b/tests/testthat/test-zzz-analyse.R @@ -39,13 +39,17 @@ test_that("analyse character vector", { residual[i] <- res_lnorm(Density[i], fit[i], exp(log_sDensity)) }" - model <- model(code = jags_template, - 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")), - new_expr = new_expr) + model <- model( + code = jags_template, + 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")), + new_expr = new_expr + ) expect_output(expect_warning(analysis <- analyse(model, data = data))) @@ -78,10 +82,14 @@ test_that("analyse character vector", { 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_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") @@ -101,9 +109,11 @@ test_that("analyse character vector", { 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(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) @@ -117,12 +127,16 @@ test_that("analyse character vector", { 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_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_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)) @@ -171,14 +185,18 @@ test_that("analyse vectorized stand alone character", { residual[i] <- res_lnorm(Density[i], fit[i], exp(log_sDensity)) }" - model <- model(code = jags_template, - 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")), - new_expr = new_expr, - new_expr_vec = TRUE) + model <- model( + code = jags_template, + 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")), + new_expr = new_expr, + new_expr_vec = TRUE + ) expect_output(expect_warning(analysis <- analyse(model, data = data))) @@ -211,10 +229,14 @@ test_that("analyse vectorized stand alone character", { 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_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") @@ -234,9 +256,11 @@ test_that("analyse vectorized stand alone character", { 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(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) @@ -250,12 +274,16 @@ test_that("analyse vectorized stand alone character", { 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_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_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)) @@ -269,7 +297,8 @@ test_that("analyse vectorized embedded expression", { data <- embr::density99 data$YearFactor <- factor(data$Year) - model <- model(code = "model{ + model <- model( + code = "model{ bIntercept ~ dnorm(0, 5^-2) bYear ~ dnorm(0, .5^-2) # bYear2 ~ dnorm(0, .5^-2) @@ -296,17 +325,20 @@ test_that("analyse vectorized embedded expression", { Density[i] ~ dlnorm(eDensity[i], sDensity^-2) } }", -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"))) + 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))) @@ -339,10 +371,14 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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") @@ -362,9 +398,11 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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(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) @@ -378,12 +416,16 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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_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)) @@ -407,7 +449,8 @@ test_that("analyse vectorized embedded nested expression", { data <- embr::density99 data$YearFactor <- factor(data$Year) - model <- model(code = "model{ + model <- model( + code = "model{ bIntercept ~ dnorm(0, 5^-2) bYear ~ dnorm(0, .5^-2) # bYear2 ~ dnorm(0, .5^-2) @@ -434,19 +477,22 @@ test_that("analyse vectorized embedded nested expression", { Density[i] ~ dlnorm(eDensity[i], sDensity^-2) } }", -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"))) + 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))) @@ -470,10 +516,14 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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") @@ -493,9 +543,11 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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(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) @@ -509,12 +561,16 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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_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)) @@ -538,46 +594,49 @@ test_that("analyse full nimble notation vectorized embedded nested expression", 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"))) + 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"))) @@ -601,10 +660,14 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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") @@ -624,9 +687,11 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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(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) @@ -640,12 +705,16 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) 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_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_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)) @@ -662,5 +731,3 @@ random_effects = list(bSiteYear = c("Site", "YearFactor"))) expect_identical(nrow(glance), 1L) expect_identical(colnames(glance), c("n", "K", "nchains", "niters", "rhat_1", "rhat_2", "rhat_all", "converged")) }) - -