Goal

Want to set a threshold for the maximum allowable rhat in the model_diagnostics() function

Context

#5 (comment)

Goal

Set up caching of vignette build. Right now, this takes a long time because it has to run the full example model.

Requirements

• if nothing that impacts the vignette gets updated, we want to cache the vignette build

Goal

Currently one long linear function that steps through each component of the model to generate simulated hospital admissions and wastewater data. It is highly complex and hard to follow. It would be easier if each individual module called a function. That way, when you want to change how a module is handled, you can just write a new function for it.

Requirements

• break down getting global and local infection dynamics, generating expected hospital admissions, and generating expected wastewater concentrations into separate modular functions
• call these in the generate_simulated_data.R function

Goal

Explore how this might work for the aspatial model, might improve identifiability

Context

Suggestion from @athowes to improve identifiability, in particular might be helpful when there are deviations that could be due to site/subpopulation level changes and also could be due to deviation in "central/global" $R(t)$
Suggests something like:
sum(subpop_level_errors[t]) ~ normal(0, <small_number> *<n_subpops>)

Requirements

• Rmd to explore implementation of different model components and compare outputs

Out of scope for this PR

• Idea of this PR is just to see how this impacts outputs/diagnostics/etc. Before adding this component to the model we would want to do some evaluation to ensure it improved performance (e.g. what we have in the pipeline repo for a single forecast date)

Goal

Would be nice to have a little model diagram. Draft is here, but we need to figure out how to show that the infections from each subpopulation sum to the state infections.

The general flow of the model can be visualized as follows:

flowchart LR;
    GlobalR --> GlobalI;
    GlobalR --> LocalR;
    LocalR --> GlobalR;
    LocalR --> LocalIX;
    LocalR --> LocalIY;
    LocalIX --> LocalWWX;
    LocalIY --> LocalWWY;
    GlobalI --> Globalhosp;

Sounds good. Maybe rename to format_data_for_stan or something? Separate issue.

Originally posted by @dylanhmorris in #5 (comment)

Would be good to make the small added on value (here 1e-8) configurable, but that can be a separate PR.

Originally posted by @dylanhmorris in #11 (comment)

This is still needed, I tried running it without it and ran into overflow issues.

I know that its pretty hard-coded as is, would you suggest we pass in the sds that we sample the from for each parameter?

Originally posted by @kaitejohnson in #5 (comment)

Let's make an issue to discuss how to make more generic.

Originally posted by @dylanhmorris in #5 (comment)

Goal

We don't need both pkg::fun() and importFrom(pkg, fn). Go through and streamline this, adding things to wwinference-package.R when its clear what package they would belong to, and then remove from the codebase.

Context

#5 (comment)

Goal

We want to show an example of how someone would use this package in a production pipeline.

Requirements

• a script that can be called from the command line with parsed args getting passed
• an example config file containing those parsed args
• a pre-processing script/set of functions from a common input dataset example
• a post-processing script/set of functions

Sure, how would you suggest fixing it? I agree its clunky. Could just unlist a list? Would prefer to not have to always call things by params$parameter_name, unless that is the best practice.

Originally posted by @kaitejohnson in #5 (comment)

Goal

Want package documentation to create a link to other function documentation

Context

#5 (comment)
See example above

Goal

Want to make sure that the input data that a user provides for hospital admissions and wastewater datasets meets requirements needed for the model functions to work correctly. If not, we want to provide them with coherent error messaging to tell them how to modify it.

Requirements

• model some checks after @zsusswein 's checkers.R https://github.com/CDCgov/cfa-gam-rt/blob/main/R/checkers.R and other validations to ensure that input data is formatted properly.

Goal

We want to set up CI s.t. all functions tests pass before merging to main. We have some of this already set up from the https://github.com/CDCgov/wastewater-informed-covid-forecasting/tree/prod so we will rely on that infrastructure to set this up.

Requirements

• a testing.yaml that uses testthat() in github actions to check the package functions

Goal

It is useful to compare the forecasts from the same model with and without the wastewater component included. Ideally, this doesn't require using the separate stan file and stan data object, but just turning off the wastewater component on the likelihood and making sure that the IHR is not time-varying

Requirements

• implement a hospital admissions only renewal fit
• demonstrate the fit to only hospital admissions data in the getting started vignette

Goal

Streamline post-processing function so that they call other functions that generate the data-spines that we will pass to the user and use to create the draws dataframe

Context

#5 (comment)

Requirements

• separate functions for each of the returned data spines in postprocess.R

Goal

We want to create STAN functions that are similar to the R functions for the spatial components of the model.

Context

The following model will be added as STAN functions :
Let $k$ represent a site, from sites $1$ to $n_k$. Let $R_k^u(t)$ be the unadjusted reproduction number at time $t$ at site $k$. So lets look at the combined "state" $R^u(t)$ that contains an AR(1) correlation structure on the first differences that is currently being used. So

$$ \log R^u(t) \sim \text{Normal}(\log R^u(t - 1) + \beta(\log R^u(t - 1) - \log R^u(t - 2)), \sigma_r), $$

or

$$ \log R^u(t) = \log R^u(t - 1) + \beta(\log R^u(t - 1) - \log R^u(t - 2)) + \sigma_r\text{Normal}(0,1). $$

Now let $\delta_k(t)$ be the time-varying subpopulation effect on $\log R^u(t)$. So we can look at the vector of all $\log R_k^u(t)$'s, so

$$ \begin{pmatrix} \log R_1^u(t)\\ \vdots\\ \log R_{n_{k}}^u(t)\\ \end{pmatrix} = \log R^u(t) \begin{pmatrix} 1\\ \vdots\\ 1\\ \end{pmatrix} + \vec{\delta}(t). $$

Where $\vec{\delta}(t) = \left(\delta_1(t),...,\delta_{n_k}(t)\right)^{T}$. We also have that $\vec{\delta}(t) = \varphi_{R(t)}\vec{\delta}(t - 1) + \vec{\epsilon}$, where $\vec{\epsilon}$ $\sim$ $\text{Normal}_{n_k}(0,\mathbf{\Sigma_{\epsilon}}).$

This adds a spatial correlation and also the current AR(1) temporal correlation to this $\vec{\delta}(\cdot)$ random variable. Like before, we have that $0 < \varphi_{R(t)} < 1$.

For $\mathbf{\Sigma_{\epsilon}}$ $=$ $\sigma_{\epsilon}^2\mathbf{\Omega}$, where $\Omega_{ij} \in \mathbf{\Omega}$ are build from the correlation functions. These include an exponential decay, Matern, spherical, rational quadratic, or, to return to the current model framework, independence correlation functions.

For simplicity, we will first assume an distance matrix is provided with the data or locations of sites where we can attain an distance matrix.

We need to make sure to take into account is the individuals that are not included in a site. So, currently, we will provide them with their own process. This separate process is as follows :

Lets denote the unadjusted reproduction number at time $t$ for individuals who are not represented by the actual sites as $R_{\text{aux}}^u(t)$.

$\log R_{\text{aux}}^u(t) = \log R^u(t) + \delta_{\text{aux}}(t),$

$\delta_{\text{aux}}(t) = \phi_{R(t)}\delta_{\text{aux}}(t - 1) + \epsilon_{\text{aux}},$

$\epsilon_{\text{aux}} \sim \text{Normal}(0,\text{ScalingFactor}\times\sigma_{\epsilon}),$

This process is very similar with the previous way site level Rt was calculated, but here we include the ScalingFactor to scale the variation of the spatial process. Also, if we have an independent structure, we will not have this ScalingFactor and will instead use the $n_k + 1$ approach currently being implemented.

Requirements

• Add spatial_rt_process.stan
• Add independence_corr_func.stan
• Add exponential_corr_func.stan
• Add aux_site_process.stan

Oh that's a great question... I can't remember why this was needed tbh, let me test out it tomorrow. It was something weird that might have had to do with the VAP, and this ensured permissions were available to create a folder...

Originally posted by @kaitejohnson in #5 (comment)

Goal

We want to use the spatial process that was recently added into the spatial-ww branch to update the package data that is used in the vignette.

Requirements

• Update data using the vignette_data.R file

Goal

We want to update the generate_simulated_data.R function in the R folder to implement a spatially correlated unadjusted reproduction number into it. We will also add necessary functions to allow different correlation structures to be used.

Context

The current model used in the generate_simulated_data.R function is listed in the model definitions in the repo. The change we are making will be at the unadjusted reproduction level. We will add the spatial correlations here with the following model :

Let $k$ represent a site, from sites $1$ to $n_k$. Let $R_k^u(t)$ be the unadjusted reproduction number at time $t$ at site $k$. So lets look at the combined "state" $R^u(t)$ that contains an AR(1) correlation structure on the first differences that is currently being used. So

$$ \log R^u(t) \sim \text{Normal}(\log R^u(t - 1) + \beta(\log R^u(t - 1) - \log R^u(t - 2)), \sigma_r), $$

or

$$ \log R^u(t) = \log R^u(t - 1) + \beta(\log R^u(t - 1) - \log R^u(t - 2)) + \sigma_r\text{Normal}(0,1). $$

Now let $\delta_k(t)$ be the time-varying subpopulation effect on $\log R^u(t)$. So we can look at the vector of all $\log R_k^u(t)$'s, so

$$ \begin{pmatrix} \log R_1^u(t)\\ \vdots\\ \log R_{n_{k}}^u(t)\\ \end{pmatrix} = \log R^u(t) \begin{pmatrix} 1\\ \vdots\\ 1\\ \end{pmatrix} + \vec{\delta}(t). $$

Where $\vec{\delta}(t) = \left(\delta_1(t),...,\delta_{n_k}(t)\right)^{T}$. We also have that $\vec{\delta}(t) = \varphi_{R(t)}\vec{\delta}(t - 1) + \vec{\epsilon}$, where $\vec{\epsilon}$ $\sim$ $\text{Normal}_{n_k}(0,\mathbf{\Sigma_{\epsilon}}).$

This adds a spatial correlation and also the current AR(1) temporal correlation to this $\vec{\delta}(\cdot)$ random variable. Like before, we have that $0 < \varphi_{R(t)} < 1$.

For $\mathbf{\Sigma_{\epsilon}}$ $=$ $\sigma_{\epsilon}^2\mathbf{\Omega}$, where $\Omega_{ij} \in \mathbf{\Omega}$ are build from the correlation functions. These include an exponential decay, Matern, spherical, rational quadratic, or, to return to the current model framework, independence correlation functions.

For simplicity, we will first assume an adjacency matrix is provided with the data or locations of sites where we can attain an adjacency matrix.

We need to make sure to take into account is the individuals that are not included in a site. So, currently, we will provide them with their own process. This separate process is as follows :

Lets denote the unadjusted reproduction number at time $t$ for individuals who are not represented by the actual sites as $R_{\text{aux}}^u(t)$.

$\log R_{\text{aux}}^u(t) = \log R^u(t) + \delta_{\text{aux}}(t),$

$\delta_{\text{aux}}(t) = \phi_{R(t)}\delta_{\text{aux}}(t - 1) + \epsilon_{\text{aux}},$

$\epsilon_{\text{aux}} \sim \text{Normal}(0,\text{ScalingFactor}\times\sigma_{\epsilon}),$

$\text{ScalingFactor} \sim \text{Log-Normal}(\mu_{\text{ScalingFactor}},\sigma_{\text{ScalingFactor}}^2).$

This process is very similar with the previous way site level Rt was calculated, but here we include the ScalingFactor to scale the variation of the spatial process. Also, if we have an independent structure, we will not have this ScalingFactor and will instead use the $n_k + 1$ approach currently being implemented.

Requirements

• Add function, named spatial_rt_process.R, to deal with the new spatial unadjusted reproduction number for a general correlation function.
• Change generate_simulated_data.R for implementation of spatial unadjusted reproduction number, but will default to the independent unadjusted reproduction number.
• Add R functions that create correlation matrices for the spatial_rt_process.R function.

Goal

We created a bunch of new functions, and we should add a testing suite for them. The functions that need to be tested include:
preprocess_ww_data preprocess_hosp_data flag_ww_outliers, indicate_ww_exclusions in preprocessing.R get_draws_df in get_draws_df.R and the rest of the functions in utils.R.

Requirements

• additional tests in tests folder for each of these functions

Goal

We want to have a full testing suite for all of the functions in the wwinference package. Some of these, particularly for the stan model, are already written in https://github.com/CDCgov/wastewater-informed-covid-forecasting/tree/prod, and so we will want to port them over, and we also want to make new tests.

Requirements:

• move over/modify if needed all existing tests
• add tests for new functions e.g. pre and post processing --> made a new issue for this #27

Goal

Tracking steps to a MVP modeling package. The goal is a minimum package that we can share with internal and external partners wanting to run the model on their own data and/or set up production pipelines.

Features

• R package infrastructure #1
• CI set up for testing suite
• Minimum pre-processing and model fitting wrapper in a "getting started" vignette #5
• Validate input data #4
• Port over existing model tests
• Add new tests for functions #27
• Make sure every exported function has complete documentation
• #3
• Minimum post-processing to make it easy for the user to access: posterior draws of global $R(t)$, site-level $R(t)$, hospital admissions calibrated and forecasted, wastewater concentrations calibrated and forecasted #5
• Updated README #7
• Updated model definition #7

Goal

Make the function more general, add the current default arguments in to the vignette_data.R to show how we generated the package data

Context

#5 (comment)