Skip to content

Commit

Permalink
Merge pull request #22 from poissonconsulting/f-style
Browse files Browse the repository at this point in the history
  • Loading branch information
joethorley authored Aug 9, 2024
2 parents bdc40e8 + cd4bdfe commit b7a996a
Show file tree
Hide file tree
Showing 8 changed files with 285 additions and 195 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -33,5 +33,5 @@ Remotes:
poissonconsulting/embr
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.2.3
RoxygenNote: 7.3.2
VignetteBuilder: knitr
43 changes: 24 additions & 19 deletions R/analyse.R
Original file line number Diff line number Diff line change
@@ -1,30 +1,32 @@
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(
jags_model <- rjags::jags.model(loaded, data, inits = inits, n.adapt = 0, quiet = quiet)
)

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")

update(jags_model, n.iter = floor(niters / 2), progress.bar = "none")

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()

Expand All @@ -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")
Expand All @@ -62,4 +68,3 @@ analyse1.jmb_model <- function(model, data, loaded, nchains, niters, nthin, quie

obj
}

4 changes: 3 additions & 1 deletion R/inits.R
Original file line number Diff line number Diff line change
@@ -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) {
Expand Down
40 changes: 24 additions & 16 deletions R/sd-priors-by.R
Original file line number Diff line number Diff line change
@@ -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)
}
Expand All @@ -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)
}
7 changes: 4 additions & 3 deletions README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -107,7 +108,7 @@ devtools::install_github("poissonconsulting/jmbr")
## Citation

```{r, comment="", echo=FALSE}
citation(package = 'jmbr')
citation(package = "jmbr")
```

## Contribution
Expand Down
13 changes: 8 additions & 5 deletions tests/testthat/test-code.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
test_that("analyse", {

template <- "model{
bIntercept ~ dlogis(0, 5^-2)
Expand Down Expand Up @@ -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}"
)
})
18 changes: 11 additions & 7 deletions tests/testthat/test-model.R
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
Loading

0 comments on commit b7a996a

Please sign in to comment.