class because the object no longer
## meets its requirements:
## ✖ Assertion on 'names(data)' failed: Names must include the elements
## {'delay_lwr','delay_upr','delay_min','relative_obs_time','pwindow','swindow','n'},
## but is missing elements {'n'}.
## ℹ Use the matching `as_epidist_*()` function to recreate the object.
```
``` r
# Run pp_check with the expanded data
pp_check(fit, newdata = data_expanded, ndraws = 100)
```
For more advanced custom visualizations, you can also use `tidybayes::add_predicted_draws()` as demonstrated in the ["Advanced features with Ebola data"](https://epidist.epinowcast.org/articles/ebola.html) vignette.
## What is an observation in the meta model?
`as_epidist_meta_model()` fits the summaries a study computed from the same delays jointly, so one row of the model data is a group of summaries rather than a single reported value.
This is the same situation as aggregated data, where one row stands for many individual delays.
`brms::log_lik()` and `loo::loo()` work at the row level, so their values are per group, and `pp_check()` needs the same care as described above.
Predictions for summary rows are the reported summary itself rather than a delay, so restrict `pp_check()` to the individual level rows of a mixed fit.
## I'd like to run a simulation study
We recommend use of the [`purrr`](https://purrr.tidyverse.org/) package for running many `epidist` models, for example as a part of a simulation study.
We particularly highlight two functions which might be useful:
1. `purrr::map()` (and other similar functions) for iterating over a list of inputs.
2. `purrr::safely()` which ensures that the function called "always succeeds". In other words, if there is an error it will be captured and output, rather than ending computation (and potentially disrupting a call to `purrr::map()`).
For an example use of these functions, have a look at the [`epidist-paper`](https://github.com/parksw3/epidist-paper) repository containing the code for @park2024estimating.
(Note that in that codebase, we use `map` as a part of a [`targets`](https://books.ropensci.org/targets/) pipeline.)
## How are the default priors for `epidist` chosen?
[`brms`](https://paulbuerkner.com/brms/) provides default priors for all parameters.
However, some of those priors do not make sense in the context of our application.
Instead, we used [prior predictive checking](https://mc-stan.org/docs/stan-users-guide/posterior-predictive-checks.html) to set `epidist`-specific default priors which produce epidemiological delay distribution mean and standard deviation parameters in a reasonable range.
For example, for the `brms::lognormal()` latent individual model, we suggest the following prior distributions for the `brms` `mu` and `sigma` intercept parameters:
``` r
# Note that we export lognormal() as part of epidist hence no need for brms::
family <- lognormal()
epidist_family <- epidist_family(data, family)
epidist_formula <- epidist_formula(
data,
family = epidist_family,
formula = mu ~ 1
)
# NULL here means no replacing priors from the user!
epidist_prior <- epidist_prior(
data = data,
family = family,
formula = epidist_formula,
prior = NULL
)
epidist_prior
```
```
## prior class coef group resp dpar nlpar lb ub tag source
## student_t(3, 5, 3) Intercept default
## student_t(3, 0, 2.5) Intercept sigma default
```
(Note that the functions `epidist_family()` and `epidist_prior()` are mostly for internal use!)
Here are the distributions on the delay distribution mean and standard deviation parameters that these prior distributions imply:
``` r
set.seed(1)
fit_ppc <- epidist(
data = data,
formula = mu ~ 1,
family = lognormal(),
sample_prior = "only",
seed = 1,
backend = "cmdstanr"
)
```
``` r
pred <- fit_ppc |>
epidist_strata() |>
add_delay_parameter_draws(fit_ppc) |>
add_summaries()
plot(pred, pars = c("mean", "sd")) +
scale_x_log10(labels = comma)
```
``` r
quantile(pred$mean, c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99))
```
```
## 1% 10% 25% 50% 75% 90% 99%
## 0.3172760 0.8667011 1.6163277 3.2461963 6.6291746 12.0852973 34.0415723
```
``` r
quantile(pred$sd, c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99))
```
```
## 1% 10% 25% 50% 75% 90% 99%
## 0.1335500 0.3950604 0.7687361 1.7107506 3.9206460 8.1670947 33.4423697
```
## How can I assess how sensitive the fitted posterior distribution is to the prior distribution used?
We recommend use of the [`priorsense`](https://github.com/stan-dev/priorsense) package [@kallioinen2024detecting] to check how sensitive the posterior distribution is to perturbations of the prior distribution and likelihood using power-scaling analysis:
``` r
library(priorsense)
powerscale_plot_dens(fit, variable = c("Intercept", "Intercept_sigma")) +
theme_minimal()
```
# What do the parameters in my model output correspond to?
The `epidist` package uses [`brms`](https://paulbuerkner.com/brms/) to fit models.
This means that the model output will include `brms`-style names for parameters.
Here, we provide a table giving the correspondence between the distributional parameter names used in `brms` and those used in standard R functions for some common likelihood families.
| Family | `brms` parameter | R parameter |
|---------------|------------------|-------------|
| `lognormal()` | `mu` | `meanlog` |
| `lognormal()` | `sigma` | `sdlog` |
| `Gamma()` | `mu` | `shape / scale` |
| `Gamma()` | `shape` | `shape` |
| `weibull()` | `mu` | `scale * gamma(1 + 1 / shape)` |
| `weibull()` | `shape` | `shape` |
| `gengamma()` | `mu` | `mu` of `flexsurv::dgengamma()` |
| `gengamma()` | `sigma` | `sigma` |
| `gengamma()` | `Q` | `Q` |
Note that all families in `brms` are parameterised with some measure of centrality `mu` as their first parameter.
This parameter does not necessarily correspond to the mean, hence the provision of a function `add_summaries()` within `epidist` to add columns containing the natural scale mean, standard deviation and quantiles to a `data.frame` of draws.
# How can I generate predictions with my fitted `epidist` model?
It is possible to generate predictions manually by working with [samples from the model output](https://epidist.epinowcast.org/articles/faq.html#i-would-like-to-work-with-the-samples-output).
However this is tricky to do, and so where possible we recommend using the [`tidybayes`](https://mjskay.github.io/tidybayes/) package.
In particular, following functions may be useful:
1. `tidybayes::add_epred_draws()` for predictions of the expected value of a delay.
2. `tidybayes::add_linpred_draws()` for predictions of the delay distributional parameter linear predictors.
3. `tidybayes::add_predicted_draws()` for predictions of the observed delay.
To see these functions demonstrated in a vignette, see ["Advanced features with Ebola data"](https://epidist.epinowcast.org/articles/ebola.html).
As a short example, to generate 4000 predictions (equal to the number of draws) of the delay that would be observed with a double censored observation process (in which the primary and secondary censoring windows are both one) then:
``` r
library(tidybayes)
draws_pmf <- tibble::tibble(
relative_obs_time = Inf, pwindow = 1, swindow = 1, delay_upr = NA,
delay_min = 0
) |>
add_predicted_draws(fit)
```
```
## Warning: Found infinite values in the data, which may cause issues for Stan.
```
``` r
ggplot(draws_pmf, aes(x = .prediction)) +
geom_bar(aes(y = after_stat(count / sum(count)))) +
labs(x = "Delay", y = "PMF") +
scale_x_continuous(limits = c(0, 30)) +
theme_minimal()
```
```
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_bar()`).
```
Importantly, this functionality is only available for `epidist` models using `brms` families that have a `log_lik_censor` method implemented internally in `brms`. If you are using another family, consider [submitting a pull request](https://github.com/epinowcast/epidist/pulls) to implement these methods!
## How can I use the `marginaleffects` package with my fitted `epidist` model?
The [`marginaleffects`](https://marginaleffects.com/) package provides tools for computing and visualising marginal effects, contrasts, and predictions from regression models. It works with `epidist` models because they are built on top of `brms`.
For `epidist` models with covariates, you can use `marginaleffects` to:
- Compute average marginal effects
- Make comparisons between different covariate values
- Visualise model predictions across the range of covariates
Here's a simple example using a model that includes location as a covariate:
``` r
library(marginaleffects)
# Compare outcomes between location categories (0 and 1)
avg_comparisons(
fit_location,
variables = list(location = c(0, 1))
)
```
```
## Warning: The training data could not be extracted reliably from the model object. Column types, factor levels, or values may have been coerced or altered. It is safer to use `set_modeldata(model, data)` to attach the training data explicitly. Original warning from `insight::get_data()`: Could not recover model data from environment. Please make sure your
## data is available in your workspace.
## Trying to retrieve data from the model frame now. Original warning from `insight::get_data()`: Following potential variables could not be found in the data: n) +
## vreal(relative_obs_time, pwindow, swindow, delay_upr, delay_min This warning appears once per session.
```
```
##
## Estimate 2.5 % 97.5 %
## 6.02 2.12 29.8
##
## Term: location
## Type: response
## Comparison: 1 - 0
```
For more details on the available functions and their use, see the [`marginaleffects` documentation](https://marginaleffects.com/).
## How can I use the `cmdstanr` backend?
The `cmdstanr` backend is typically more performant than the default `rstan` backend.
To use the `cmdstanr` backend, we first need to install CmdStan (see the README for more details). We can check we have everything we need as follows:
``` r
cmdstanr::cmdstan_version()
```
```
## [1] "2.40.0"
```
## References
