diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a33aea8cc..d2ad27ffc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -21,7 +21,6 @@ repos: - id: check-added-large-files args: ['--maxkb=200'] - id: end-of-file-fixer - exclude: ^tests/testthat/_snaps - repo: local hooks: - id: forbid-to-commit diff --git a/NAMESPACE b/NAMESPACE index a03b623c0..93ba3ff12 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,13 +13,13 @@ S3method(assert_epidist,epidist_linelist_data) S3method(assert_epidist,epidist_naive_model) S3method(epidist_family_model,default) S3method(epidist_family_model,epidist_latent_model) +S3method(epidist_family_param,default) S3method(epidist_family_prior,default) S3method(epidist_family_prior,lognormal) -S3method(epidist_family_reparam,default) -S3method(epidist_family_reparam,gamma) S3method(epidist_formula_model,default) S3method(epidist_formula_model,epidist_latent_model) S3method(epidist_model_prior,default) +S3method(epidist_model_prior,epidist_latent_model) S3method(epidist_stancode,default) S3method(epidist_stancode,epidist_latent_model) export(Gamma) @@ -33,8 +33,8 @@ export(epidist) export(epidist_diagnostics) export(epidist_family) export(epidist_family_model) +export(epidist_family_param) export(epidist_family_prior) -export(epidist_family_reparam) export(epidist_formula) export(epidist_formula_model) export(epidist_gen_posterior_epred) @@ -57,9 +57,12 @@ export(simulate_secondary) export(simulate_uniform_cases) export(weibull) import(ggplot2) +importFrom(brms,as.brmsprior) importFrom(brms,bf) importFrom(brms,lognormal) +importFrom(brms,make_stancode) importFrom(brms,prior) +importFrom(brms,set_prior) importFrom(brms,stanvar) importFrom(brms,weibull) importFrom(checkmate,assert_class) @@ -67,6 +70,7 @@ importFrom(checkmate,assert_data_frame) importFrom(checkmate,assert_date) importFrom(checkmate,assert_factor) importFrom(checkmate,assert_integer) +importFrom(checkmate,assert_integerish) importFrom(checkmate,assert_names) importFrom(checkmate,assert_numeric) importFrom(checkmate,assert_true) @@ -75,12 +79,14 @@ importFrom(cli,cli_alert_info) importFrom(cli,cli_inform) importFrom(cli,cli_warn) importFrom(dplyr,bind_cols) +importFrom(dplyr,bind_rows) importFrom(dplyr,filter) importFrom(dplyr,full_join) importFrom(dplyr,mutate) importFrom(dplyr,select) importFrom(lubridate,days) importFrom(lubridate,is.timepoint) +importFrom(purrr,map_dbl) importFrom(stats,Gamma) importFrom(stats,as.formula) importFrom(stats,setNames) diff --git a/NEWS.md b/NEWS.md index 2153d2507..63d0ca7c8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ Development version of `epidist`. ## Package - Remove the default method for `epidist()`. See #473. +- Added `enforce_presence` argument to `epidist_prior()` to allow for priors to be + specified if they do not match existing parameters. See #474. +- Added a `merge` argument to `epidist_prior()` to allow for not merging user and package priors. See #474. +- Added user settable primary event priors to the latent model. See #474. +- Added a marginalised likelihood to the latent model. See #474. +- Generalised the stan reparametrisation feature to work across all distributions without manual specification by generating stan code with `brms` and then extracting the reparameterisation. See #474. ## Documentation diff --git a/R/epidist.R b/R/epidist.R index 7241f570d..78db416a8 100644 --- a/R/epidist.R +++ b/R/epidist.R @@ -14,7 +14,13 @@ #' reexported as part of `epidist`. #' @param prior One or more `brmsprior` objects created by [brms::set_prior()] #' or related functions. These priors are passed to [epidist_prior()] in the -#' `prior` argument. +#' `prior` argument. Some models have default priors that are automatically +#' added (see [epidist_model_prior()]). These can be merged with user-provided +#' priors using the `merge_priors` argument. +#' @param merge_priors If `TRUE` then merge user priors with default priors, if +#' `FALSE` only use user priors. Defaults to `TRUE`. This may be useful if +#' the built in approaches for merging priors are not flexible enough for a +#' particular use case. #' @param fn The internal function to be called. By default this is #' [brms::brm()] which performs inference for the specified model. Other options #' are [brms::make_stancode()] which returns the Stan code for the specified @@ -25,6 +31,7 @@ #' @export epidist <- function(data, formula = mu ~ 1, family = lognormal(), prior = NULL, + merge_priors = TRUE, fn = brms::brm, ...) { assert_epidist(data) epidist_family <- epidist_family(data, family) @@ -32,7 +39,8 @@ epidist <- function(data, formula = mu ~ 1, data = data, family = epidist_family, formula = formula ) epidist_prior <- epidist_prior( - data = data, family = epidist_family, formula = epidist_formula, prior + data = data, family = epidist_family, formula = epidist_formula, prior, + merge = merge_priors ) epidist_stancode <- epidist_stancode( data = data, family = epidist_family, formula = epidist_formula diff --git a/R/family.R b/R/family.R index 9512b9019..583100f0b 100644 --- a/R/family.R +++ b/R/family.R @@ -15,7 +15,7 @@ epidist_family <- function(data, family = lognormal(), ...) { family <- .add_dpar_info(family) custom_family <- epidist_family_model(data, family, ...) class(custom_family) <- c(family$family, class(custom_family)) - custom_family <- epidist_family_reparam(custom_family) + custom_family <- epidist_family_param(custom_family) return(custom_family) } @@ -43,29 +43,73 @@ epidist_family_model.default <- function(data, family, ...) { #' Reparameterise an `epidist` family to align `brms` and Stan #' #' @inheritParams epidist_family -#' @rdname epidist_family_reparam +#' @rdname epidist_family_param #' @family family #' @export -epidist_family_reparam <- function(family, ...) { - UseMethod("epidist_family_reparam") +epidist_family_param <- function(family, ...) { + UseMethod("epidist_family_param") } #' Default method for families which do not require a reparameterisation #' -#' @inheritParams epidist_family -#' @family family -#' @export -epidist_family_reparam.default <- function(family, ...) { - family$reparam <- family$dpars - return(family) -} - -#' Reparameterisation for the gamma family +#' This function extracts the Stan parameterisation for a given brms family by +#' creating a dummy model and parsing its Stan code. It looks for the log +#' probability density function (lpdf) call in the Stan code and extracts the +#' parameter order used by Stan. This is needed because brms and Stan may use +#' different parameter orderings for the same distribution. +#' +#' @param family A brms family object containing at minimum a `family` element +#' specifying the distribution family name +#' @param ... Additional arguments passed to methods (not used) +#' +#' @details +#' The function works by: +#' 1. Creating a minimal dummy model using the specified family +#' 2. Extracting the Stan code for this model +#' 3. Finding the lpdf function call for the family +#' 4. Parsing out the parameter ordering used in Stan +#' 5. Adding this as the `param` element to the family object +#' +#' @return The input family object with an additional `param` element containing +#' the Stan parameter ordering as a string #' -#' @inheritParams epidist_family #' @family family +#' @importFrom brms make_stancode +#' @importFrom cli cli_abort #' @export -epidist_family_reparam.gamma <- function(family, ...) { - family$reparam <- c("shape", "shape ./ mu") # nolint +epidist_family_param.default <- function(family, ...) { + df <- data.frame(y = c(1, 2)) + dummy_mdl <- make_stancode(y ~ 1, data = df, family = class(family)[1]) + + # get the lowered family name + family_name <- tolower(class(family)[1]) + + # Extract the Stan parameterisation from the dummy model code + lpdf_pattern <- paste0( + "target \\+= ", family_name, "_(lpdf|lpmf)\\(Y \\| ([^)]+)\\)" # nolint + ) + lpdf_match <- regexpr(lpdf_pattern, dummy_mdl) + reparam <- if (lpdf_match > 0) { + matches <- unlist(regmatches(dummy_mdl, lpdf_match)) + mu_matches <- matches[grepl("mu", matches, fixed = TRUE)] + if (length(mu_matches) > 1) { + cli_abort("Multiple Stan parameterisations found with 'mu' parameter.") + } else if (length(mu_matches) == 0) { + cli_abort("No Stan parameterisation found with 'mu' parameter.") + } + match_str <- mu_matches[1] + param <- sub( + paste0( + "target \\+= ", family_name, "_(lpdf|lpmf)\\(Y \\| " # nolint + ), "", + match_str + ) + param <- sub(")", "", param, fixed = TRUE) + family$param <- param + } else { + cli_abort( + "Unable to extract Stan parameterisation for {family_name}." + ) + } return(family) } diff --git a/R/gen.R b/R/gen.R index ebef86b3b..df7040d1b 100644 --- a/R/gen.R +++ b/R/gen.R @@ -1,5 +1,75 @@ +#' Create a function to calculate the marginalised log likelihood for double +#' censored and truncated delay distributions +#' +#' This function creates a log likelihood function that calculates the marginal +#' likelihood for a single observation by integrating over the latent primary +#' and secondary event windows. The integration is performed numerically using +#' [primarycensored::dpcens()] which handles the double censoring and truncation +#' of the delay distribution. +#' +#' The marginal likelihood accounts for uncertainty in both the primary and +#' secondary event windows by integrating over their possible values, weighted +#' by their respective uniform distributions. +#' +#' @seealso [brms::log_lik()] for details on the brms log likelihood interface. +#' +#' @inheritParams epidist_family +#' +#' @return A function that calculates the marginal log likelihood for a single +#' observation. The prep object must have the following variables: +#' * `vreal1`: relative observation time +#' * `vreal2`: primary event window +#' * `vreal3`: secondary event window +#' +#' @family gen +#' @autoglobal +#' @importFrom purrr map_dbl +epidist_gen_log_lik <- function(family) { + # Get internal brms log_lik function + log_lik_brms <- .get_brms_fn("log_lik", family) + + .log_lik <- function(i, prep) { + y <- prep$data$Y[i] + relative_obs_time <- prep$data$vreal1[i] + pwindow <- prep$data$vreal2[i] + swindow <- prep$data$vreal3[i] + + # make the prep object censored + # -1 here is equivalent to right censored in brms + prep$data$cens <- -1 + + # Calculate density for each draw using primarycensored::dpcens() + lpdf <- purrr::map_dbl(seq_len(prep$ndraws), function(draw) { + # Define pdist function that filters to current draw + pdist_draw <- function(q, i, prep, ...) { + purrr::map_dbl(q, function(x) { + prep$data$Y <- rep(x, length(prep$data$Y)) + ll <- exp(log_lik_brms(i, prep)[draw]) + return(ll) + }) + } + + primarycensored::dpcens( + x = y, + pdist = pdist_draw, + i = i, + prep = prep, + pwindow = pwindow, + swindow = swindow, + D = relative_obs_time, + dprimary = stats::dunif, + log = TRUE + ) + }) + lpdf <- brms:::log_lik_weight(lpdf, i = i, prep = prep) # nolint + return(lpdf) + } + + return(.log_lik) +} + #' Create a function to draw from the posterior predictive distribution for a -#' latent model +#' double censored and truncated delay distribution #' #' This function creates a function that draws from the posterior predictive #' distribution for a latent model using [primarycensored::rpcens()] to handle @@ -49,7 +119,7 @@ epidist_gen_posterior_predict <- function(family) { } #' Create a function to draw from the expected value of the posterior predictive -#' distribution for a latent model +#' distribution for a model #' #' This function creates a function that calculates the expected value of the #' posterior predictive distribution for a latent model. The returned function diff --git a/R/globals.R b/R/globals.R index bfa8d919d..6c525d44b 100644 --- a/R/globals.R +++ b/R/globals.R @@ -4,6 +4,7 @@ utils::globalVariables(c( "samples", # "woverlap", # "rlnorm", # + "fix", # <.replace_prior> "prior_new", # <.replace_prior> "source_new", # <.replace_prior> NULL diff --git a/R/latent_model.R b/R/latent_model.R index 96215b971..3d5c72dbb 100644 --- a/R/latent_model.R +++ b/R/latent_model.R @@ -57,6 +57,7 @@ is_epidist_latent_model <- function(data) { #' @method assert_epidist epidist_latent_model #' @family latent_model +#' @importFrom checkmate assert_names assert_numeric assert_integerish #' @export assert_epidist.epidist_latent_model <- function(data, ...) { col_names <- c( @@ -69,6 +70,7 @@ assert_epidist.epidist_latent_model <- function(data, ...) { assert_numeric(data$woverlap, lower = 0) assert_numeric(data$swindow, lower = 0) assert_numeric(data$delay, lower = 0) + assert_integerish(data$.row_id, lower = 1) } #' Create the model-specific component of an `epidist` custom family @@ -88,9 +90,12 @@ epidist_family_model.epidist_latent_model <- function( lb = c(NA, as.numeric(lapply(family$other_bounds, "[[", "lb"))), ub = c(NA, as.numeric(lapply(family$other_bounds, "[[", "ub"))), type = family$type, - vars = c("pwindow", "swindow", "vreal1"), + vars = c( + "vreal1", "vreal2", "vreal3", "pwindow_raw", "swindow_raw", + "woverlap", "wN" + ), loop = FALSE, - log_lik = epidist_gen_log_lik_latent(family), + log_lik = epidist_gen_log_lik(family), posterior_predict = epidist_gen_posterior_predict(family), posterior_epred = epidist_gen_posterior_epred(family) ) @@ -98,67 +103,6 @@ epidist_family_model.epidist_latent_model <- function( return(custom_family) } -#' Create a function to calculate the pointwise log likelihood of the latent -#' model -#' -#' This function creates a log likelihood function that accounts for the latent -#' variables in the model, including primary and secondary event windows and -#' their overlap. The returned function calculates the log likelihood for a -#' single observation by augmenting the data with the latent variables and -#' using the underlying brms log likelihood function. -#' -#' @seealso [brms::log_lik()] for details on the brms log likelihood interface. -#' -#' @inheritParams epidist_family -#' -#' @return A function that calculates the log likelihood for a single -#' observation. The prep object must have the following variables: -#' * `vreal1`: relative observation time -#' * `vreal2`: primary event window -#' * `vreal3`: secondary event window -#' -#' @family latent_model -#' @autoglobal -epidist_gen_log_lik_latent <- function(family) { - # Get internal brms log_lik function - log_lik_brms <- .get_brms_fn("log_lik", family) - - .log_lik <- function(i, prep) { - y <- prep$data$Y[i] - relative_obs_time <- prep$data$vreal1[i] - pwindow <- prep$data$vreal2[i] - swindow <- prep$data$vreal3[i] - - # Generates values of the swindow_raw and pwindow_raw, but really these - # should be extracted from prep or the fitted raws somehow. See: - # https://github.com/epinowcast/epidist/issues/267 - swindow_raw <- stats::runif(prep$ndraws) - pwindow_raw <- stats::runif(prep$ndraws) - - swindow <- swindow_raw * swindow - - # For no overlap calculate as usual, for overlap ensure pwindow < swindow - if (i %in% prep$data$noverlap) { - pwindow <- pwindow_raw * pwindow - } else { - pwindow <- pwindow_raw * swindow - } - - d <- y - pwindow + swindow - obs_time <- relative_obs_time - pwindow - # Create brms truncation upper bound - prep$data$ub <- rep(obs_time, length(prep$data$Y)) - # Update augmented data - prep$data$Y <- rep(d, length(prep$data$Y)) - - # Call internal brms log_lik function with augmented data - lpdf <- log_lik_brms(i, prep) - return(lpdf) - } - - return(.log_lik) -} - #' Define the model-specific component of an `epidist` custom formula #' #' @inheritParams epidist_formula_model @@ -168,13 +112,35 @@ epidist_gen_log_lik_latent <- function(family) { #' @export epidist_formula_model.epidist_latent_model <- function( data, formula, ...) { - # data is only used to dispatch on formula <- stats::update( formula, delay | vreal(relative_obs_time, pwindow, swindow) ~ . ) return(formula) } +#' Model specific prior distributions for latent models +#' +#' Defines prior distributions for the latent model parameters `pwindow_raw` and +#' `swindow_raw` which control the width of the observation windows. +#' +#' @inheritParams epidist +#' @importFrom brms set_prior +#' @family latent_model +#' @export +epidist_model_prior.epidist_latent_model <- function(data, formula, ...) { + priors <- prior( + "pwindow_raw ~ uniform(0, 1);", + dpar = "pwindow_raw", + check = FALSE + ) + + prior( + "swindow_raw ~ uniform(0, 1);", + dpar = "swindow_raw", + check = FALSE + ) + return(priors) +} + #' @method epidist_stancode epidist_latent_model #' @importFrom brms stanvar #' @family latent_model @@ -212,12 +178,8 @@ epidist_stancode.epidist_latent_model <- function( fixed = TRUE ) - # dpars_B refers to the insertion into the lpdf call - # For some families, we tranform brms dpars to match Stan parameterisation stanvars_functions[[1]]$scode <- gsub( - "dpars_B", - toString(family$reparam), - stanvars_functions[[1]]$scode, + "dpars_B", family$param, stanvars_functions[[1]]$scode, fixed = TRUE ) @@ -227,12 +189,6 @@ epidist_stancode.epidist_latent_model <- function( x = nrow(filter(data, woverlap > 0)), name = "wN" ) + - stanvar( - block = "data", - scode = "array[N - wN] int noverlap;", - x = filter(data, woverlap == 0)$.row_id, - name = "noverlap" - ) + stanvar( block = "data", scode = "array[wN] int woverlap;", @@ -242,21 +198,15 @@ epidist_stancode.epidist_latent_model <- function( stanvars_parameters <- stanvar( block = "parameters", - scode = .stan_chunk(file.path("latent_model", "parameters.stan")) - ) - - stanvars_tparameters <- stanvar( - block = "tparameters", - scode = .stan_chunk(file.path("latent_model", "tparameters.stan")) - ) - - stanvars_priors <- stanvar( - block = "model", - scode = .stan_chunk(file.path("latent_model", "priors.stan")) - ) + scode = "vector[N] pwindow_raw;" + ) + + stanvar( + block = "parameters", + scode = "vector[N] swindow_raw;" + ) stanvars_all <- stanvars_version + stanvars_functions + stanvars_data + - stanvars_parameters + stanvars_tparameters + stanvars_priors + stanvars_parameters return(stanvars_all) } diff --git a/R/prior.R b/R/prior.R index 015b16903..1a104b410 100644 --- a/R/prior.R +++ b/R/prior.R @@ -1,37 +1,58 @@ -#' Define prior distributions using `brms` defaults, model specific priors, -#' family specific priors, and user provided priors +#' Define custom prior distributions for epidist models #' -#' This function obtains the `brms` default prior distributions for a particular -#' model, then replaces these prior distributions using: -#' 1. Model specific prior distributions from [epidist_model_prior()] -#' 2. Family specific prior distributions from [epidist_family_prior()] -#' 3. User provided prior distributions -#' Each element of this list overwrites previous elements, such that user -#' provided prior distribution have the highest priority. At the third stage, -#' if a prior distribution is provided which is not included in the model, then -#' a warning will be shown. To prevent this warning, do not pass prior -#' distributions for parameters which are not in the model. +#' This function combines model specific prior distributions from +#' [epidist_model_prior()], family specific prior distributions from +#' [epidist_family_prior()], and user provided prior distributions into a single +#' set of custom priors. Each element overwrites previous elements, such that +#' user provided prior distributions have the highest priority. If a user prior +#' distribution is provided which is not included in the model, a warning will +#' be shown. +#' +#' Note that the matching of priors is imperfect as it does not use brms' +#' internal prior matching functionality. For example, it cannot distinguish +#' between a prior for all coefficients (class = "b") and a prior for a +#' specific coefficient (class = "b" and coef specified). #' #' @inheritParams epidist #' @param family A description of the response distribution and link function to #' be used in the model created using [epidist_family()]. #' @param formula A symbolic description of the model to be fitted created using #' [epidist_formula()]. +#' @param merge If `TRUE` then merge new priors with existing ones, if `FALSE` +#' only use new priors. Defaults to `TRUE`. This may be useful if the built in +#' approaches for merging priors are not flexible enough for a particular use +#' case. +#' @param enforce_presence If `TRUE` then only allow user priors that match +#' existing default priors. If `FALSE` then allow user priors that are not +#' present in the default set. Defaults to `FALSE`. +#' @return A `brmsprior` object containing the combined custom prior +#' distributions. #' @rdname epidist_prior #' @family prior #' @export -epidist_prior <- function(data, family, formula, prior) { +epidist_prior <- function(data, family, formula, prior, merge = TRUE, + enforce_presence = FALSE) { assert_epidist(data) default <- brms::default_prior(formula, data = data) model <- epidist_model_prior(data, formula) + if (!is.null(model)) { + model$source <- "model" + } family <- epidist_family_prior(family, formula) if (!is.null(family)) { family$source <- "family" - family[is.na(family)] <- "" # brms likes empty over NA - family[family == "NA"] <- NA # To keep particular NA } - internal_prior <- Reduce(.replace_prior, list(default, model, family)) - prior <- .replace_prior(internal_prior, prior, warn = TRUE) + custom <- .replace_prior( + family, model, + merge = TRUE, enforce_presence = FALSE + ) + internal <- .replace_prior(default, custom, merge = TRUE) + prior <- .replace_prior( + internal, prior, + warn = TRUE, merge = merge, + enforce_presence = enforce_presence + ) + return(prior) } diff --git a/R/utils.R b/R/utils.R index 9fd6a419e..826beb3e9 100644 --- a/R/utils.R +++ b/R/utils.R @@ -44,49 +44,125 @@ #' Replace `brms` prior distributions #' -#' This function takes `old_prior` and replaces any prior distributions -#' contained in it by the corresponding prior distribution in `prior`. If there -#' is a prior distribution in `prior` with no match in `old_prior` then this -#' function can optionally give a warning. -#' -#' @param old_prior One or more prior distributions in the class `brmsprior` +#' This function takes an existing set of prior distributions and updates them +#' with new prior specifications. It matches priors based on their parameter +#' class, coefficient, group, response, distributional parameter, and non-linear +#' parameter. Any new priors that don't match existing ones can optionally +#' trigger a warning. +#' +#' Prior distributions can be specified in two ways: +#' 1. Using the standard `brms` prior specification format. These priors are +#' replaced based on matching parameter metadata (class, coefficient, group, +#' etc.). +#' 2. Using custom set priors with the syntax `parameter ~ distribution`. These +#' will only remove existing custom priors for the same parameter name but +#' will not affect priors set via the standard `brms` specification format. +#' Custom priors are excluded from the metadata-based joining process. +#' +#' @param old_prior One or more prior distributions in the class `brmsprior` to +#' be updated #' @param prior One or more prior distributions in the class `brmsprior` -#' @param warn If `TRUE` then a warning will be displayed if a `new_prior` is -#' provided for which there is no matching `old_prior`. Defaults to `FALSE` +#' containing the new specifications. Can include custom set priors using the +#' syntax `parameter ~ distribution` +#' @param warn If `TRUE` then a warning will be displayed if a prior in `prior` +#' has no match in `old_prior`. Defaults to `FALSE` +#' @param merge If `TRUE` then merge new priors with existing ones, if `FALSE` +#' only use new priors. Defaults to `TRUE` +#' @param enforce_presence If `TRUE` then only keep rows that have both old and +#' new priors. If `FALSE` then keep all rows but use new priors where +#' available, otherwise keep old priors. Defaults to `TRUE` #' @autoglobal -#' @importFrom dplyr full_join filter select +#' @importFrom dplyr full_join filter select mutate bind_rows +#' @importFrom brms as.brmsprior #' @keywords internal -.replace_prior <- function(old_prior, prior, warn = FALSE) { +.replace_prior <- function(old_prior, prior, warn = FALSE, merge = TRUE, + enforce_presence = TRUE) { + if (!isTRUE(merge)) { + return(prior) + } + if (is.null(prior)) { return(old_prior) } - cols <- c("class", "coef", "group", "resp", "dpar", "nlpar", "lb", "ub") - prior <- dplyr::full_join( - old_prior, prior, - by = cols, suffix = c("_old", "_new") - ) - if (anyNA(prior$prior_old)) { - missing_prior <- utils::capture.output(print( - prior |> - filter(is.na(.data$prior_old)) |> - select(prior = prior_new, dplyr::all_of(cols), source = source_new) - )) - if (warn) { - msg <- c( - "!" = "One or more priors have no match in existing parameters:", - missing_prior, - "i" = "To remove this warning consider changing prior specification." # nolint - ) - cli_warn(message = msg) - } + if (is.null(old_prior)) { + return(prior) + } + + # Find priors defined with ~ in prior column + tilde_priors <- prior[grepl("~", prior$prior, fix, fixed = TRUE), ] + if (nrow(tilde_priors) > 0) { + # Extract parameter names from left side of ~ + param_names <- gsub("\\s*~.*$", "", tilde_priors$prior) + + # Remove matching parameter priors from old_prior + old_prior <- old_prior[ + !grepl(paste(param_names, collapse = "|"), old_prior$prior), + ] } - prior <- prior |> - filter(!is.na(.data$prior_old), !is.na(.data$prior_new)) |> - select(prior = prior_new, dplyr::all_of(cols), source = source_new) + # Hold out manual priors + hold_prior <- prior[grepl("~", prior$prior, fixed = TRUE), ] + hold_prior_old <- old_prior[grepl("~", old_prior$prior, fixed = TRUE), ] - return(prior) + prior <- prior[!grepl("~", prior$prior, fixed = TRUE), ] + old_prior <- old_prior[!grepl("~", old_prior$prior, fixed = TRUE), ] + + join_cols <- c("class", "coef", "group", "resp", "dpar", "nlpar") + cols <- c(join_cols, "lb", "ub") + + if (nrow(prior) == 0 && nrow(old_prior) == 0) { + # If no non-manual priors found, just combine manual priors + cli::cli_inform("No non-manual priors found, only combining manual priors") + prior <- bind_rows(hold_prior, hold_prior_old) + } else { + prior <- dplyr::full_join( + old_prior, prior, + by = join_cols, suffix = c("_old", "_new") + ) + + if (anyNA(prior$prior_old)) { + missing_prior <- utils::capture.output(print( + prior |> + filter(is.na(.data$prior_old)) |> + as.data.frame() + )) + if (warn) { + msg <- c( + "!" = "One or more priors have no match in existing parameters:", + missing_prior, + "i" = "To remove this warning consider changing prior specification." # nolint + ) + cli_warn(message = msg) + } + } + + # use new lb and ub if prior_new is present + prior <- prior |> + mutate( + lb = ifelse(!is.na(.data$prior_new), .data$lb_new, .data$lb_old), + ub = ifelse(!is.na(.data$prior_new), .data$ub_new, .data$ub_old) + ) + + # only keep rows that have both old and new priors + if (isTRUE(enforce_presence)) { + prior <- prior |> + filter(!is.na(.data$prior_old), !is.na(.data$prior_new)) |> + select(prior = prior_new, dplyr::all_of(cols), source = source_new) + } else { + prior <- prior |> + mutate(prior = ifelse( + !is.na(.data$prior_new), .data$prior_new, .data$prior_old + )) |> + mutate(source = ifelse( + !is.na(.data$prior_new), .data$source_new, .data$source_old + )) |> + select(prior, dplyr::all_of(cols), source) + } + + prior <- bind_rows(prior, hold_prior, hold_prior_old) + } + return(as.brmsprior(prior)) } #' Additional distributional parameter information for `brms` families diff --git a/inst/stan/latent_model/functions.stan b/inst/stan/latent_model/functions.stan index 6aaa1dce6..3134f7f95 100644 --- a/inst/stan/latent_model/functions.stan +++ b/inst/stan/latent_model/functions.stan @@ -1,12 +1,39 @@ -// Here the strings -// * family -// * dpars_A -// * dpars_B -// are/have been replaced using regex - -real latent_family_lpdf(vector y, dpars_A, vector pwindow, - vector swindow, array[] real relative_obs_t) { +/** + * Compute the log probability density function for a latent model with censoring + * + * This function is designed to be read into R where: + * - 'family' is replaced with the target distribution (e.g., 'lognormal') + * - 'dpars_A' is replaced with multiple parameters in the format + * "vector|real paramname1, vector|real paramname2, ..." depending on whether + * each parameter has a model. This includes distribution parameters. + * - 'dpars_B' is replaced with the same parameters as dpars_A but with window + * indices removed. + * + * @param y Vector of observed values (delays) + * @param dpars_A Distribution parameters (replaced via regex) + * @param relative_obs_t Array of observation times relative to primary window + * start + * @param pwindow_width Array of primary window widths (actual time scale) + * @param swindow_width Array of secondary window widths (actual time scale) + * @param pwindow_raw Vector of primary window positions (0-1 scale) + * @param swindow_raw Vector of secondary window positions (0-1 scale) + * @param woverlap Array of indices for overlapping windows + * @param wN Number of overlapping windows + * + * @return Log probability density with censoring adjustment for latent model + */ +real latent_family_lpdf(vector y, dpars_A, + array[] real relative_obs_t, + array[] real pwindow_width, array[] real swindow_width, + vector pwindow_raw, vector swindow_raw, + array[] int woverlap, int wN) { int n = num_elements(y); + vector[n] pwindow = to_vector(pwindow_width) .* pwindow_raw; + vector[n] swindow = to_vector(swindow_width) .* swindow_raw; + + if (wN) { + pwindow[woverlap] = swindow[woverlap] .* pwindow_raw[woverlap]; + } vector[n] d = y - pwindow + swindow; vector[n] obs_time = to_vector(relative_obs_t) - pwindow; return family_lpdf(d | dpars_B) - family_lcdf(obs_time | dpars_B); diff --git a/inst/stan/latent_model/parameters.stan b/inst/stan/latent_model/parameters.stan deleted file mode 100644 index 066703030..000000000 --- a/inst/stan/latent_model/parameters.stan +++ /dev/null @@ -1,2 +0,0 @@ -vector[N] swindow_raw; -vector[N] pwindow_raw; diff --git a/inst/stan/latent_model/priors.stan b/inst/stan/latent_model/priors.stan deleted file mode 100644 index 6b141121c..000000000 --- a/inst/stan/latent_model/priors.stan +++ /dev/null @@ -1,2 +0,0 @@ -swindow_raw ~ uniform(0, 1); -pwindow_raw ~ uniform(0, 1); diff --git a/inst/stan/latent_model/tparameters.stan b/inst/stan/latent_model/tparameters.stan deleted file mode 100644 index d7acadd8f..000000000 --- a/inst/stan/latent_model/tparameters.stan +++ /dev/null @@ -1,7 +0,0 @@ -vector[N] pwindow; -vector[N] swindow; -swindow = to_vector(vreal3) .* swindow_raw; -pwindow[noverlap] = to_vector(vreal2[noverlap]) .* pwindow_raw[noverlap]; -if (wN) { - pwindow[woverlap] = swindow[woverlap] .* pwindow_raw[woverlap]; -} diff --git a/man/as_epidist_latent_model.Rd b/man/as_epidist_latent_model.Rd index 2aa6e600b..9f9a419aa 100644 --- a/man/as_epidist_latent_model.Rd +++ b/man/as_epidist_latent_model.Rd @@ -17,7 +17,7 @@ Other latent_model: \code{\link{as_epidist_latent_model.epidist_linelist_data}()}, \code{\link{epidist_family_model.epidist_latent_model}()}, \code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{is_epidist_latent_model}()}, \code{\link{new_epidist_latent_model}()} } diff --git a/man/as_epidist_latent_model.epidist_linelist_data.Rd b/man/as_epidist_latent_model.epidist_linelist_data.Rd index 03e839d66..3b91956df 100644 --- a/man/as_epidist_latent_model.epidist_linelist_data.Rd +++ b/man/as_epidist_latent_model.epidist_linelist_data.Rd @@ -17,7 +17,7 @@ Other latent_model: \code{\link{as_epidist_latent_model}()}, \code{\link{epidist_family_model.epidist_latent_model}()}, \code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{is_epidist_latent_model}()}, \code{\link{new_epidist_latent_model}()} } diff --git a/man/dot-replace_prior.Rd b/man/dot-replace_prior.Rd index 3c70a95de..9b9d576ca 100644 --- a/man/dot-replace_prior.Rd +++ b/man/dot-replace_prior.Rd @@ -4,20 +4,49 @@ \alias{.replace_prior} \title{Replace \code{brms} prior distributions} \usage{ -.replace_prior(old_prior, prior, warn = FALSE) +.replace_prior( + old_prior, + prior, + warn = FALSE, + merge = TRUE, + enforce_presence = TRUE +) } \arguments{ -\item{old_prior}{One or more prior distributions in the class \code{brmsprior}} +\item{old_prior}{One or more prior distributions in the class \code{brmsprior} to +be updated} -\item{prior}{One or more prior distributions in the class \code{brmsprior}} +\item{prior}{One or more prior distributions in the class \code{brmsprior} +containing the new specifications. Can include custom set priors using the +syntax \code{parameter ~ distribution}} -\item{warn}{If \code{TRUE} then a warning will be displayed if a \code{new_prior} is -provided for which there is no matching \code{old_prior}. Defaults to \code{FALSE}} +\item{warn}{If \code{TRUE} then a warning will be displayed if a prior in \code{prior} +has no match in \code{old_prior}. Defaults to \code{FALSE}} + +\item{merge}{If \code{TRUE} then merge new priors with existing ones, if \code{FALSE} +only use new priors. Defaults to \code{TRUE}} + +\item{enforce_presence}{If \code{TRUE} then only keep rows that have both old and +new priors. If \code{FALSE} then keep all rows but use new priors where +available, otherwise keep old priors. Defaults to \code{TRUE}} } \description{ -This function takes \code{old_prior} and replaces any prior distributions -contained in it by the corresponding prior distribution in \code{prior}. If there -is a prior distribution in \code{prior} with no match in \code{old_prior} then this -function can optionally give a warning. +This function takes an existing set of prior distributions and updates them +with new prior specifications. It matches priors based on their parameter +class, coefficient, group, response, distributional parameter, and non-linear +parameter. Any new priors that don't match existing ones can optionally +trigger a warning. +} +\details{ +Prior distributions can be specified in two ways: +\enumerate{ +\item Using the standard \code{brms} prior specification format. These priors are +replaced based on matching parameter metadata (class, coefficient, group, +etc.). +\item Using custom set priors with the syntax \code{parameter ~ distribution}. These +will only remove existing custom priors for the same parameter name but +will not affect priors set via the standard \code{brms} specification format. +Custom priors are excluded from the metadata-based joining process. +} } \keyword{internal} diff --git a/man/epidist.Rd b/man/epidist.Rd index 27e754a90..9a7d6b374 100644 --- a/man/epidist.Rd +++ b/man/epidist.Rd @@ -9,6 +9,7 @@ epidist( formula = mu ~ 1, family = lognormal(), prior = NULL, + merge_priors = TRUE, fn = brms::brm, ... ) @@ -31,7 +32,14 @@ reexported as part of \code{epidist}.} \item{prior}{One or more \code{brmsprior} objects created by \code{\link[brms:set_prior]{brms::set_prior()}} or related functions. These priors are passed to \code{\link[=epidist_prior]{epidist_prior()}} in the -\code{prior} argument.} +\code{prior} argument. Some models have default priors that are automatically +added (see \code{\link[=epidist_model_prior]{epidist_model_prior()}}). These can be merged with user-provided +priors using the \code{merge_priors} argument.} + +\item{merge_priors}{If \code{TRUE} then merge user priors with default priors, if +\code{FALSE} only use user priors. Defaults to \code{TRUE}. This may be useful if +the built in approaches for merging priors are not flexible enough for a +particular use case.} \item{fn}{The internal function to be called. By default this is \code{\link[brms:brm]{brms::brm()}} which performs inference for the specified model. Other options diff --git a/man/epidist_family.Rd b/man/epidist_family.Rd index e7272bec3..1f0ebec9b 100644 --- a/man/epidist_family.Rd +++ b/man/epidist_family.Rd @@ -28,8 +28,7 @@ be transparent about what happens inside of a call to \code{\link[=epidist]{epid Other family: \code{\link{epidist_family_model}()}, \code{\link{epidist_family_model.default}()}, -\code{\link{epidist_family_reparam}()}, -\code{\link{epidist_family_reparam.default}()}, -\code{\link{epidist_family_reparam.gamma}()} +\code{\link{epidist_family_param}()}, +\code{\link{epidist_family_param.default}()} } \concept{family} diff --git a/man/epidist_family_model.Rd b/man/epidist_family_model.Rd index 708d1e4db..d616f9364 100644 --- a/man/epidist_family_model.Rd +++ b/man/epidist_family_model.Rd @@ -32,9 +32,8 @@ The model-specific parts of an \code{epidist_formula()} call Other family: \code{\link{epidist_family}()}, \code{\link{epidist_family_model.default}()}, -\code{\link{epidist_family_reparam}()}, -\code{\link{epidist_family_reparam.default}()}, -\code{\link{epidist_family_reparam.gamma}()} +\code{\link{epidist_family_param}()}, +\code{\link{epidist_family_param.default}()} Other formula: \code{\link{epidist_formula}()}, diff --git a/man/epidist_family_model.default.Rd b/man/epidist_family_model.default.Rd index 13a036730..4a8b37924 100644 --- a/man/epidist_family_model.default.Rd +++ b/man/epidist_family_model.default.Rd @@ -21,8 +21,7 @@ Default method for defining a model specific family Other family: \code{\link{epidist_family}()}, \code{\link{epidist_family_model}()}, -\code{\link{epidist_family_reparam}()}, -\code{\link{epidist_family_reparam.default}()}, -\code{\link{epidist_family_reparam.gamma}()} +\code{\link{epidist_family_param}()}, +\code{\link{epidist_family_param.default}()} } \concept{family} diff --git a/man/epidist_family_model.epidist_latent_model.Rd b/man/epidist_family_model.epidist_latent_model.Rd index 6794976a2..290289960 100644 --- a/man/epidist_family_model.epidist_latent_model.Rd +++ b/man/epidist_family_model.epidist_latent_model.Rd @@ -22,7 +22,7 @@ Other latent_model: \code{\link{as_epidist_latent_model}()}, \code{\link{as_epidist_latent_model.epidist_linelist_data}()}, \code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{is_epidist_latent_model}()}, \code{\link{new_epidist_latent_model}()} } diff --git a/man/epidist_family_reparam.Rd b/man/epidist_family_param.Rd similarity index 83% rename from man/epidist_family_reparam.Rd rename to man/epidist_family_param.Rd index 13a39547b..0e079a555 100644 --- a/man/epidist_family_reparam.Rd +++ b/man/epidist_family_param.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/family.R -\name{epidist_family_reparam} -\alias{epidist_family_reparam} +\name{epidist_family_param} +\alias{epidist_family_param} \title{Reparameterise an \code{epidist} family to align \code{brms} and Stan} \usage{ -epidist_family_reparam(family, ...) +epidist_family_param(family, ...) } \arguments{ \item{family}{A description of the response distribution and link function to @@ -24,7 +24,6 @@ Other family: \code{\link{epidist_family}()}, \code{\link{epidist_family_model}()}, \code{\link{epidist_family_model.default}()}, -\code{\link{epidist_family_reparam.default}()}, -\code{\link{epidist_family_reparam.gamma}()} +\code{\link{epidist_family_param.default}()} } \concept{family} diff --git a/man/epidist_family_param.default.Rd b/man/epidist_family_param.default.Rd new file mode 100644 index 000000000..6959ae776 --- /dev/null +++ b/man/epidist_family_param.default.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/family.R +\name{epidist_family_param.default} +\alias{epidist_family_param.default} +\title{Default method for families which do not require a reparameterisation} +\usage{ +\method{epidist_family_param}{default}(family, ...) +} +\arguments{ +\item{family}{A brms family object containing at minimum a \code{family} element +specifying the distribution family name} + +\item{...}{Additional arguments passed to methods (not used)} +} +\value{ +The input family object with an additional \code{param} element containing +the Stan parameter ordering as a string +} +\description{ +This function extracts the Stan parameterisation for a given brms family by +creating a dummy model and parsing its Stan code. It looks for the log +probability density function (lpdf) call in the Stan code and extracts the +parameter order used by Stan. This is needed because brms and Stan may use +different parameter orderings for the same distribution. +} +\details{ +The function works by: +\enumerate{ +\item Creating a minimal dummy model using the specified family +\item Extracting the Stan code for this model +\item Finding the lpdf function call for the family +\item Parsing out the parameter ordering used in Stan +\item Adding this as the \code{param} element to the family object +} +} +\seealso{ +Other family: +\code{\link{epidist_family}()}, +\code{\link{epidist_family_model}()}, +\code{\link{epidist_family_model.default}()}, +\code{\link{epidist_family_param}()} +} +\concept{family} diff --git a/man/epidist_family_reparam.default.Rd b/man/epidist_family_reparam.default.Rd deleted file mode 100644 index 8ab4d39ce..000000000 --- a/man/epidist_family_reparam.default.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/family.R -\name{epidist_family_reparam.default} -\alias{epidist_family_reparam.default} -\title{Default method for families which do not require a reparameterisation} -\usage{ -\method{epidist_family_reparam}{default}(family, ...) -} -\arguments{ -\item{family}{A description of the response distribution and link function to -be used in the model. Every family function has a link argument allowing -users to specify the link function to be applied on the response variable. -If not specified, default links are used. For details of all supported -families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also -reexported as part of \code{epidist}.} - -\item{...}{Additional arguments passed to \code{fn} method.} -} -\description{ -Default method for families which do not require a reparameterisation -} -\seealso{ -Other family: -\code{\link{epidist_family}()}, -\code{\link{epidist_family_model}()}, -\code{\link{epidist_family_model.default}()}, -\code{\link{epidist_family_reparam}()}, -\code{\link{epidist_family_reparam.gamma}()} -} -\concept{family} diff --git a/man/epidist_family_reparam.gamma.Rd b/man/epidist_family_reparam.gamma.Rd deleted file mode 100644 index 090cae603..000000000 --- a/man/epidist_family_reparam.gamma.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/family.R -\name{epidist_family_reparam.gamma} -\alias{epidist_family_reparam.gamma} -\title{Reparameterisation for the gamma family} -\usage{ -\method{epidist_family_reparam}{gamma}(family, ...) -} -\arguments{ -\item{family}{A description of the response distribution and link function to -be used in the model. Every family function has a link argument allowing -users to specify the link function to be applied on the response variable. -If not specified, default links are used. For details of all supported -families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also -reexported as part of \code{epidist}.} - -\item{...}{Additional arguments passed to \code{fn} method.} -} -\description{ -Reparameterisation for the gamma family -} -\seealso{ -Other family: -\code{\link{epidist_family}()}, -\code{\link{epidist_family_model}()}, -\code{\link{epidist_family_model.default}()}, -\code{\link{epidist_family_reparam}()}, -\code{\link{epidist_family_reparam.default}()} -} -\concept{family} diff --git a/man/epidist_formula_model.epidist_latent_model.Rd b/man/epidist_formula_model.epidist_latent_model.Rd index 6a5cdfc46..40314245e 100644 --- a/man/epidist_formula_model.epidist_latent_model.Rd +++ b/man/epidist_formula_model.epidist_latent_model.Rd @@ -25,7 +25,7 @@ Other latent_model: \code{\link{as_epidist_latent_model}()}, \code{\link{as_epidist_latent_model.epidist_linelist_data}()}, \code{\link{epidist_family_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{is_epidist_latent_model}()}, \code{\link{new_epidist_latent_model}()} } diff --git a/man/epidist_gen_log_lik.Rd b/man/epidist_gen_log_lik.Rd new file mode 100644 index 000000000..cc989ff12 --- /dev/null +++ b/man/epidist_gen_log_lik.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gen.R +\name{epidist_gen_log_lik} +\alias{epidist_gen_log_lik} +\title{Create a function to calculate the marginalised log likelihood for double +censored and truncated delay distributions} +\usage{ +epidist_gen_log_lik(family) +} +\arguments{ +\item{family}{A description of the response distribution and link function to +be used in the model. Every family function has a link argument allowing +users to specify the link function to be applied on the response variable. +If not specified, default links are used. For details of all supported +families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also +reexported as part of \code{epidist}.} +} +\value{ +A function that calculates the marginal log likelihood for a single +observation. The prep object must have the following variables: +\itemize{ +\item \code{vreal1}: relative observation time +\item \code{vreal2}: primary event window +\item \code{vreal3}: secondary event window +} +} +\description{ +This function creates a log likelihood function that calculates the marginal +likelihood for a single observation by integrating over the latent primary +and secondary event windows. The integration is performed numerically using +\code{\link[primarycensored:dprimarycensored]{primarycensored::dpcens()}} which handles the double censoring and truncation +of the delay distribution. +} +\details{ +The marginal likelihood accounts for uncertainty in both the primary and +secondary event windows by integrating over their possible values, weighted +by their respective uniform distributions. +} +\seealso{ +\code{\link[brms:log_lik.brmsfit]{brms::log_lik()}} for details on the brms log likelihood interface. + +Other gen: +\code{\link{epidist_gen_posterior_epred}()}, +\code{\link{epidist_gen_posterior_predict}()} +} +\concept{gen} diff --git a/man/epidist_gen_log_lik_latent.Rd b/man/epidist_gen_log_lik_latent.Rd deleted file mode 100644 index 5fc1e1dee..000000000 --- a/man/epidist_gen_log_lik_latent.Rd +++ /dev/null @@ -1,45 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/latent_model.R -\name{epidist_gen_log_lik_latent} -\alias{epidist_gen_log_lik_latent} -\title{Create a function to calculate the pointwise log likelihood of the latent -model} -\usage{ -epidist_gen_log_lik_latent(family) -} -\arguments{ -\item{family}{A description of the response distribution and link function to -be used in the model. Every family function has a link argument allowing -users to specify the link function to be applied on the response variable. -If not specified, default links are used. For details of all supported -families see \code{\link[=brmsfamily]{brmsfamily()}}. Commonly used, such as \code{\link[=lognormal]{lognormal()}}, are also -reexported as part of \code{epidist}.} -} -\value{ -A function that calculates the log likelihood for a single -observation. The prep object must have the following variables: -\itemize{ -\item \code{vreal1}: relative observation time -\item \code{vreal2}: primary event window -\item \code{vreal3}: secondary event window -} -} -\description{ -This function creates a log likelihood function that accounts for the latent -variables in the model, including primary and secondary event windows and -their overlap. The returned function calculates the log likelihood for a -single observation by augmenting the data with the latent variables and -using the underlying brms log likelihood function. -} -\seealso{ -\code{\link[brms:log_lik.brmsfit]{brms::log_lik()}} for details on the brms log likelihood interface. - -Other latent_model: -\code{\link{as_epidist_latent_model}()}, -\code{\link{as_epidist_latent_model.epidist_linelist_data}()}, -\code{\link{epidist_family_model.epidist_latent_model}()}, -\code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{is_epidist_latent_model}()}, -\code{\link{new_epidist_latent_model}()} -} -\concept{latent_model} diff --git a/man/epidist_gen_posterior_epred.Rd b/man/epidist_gen_posterior_epred.Rd index 05d072457..87f66488e 100644 --- a/man/epidist_gen_posterior_epred.Rd +++ b/man/epidist_gen_posterior_epred.Rd @@ -3,7 +3,7 @@ \name{epidist_gen_posterior_epred} \alias{epidist_gen_posterior_epred} \title{Create a function to draw from the expected value of the posterior predictive -distribution for a latent model} +distribution for a model} \usage{ epidist_gen_posterior_epred(family) } @@ -32,6 +32,7 @@ values for latent models. \code{brms}. Other gen: +\code{\link{epidist_gen_log_lik}()}, \code{\link{epidist_gen_posterior_predict}()} } \concept{gen} diff --git a/man/epidist_gen_posterior_predict.Rd b/man/epidist_gen_posterior_predict.Rd index b572d4c6e..0f349a3a2 100644 --- a/man/epidist_gen_posterior_predict.Rd +++ b/man/epidist_gen_posterior_predict.Rd @@ -3,7 +3,7 @@ \name{epidist_gen_posterior_predict} \alias{epidist_gen_posterior_predict} \title{Create a function to draw from the posterior predictive distribution for a -latent model} +double censored and truncated delay distribution} \usage{ epidist_gen_posterior_predict(family) } @@ -37,6 +37,7 @@ censoring and truncation. The returned function takes a \code{prep} argument fro \code{brms}, \code{\link[primarycensored:rprimarycensored]{primarycensored::rpcens()}} for details on the censoring approach Other gen: +\code{\link{epidist_gen_log_lik}()}, \code{\link{epidist_gen_posterior_epred}()} } \concept{gen} diff --git a/man/epidist_model_prior.epidist_latent_model.Rd b/man/epidist_model_prior.epidist_latent_model.Rd new file mode 100644 index 000000000..dab935439 --- /dev/null +++ b/man/epidist_model_prior.epidist_latent_model.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/latent_model.R +\name{epidist_model_prior.epidist_latent_model} +\alias{epidist_model_prior.epidist_latent_model} +\title{Model specific prior distributions for latent models} +\usage{ +\method{epidist_model_prior}{epidist_latent_model}(data, formula, ...) +} +\arguments{ +\item{data}{An object with class corresponding to an implemented model.} + +\item{formula}{An object of class \link[stats:formula]{stats::formula} or \link[brms:brmsformula]{brms::brmsformula} +(or one that can be coerced to those classes). A symbolic description of the +model to be fitted. A formula must be provided for the distributional +parameter \code{mu}, and may optionally be provided for other distributional +parameters.} + +\item{...}{Additional arguments passed to \code{fn} method.} +} +\description{ +Defines prior distributions for the latent model parameters \code{pwindow_raw} and +\code{swindow_raw} which control the width of the observation windows. +} +\seealso{ +Other latent_model: +\code{\link{as_epidist_latent_model}()}, +\code{\link{as_epidist_latent_model.epidist_linelist_data}()}, +\code{\link{epidist_family_model.epidist_latent_model}()}, +\code{\link{epidist_formula_model.epidist_latent_model}()}, +\code{\link{is_epidist_latent_model}()}, +\code{\link{new_epidist_latent_model}()} +} +\concept{latent_model} diff --git a/man/epidist_prior.Rd b/man/epidist_prior.Rd index 7863ed3b7..4ac2dec5e 100644 --- a/man/epidist_prior.Rd +++ b/man/epidist_prior.Rd @@ -2,10 +2,16 @@ % Please edit documentation in R/prior.R \name{epidist_prior} \alias{epidist_prior} -\title{Define prior distributions using \code{brms} defaults, model specific priors, -family specific priors, and user provided priors} +\title{Define custom prior distributions for epidist models} \usage{ -epidist_prior(data, family, formula, prior) +epidist_prior( + data, + family, + formula, + prior, + merge = TRUE, + enforce_presence = FALSE +) } \arguments{ \item{data}{An object with class corresponding to an implemented model.} @@ -18,21 +24,37 @@ be used in the model created using \code{\link[=epidist_family]{epidist_family() \item{prior}{One or more \code{brmsprior} objects created by \code{\link[brms:set_prior]{brms::set_prior()}} or related functions. These priors are passed to \code{\link[=epidist_prior]{epidist_prior()}} in the -\code{prior} argument.} +\code{prior} argument. Some models have default priors that are automatically +added (see \code{\link[=epidist_model_prior]{epidist_model_prior()}}). These can be merged with user-provided +priors using the \code{merge_priors} argument.} + +\item{merge}{If \code{TRUE} then merge new priors with existing ones, if \code{FALSE} +only use new priors. Defaults to \code{TRUE}. This may be useful if the built in +approaches for merging priors are not flexible enough for a particular use +case.} + +\item{enforce_presence}{If \code{TRUE} then only allow user priors that match +existing default priors. If \code{FALSE} then allow user priors that are not +present in the default set. Defaults to \code{FALSE}.} +} +\value{ +A \code{brmsprior} object containing the combined custom prior +distributions. } \description{ -This function obtains the \code{brms} default prior distributions for a particular -model, then replaces these prior distributions using: -\enumerate{ -\item Model specific prior distributions from \code{\link[=epidist_model_prior]{epidist_model_prior()}} -\item Family specific prior distributions from \code{\link[=epidist_family_prior]{epidist_family_prior()}} -\item User provided prior distributions -Each element of this list overwrites previous elements, such that user -provided prior distribution have the highest priority. At the third stage, -if a prior distribution is provided which is not included in the model, then -a warning will be shown. To prevent this warning, do not pass prior -distributions for parameters which are not in the model. +This function combines model specific prior distributions from +\code{\link[=epidist_model_prior]{epidist_model_prior()}}, family specific prior distributions from +\code{\link[=epidist_family_prior]{epidist_family_prior()}}, and user provided prior distributions into a single +set of custom priors. Each element overwrites previous elements, such that +user provided prior distributions have the highest priority. If a user prior +distribution is provided which is not included in the model, a warning will +be shown. } +\details{ +Note that the matching of priors is imperfect as it does not use brms' +internal prior matching functionality. For example, it cannot distinguish +between a prior for all coefficients (class = "b") and a prior for a +specific coefficient (class = "b" and coef specified). } \seealso{ Other prior: diff --git a/man/is_epidist_latent_model.Rd b/man/is_epidist_latent_model.Rd index 647725836..c3d332be1 100644 --- a/man/is_epidist_latent_model.Rd +++ b/man/is_epidist_latent_model.Rd @@ -18,7 +18,7 @@ Other latent_model: \code{\link{as_epidist_latent_model.epidist_linelist_data}()}, \code{\link{epidist_family_model.epidist_latent_model}()}, \code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{new_epidist_latent_model}()} } \concept{latent_model} diff --git a/man/new_epidist_latent_model.Rd b/man/new_epidist_latent_model.Rd index fe5e672a7..8658161c4 100644 --- a/man/new_epidist_latent_model.Rd +++ b/man/new_epidist_latent_model.Rd @@ -21,7 +21,7 @@ Other latent_model: \code{\link{as_epidist_latent_model.epidist_linelist_data}()}, \code{\link{epidist_family_model.epidist_latent_model}()}, \code{\link{epidist_formula_model.epidist_latent_model}()}, -\code{\link{epidist_gen_log_lik_latent}()}, +\code{\link{epidist_model_prior.epidist_latent_model}()}, \code{\link{is_epidist_latent_model}()} } \concept{latent_model} diff --git a/tests/testthat/_snaps/utils.md b/tests/testthat/_snaps/utils.md deleted file mode 100644 index fee838875..000000000 --- a/tests/testthat/_snaps/utils.md +++ /dev/null @@ -1,56 +0,0 @@ -# .replace_prior errors when passed a new prior without a match in old_prior [plain] - - Code - .replace_prior(old_prior, new_prior, warn = TRUE) - Condition - Warning: - ! One or more priors have no match in existing parameters: - Intercept_shape ~ normal(0, 5) - i To remove this warning consider changing prior specification. - Output - prior class coef group resp dpar nlpar lb ub source - normal(0, 5) Intercept user - normal(0, 5) Intercept sigma user - -# .replace_prior errors when passed a new prior without a match in old_prior [ansi] - - Code - .replace_prior(old_prior, new_prior, warn = TRUE) - Condition - Warning: - ! One or more priors have no match in existing parameters: - Intercept_shape ~ normal(0, 5) - i To remove this warning consider changing prior specification. - Output - prior class coef group resp dpar nlpar lb ub source - normal(0, 5) Intercept user - normal(0, 5) Intercept sigma user - -# .replace_prior errors when passed a new prior without a match in old_prior [unicode] - - Code - .replace_prior(old_prior, new_prior, warn = TRUE) - Condition - Warning: - ! One or more priors have no match in existing parameters: - Intercept_shape ~ normal(0, 5) - ℹ To remove this warning consider changing prior specification. - Output - prior class coef group resp dpar nlpar lb ub source - normal(0, 5) Intercept user - normal(0, 5) Intercept sigma user - -# .replace_prior errors when passed a new prior without a match in old_prior [fancy] - - Code - .replace_prior(old_prior, new_prior, warn = TRUE) - Condition - Warning: - ! One or more priors have no match in existing parameters: - Intercept_shape ~ normal(0, 5) - ℹ To remove this warning consider changing prior specification. - Output - prior class coef group resp dpar nlpar lb ub source - normal(0, 5) Intercept user - normal(0, 5) Intercept sigma user - diff --git a/tests/testthat/test-family.R b/tests/testthat/test-family.R index adeb50443..d7bdb8e66 100644 --- a/tests/testthat/test-family.R +++ b/tests/testthat/test-family.R @@ -21,7 +21,7 @@ test_that("the family argument in epidist_family passes as expected for brms and test_that("epidist_family contains the correct reparameterisations for lognormal (no change) and gamma (a change)", { # nolint: line_length_linter. family_lognormal <- epidist_family(prep_obs, family = "lognormal") - expect_identical(family_lognormal$reparam, c("mu", "sigma")) + expect_identical(family_lognormal$param, "mu, sigma") # nolint family_gamma <- epidist_family(prep_obs, family = Gamma(link = "log")) - expect_identical(family_gamma$reparam, c("shape", "shape ./ mu")) # nolint + expect_identical(family_gamma$param, "shape, shape ./ mu") # nolint }) diff --git a/tests/testthat/test-gen.R b/tests/testthat/test-gen.R index 5dc4c5686..37b7d7bf3 100644 --- a/tests/testthat/test-gen.R +++ b/tests/testthat/test-gen.R @@ -115,3 +115,25 @@ test_that("epidist_gen_posterior_epred returns a function that creates arrays wi expect_equal(mean(epred_gamma$.epred), 6.56, tolerance = 0.1) expect_gte(min(epred_gamma$.epred), 0) }) + +test_that("epidist_gen_log_lik returns a function that produces valid log likelihoods", { # nolint: line_length_linter. + skip_on_cran() + # Test lognormal + prep <- brms::prepare_predictions(fit) + prep$ndraws <- 10 + i <- 1 + log_lik_fn <- epidist_gen_log_lik(lognormal()) + log_lik <- log_lik_fn(i = i, prep) + expect_length(log_lik, prep$ndraws) + expect_false(anyNA(log_lik)) + expect_true(all(is.finite(log_lik))) + + # Test gamma + prep_gamma <- brms::prepare_predictions(fit_gamma) + prep$ndraws <- 10 + log_lik_fn_gamma <- epidist_gen_log_lik(Gamma()) + log_lik_gamma <- log_lik_fn_gamma(i = i, prep_gamma) + expect_length(log_lik_gamma, prep_gamma$ndraws) + expect_false(anyNA(log_lik_gamma)) + expect_true(all(is.finite(log_lik_gamma))) +}) diff --git a/tests/testthat/test-latent_model.R b/tests/testthat/test-latent_model.R index c1a64560f..579af69fc 100644 --- a/tests/testthat/test-latent_model.R +++ b/tests/testthat/test-latent_model.R @@ -56,23 +56,3 @@ test_that("epidist_stancode.epidist_latent_model produces valid stanvars", { # n ) expect_s3_class(stancode, "stanvars") }) - -test_that("epidist_gen_log_lik_latent returns a function that produces valid log likelihoods", { # nolint: line_length_linter. - skip_on_cran() - # Test lognormal - prep <- brms::prepare_predictions(fit) - i <- 1 - log_lik_fn <- epidist_gen_log_lik_latent(lognormal()) - log_lik <- log_lik_fn(i = i, prep) - expect_length(log_lik, prep$ndraws) - expect_false(anyNA(log_lik)) - expect_true(all(is.finite(log_lik))) - - # Test gamma - prep_gamma <- brms::prepare_predictions(fit_gamma) - log_lik_fn_gamma <- epidist_gen_log_lik_latent(Gamma()) - log_lik_gamma <- log_lik_fn_gamma(i = i, prep_gamma) - expect_length(log_lik_gamma, prep_gamma$ndraws) - expect_false(anyNA(log_lik_gamma)) - expect_true(all(is.finite(log_lik_gamma))) -}) diff --git a/tests/testthat/test-prior.R b/tests/testthat/test-prior.R index aab1d08f8..053b8d8cb 100644 --- a/tests/testthat/test-prior.R +++ b/tests/testthat/test-prior.R @@ -10,3 +10,61 @@ test_that("epidist_prior with default settings produces an object of the right c expect_s3_class(prior, "brmsprior") expect_s3_class(prior, "data.frame") }) + +test_that("epidist_prior correctly handles user-provided priors", { + data <- as_epidist_latent_model(sim_obs) + family <- lognormal() + formula <- bf(mu ~ 1, sigma ~ 1) + epidist_family <- epidist_family(data, family) + epidist_formula <- epidist_formula( + data = data, family = epidist_family, formula = formula + ) + + user_prior <- prior("normal(0,1)", class = "Intercept") + prior <- epidist_prior( + data, epidist_family, epidist_formula, + prior = user_prior + ) + + expect_identical( + prior$prior[1], + "normal(0,1)" + ) +}) + +test_that("epidist_prior warns about invalid user priors", { + data <- as_epidist_latent_model(sim_obs) + family <- lognormal() + formula <- bf(mu ~ 1, sigma ~ 1) + epidist_family <- epidist_family(data, family) + epidist_formula <- epidist_formula( + data = data, family = epidist_family, formula = formula + ) + + invalid_prior <- prior("normal(0,1)", class = "InvalidClass") + expect_warning( + epidist_prior(data, epidist_family, epidist_formula, prior = invalid_prior), + "One or more priors have no match in existing parameters" + ) +}) + +test_that("epidist_prior correctly applies family-specific priors", { + data <- as_epidist_latent_model(sim_obs) + family <- lognormal() + formula <- bf(mu ~ 1, sigma ~ 1) + epidist_family <- epidist_family(data, family) + epidist_formula <- epidist_formula( + data = data, family = epidist_family, formula = formula + ) + + prior <- epidist_prior(data, epidist_family, epidist_formula, prior = NULL) + + expect_identical( + prior$prior[1], + "normal(1, 1)" + ) + expect_identical( + prior$prior[2], + "normal(-0.7, 0.4)" + ) +}) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index f0012d6ed..1fb3add13 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -9,24 +9,52 @@ test_that(".replace_prior successfully replaces priors", { # nolint: line_length old_prior <- brms::prior("normal(0, 10)", class = "Intercept") + brms::prior("normal(0, 10)", class = "Intercept", dpar = "sigma") new_prior <- brms::prior("normal(0, 5)", class = "Intercept") + - brms::prior("normal(0, 5)", class = "Intercept", dpar = "sigma") + brms::prior("normal(0, 2)", class = "Intercept", dpar = "sigma") prior <- .replace_prior(old_prior, new_prior) - expect_identical(prior$prior, c("normal(0, 5)", "normal(0, 5)")) + expect_identical(prior$prior, c("normal(0, 5)", "normal(0, 2)")) expect_identical(as.double(nrow(prior)), 2) expect_s3_class(prior, "brmsprior") expect_s3_class(prior, "data.frame") }) -cli::test_that_cli(".replace_prior errors when passed a new prior without a match in old_prior", { # nolint: line_length_linter. +cli::test_that_cli(".replace_prior warns when passed a new prior without a match in old_prior", { # nolint: line_length_linter. old_prior <- brms::prior("normal(0, 10)", class = "Intercept") + brms::prior("normal(0, 10)", class = "Intercept", dpar = "sigma") new_prior <- brms::prior("normal(0, 5)", class = "Intercept") + brms::prior("normal(0, 5)", class = "Intercept", dpar = "sigma") + brms::prior("normal(0, 5)", class = "Intercept", dpar = "shape") - expect_snapshot({ - .replace_prior(old_prior, new_prior, warn = TRUE) - }) + expect_warning( + .replace_prior(old_prior, new_prior, warn = TRUE), + "One or more priors have no match in existing parameters" + ) +}) + +test_that(".replace_prior handles custom ~ priors correctly", { + # Create old priors with mix of standard and ~ syntax + old_prior <- brms::prior("mu ~ normal(0, 10)", check = FALSE) + + brms::prior("normal(0, 10)", dpar = "sigma") + + brms::prior("beta ~ normal(0, 1)", check = FALSE) + + # Create new priors with ~ syntax + new_prior <- brms::prior("mu ~ normal(0, 5)", check = FALSE) + + brms::prior("gamma ~ normal(0, 2)", check = FALSE) + + # Test that only old priors with matching ~ parameter names are removed + prior <- .replace_prior(old_prior, new_prior, enforce_presence = FALSE) + + # Should keep sigma prior, replace mu prior, remove beta prior, add gamma + # prior + expect_identical( + prior$prior, + c( + "normal(0, 10)", "mu ~ normal(0, 5)", "gamma ~ normal(0, 2)", + "beta ~ normal(0, 1)" + ) + ) + expect_identical(as.double(nrow(prior)), 4) + expect_s3_class(prior, "brmsprior") + expect_s3_class(prior, "data.frame") }) test_that(".add_dpar_info works as expected for the lognormal and gamma families", { # nolint: line_length_linter. @@ -54,7 +82,8 @@ test_that(".make_intercepts_explicit creates a formula which is the same as if i ) attr(formula$pforms$sigma, ".Environment") <- NULL attr(formula_explicit$pforms$sigma, ".Environment") <- NULL - expect_identical(formula, formula_explicit) + expect_identical(formula$pforms$mu, formula_explicit$pforms$mu) + expect_identical(formula$pforms$sigma, formula_explicit$pforms$sigma) }) test_that(".make_intercepts_explicit does not add an intercept if the distributional parameter is set to be fixed", { # nolint: line_length_linter. @@ -66,7 +95,6 @@ test_that(".make_intercepts_explicit does not add an intercept if the distributi family = epidist_family ) formula_updated <- .make_intercepts_explicit(formula) - attr(formula$pforms$sigma, ".Environment") <- NULL - attr(formula_updated$pforms$sigma, ".Environment") <- NULL - expect_identical(formula, formula_updated) + expect_identical(formula$pforms$mu, formula_updated$pforms$mu) + expect_identical(formula$pforms$sigma, formula_updated$pforms$sigma) }) diff --git a/vignettes/epidist.Rmd b/vignettes/epidist.Rmd index 399652a84..e7c8e807c 100644 --- a/vignettes/epidist.Rmd +++ b/vignettes/epidist.Rmd @@ -173,12 +173,12 @@ obs_cens |> ) + # The primary event censoring intervals geom_errorbarh( - aes(xmin = ptime_lwr, xmax = ptime_upr, y = case, yend = case), + aes(xmin = ptime_lwr, xmax = ptime_upr, y = case), col = "#56B4E9", height = 5 ) + # The secondary event censoring intervals geom_errorbarh( - aes(xmin = stime_lwr, xmax = stime_upr, y = case, yend = case), + aes(xmin = stime_lwr, xmax = stime_upr, y = case), col = "#009E73", height = 5 ) + geom_point(aes(x = ptime), fill = "white", col = "#56B4E9", shape = 21) + @@ -190,7 +190,6 @@ obs_cens |> During an outbreak we will usually be estimating delays in real time. The result is that only those cases with a secondary event occurring before some time will be observed. This is called (right) truncation, and biases the observation process towards shorter delays: - ```{r} obs_time <- 25 diff --git a/vignettes/faq.Rmd b/vignettes/faq.Rmd index 82f028cdb..3343eaf24 100644 --- a/vignettes/faq.Rmd +++ b/vignettes/faq.Rmd @@ -181,7 +181,7 @@ quantile(pred$sd, c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)) We recommend use of the [`priorsense`](https://github.com/n-kall/priorsense) package [@kallioinen2024detecting] to check how sensitive the posterior distribution is to perturbations of the prior distribution and likelihood using power-scaling analysis: -```{r} +```{r, eval = FALSE} library(priorsense) powerscale_plot_dens(fit, variable = c("Intercept", "Intercept_sigma")) + theme_minimal()