Skip to content

Commit

Permalink
Issue #266: Include intercept when no model provided (#370)
Browse files Browse the repository at this point in the history
* Scope notes on forcing a ~ 1 when noting provided

* Update globals.R

* Run document()

* The structure for the .make_intercepts_explicit helper function insertion

* Set-up test and documentation for helper function

* Return formula in placeholder

* First draft of "add ~ 1" functionality

* Trying to test by comparison to explicitly created formula

* Missing bracket

* Don't change the formula for parameters set to be a constant

* Simplify prior structure

* Set environments to NULL such that formulas are equal

* Fix to bug in prior checking test
  • Loading branch information
athowes authored Oct 10, 2024
1 parent 94812c2 commit 4150efa
Show file tree
Hide file tree
Showing 6 changed files with 75 additions and 10 deletions.
2 changes: 1 addition & 1 deletion R/latent_individual.R
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ epidist_formula.epidist_latent_individual <- function(data, family, formula,
...) {
epidist_validate(data)
formula <- brms:::validate_formula(formula, data = data, family = family)

formula <- .make_intercepts_explicit(formula)
formula <- stats::update(
formula, delay | vreal(relative_obs_time, pwindow, swindow) ~ .
)
Expand Down
11 changes: 3 additions & 8 deletions R/prior.R
Original file line number Diff line number Diff line change
Expand Up @@ -90,18 +90,13 @@ epidist_family_prior.default <- function(family, formula, ...) {
#' @export
epidist_family_prior.lognormal <- function(family, formula, ...) {
prior <- prior("normal(1, 1)", class = "Intercept")
if ("sigma" %in% names(formula$pforms)) {
# Case with a model on sigma
sigma_prior <- prior(
"normal(-0.7, 0.4)", class = "Intercept", dpar = "sigma"
)
} else if ("sigma" %in% names(formula$pfix)) {
if ("sigma" %in% names(formula$pfix)) {
# Case with sigma fixed to a constant
sigma_prior <- NULL
} else {
# Case with no model on sigma
# Case with a model on sigma
sigma_prior <- prior(
"lognormal(-0.7, 0.4)", class = "sigma", lb = 0, ub = "NA"
"normal(-0.7, 0.4)", class = "Intercept", dpar = "sigma"
)
}
prior <- prior + sigma_prior
Expand Down
20 changes: 20 additions & 0 deletions R/utils.R
Original file line number Diff line number Diff line change
Expand Up @@ -102,3 +102,23 @@
family$other_bounds <- other_bounds
return(family)
}

#' Include implicit intercepts in `brms` formula as explicit
#'
#' This function detects the distributional parameters in a `brms` formula
#' object, and alters to formula to include explicit intercept parameters for
#' them i.e. `~ 1`.
#'
#' @param formula ...
#' @keywords internal
.make_intercepts_explicit <- function(formula) {
other_dpars <- setdiff(formula$family$dpars, "mu")
fixed_dpars <- names(formula$pfix)
formula_dpars <- names(formula$pforms)
replace_dpars <- setdiff(other_dpars, c(fixed_dpars, formula_dpars))
for (dpar in replace_dpars) {
new_formula <- as.formula(paste0(dpar, " ~ 1"))
formula$pforms[[dpar]] <- new_formula
}
return(formula)
}
17 changes: 17 additions & 0 deletions man/dot-make_intercepts_explicit.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion tests/testthat/test-int-latent_individual.R
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ test_that("epidist.epidist_latent_individual samples from the prior according to
param1 <- extract_normal_parameters_brms(epidist_prior[1, ])
param2 <- extract_normal_parameters_brms(epidist_prior[2, ])
samples1 <- rnorm(1000, mean = param1$mean, sd = param1$sd)
samples2 <- exp(rnorm(1000, mean = param2$mean, sd = param2$sd))
samples2 <- rnorm(1000, mean = param2$mean, sd = param2$sd)
# suppressWarnings here used to prevent warnings about ties
ks1 <- suppressWarnings(stats::ks.test(pred$mu, samples1))
ks2 <- suppressWarnings(stats::ks.test(pred$sigma, samples2))
Expand Down
33 changes: 33 additions & 0 deletions tests/testthat/test-utils.R
Original file line number Diff line number Diff line change
Expand Up @@ -34,3 +34,36 @@ test_that(".add_dpar_info works as expected for the lognormal and gamma families
expect_equal(gamma_extra$other_links, NULL)
expect_equal(gamma_extra$other_bounds, list(list("lb" = "0", ub = "")))
})

test_that(".make_intercepts_explicit creates a formula which is the same as if it had been explicitly created", { # nolint: line_length_linter.
prep_obs <- as_latent_individual(sim_obs)
epidist_family <- epidist_family(prep_obs, family = "lognormal")
formula <- brms:::validate_formula(
formula = brms::bf(mu ~ 1),
data = prep_obs,
family = epidist_family
)
formula <- .make_intercepts_explicit(formula)
formula_explicit <- brms:::validate_formula(
formula = brms::bf(mu ~ 1, sigma ~ 1),
data = prep_obs,
family = epidist_family
)
attr(formula$pforms$sigma, ".Environment") <- NULL
attr(formula_explicit$pforms$sigma, ".Environment") <- NULL
expect_equal(formula, formula_explicit)
})

test_that(".make_intercepts_explicit does not add an intercept if the distributional parameter is set to be fixed", { # nolint: line_length_linter.
prep_obs <- as_latent_individual(sim_obs)
epidist_family <- epidist_family(prep_obs, family = "lognormal")
formula <- brms:::validate_formula(
formula = brms::bf(mu ~ 1, sigma = 1),
data = prep_obs,
family = epidist_family
)
formula_updated <- .make_intercepts_explicit(formula)
attr(formula$pforms$sigma, ".Environment") <- NULL
attr(formula_updated$pforms$sigma, ".Environment") <- NULL
expect_equal(formula, formula_updated)
})

0 comments on commit 4150efa

Please sign in to comment.