This package has a way to sample from a posterior or prior that is defined in the future paper using the uciwweihr_fit.jl
and uciwweihr_model.jl
. We can then generate desired quantities and forecast for a given time period with the posterior predictive distribution, using uciwweihr_gq_pp.jl
. We first generate data using the generate_simulation_data_uciwweihr
function which is a non-mispecified version of the model, we will also be using prespecified effective reporduction curves and prespecified hospitalization probability curves.
using UCIWWEIHR
# Running simulation function with presets
rt_custom = vcat(
range(1, stop=1.8, length=7*2),
fill(1.8, 7*1),
range(1.8, stop=1, length=7*4),
range(0.98, stop=0.8, length=7*1),
range(0.8, stop=1.1, length=7*3),
range(1.1, stop=0.97, length=7*1)
)
w_custom = vcat(
range(0.3, stop=0.38, length=7*2),
fill(0.38, 7*1),
range(0.38, stop=0.25, length=7*4),
range(0.25, stop=0.28, length=7*1),
range(0.28, stop=0.34, length=7*3),
range(0.34, stop=0.28, length=7*1)
)
params = create_uciwweihr_sim_params(
time_points = length(rt_custom),
Rt = rt_custom,
w = w_custom
)
df = generate_simulation_data_uciwweihr(params)
first(df, 5)
1 | 1 | 0.0882519 | 17 | 1.0 | 0.3 | 200.0 | 100.0 | 20.0 |
2 | 2 | 0.407232 | 31 | 1.06154 | 0.306154 | 170.481 | 129.534 | 20.8906 |
3 | 3 | 0.417484 | 17 | 1.12308 | 0.312308 | 151.641 | 149.61 | 22.6748 |
4 | 4 | 0.833164 | 26 | 1.18462 | 0.318462 | 140.427 | 163.584 | 24.9225 |
5 | 5 | 0.764396 | 24 | 1.24615 | 0.324615 | 134.695 | 173.776 | 27.3775 |
Here we sample from the posterior distribution using the uciwweihr_fit.jl
function. First, we setup some presets, where we need to use create_uciwweihr_model_params()
to get default parameters for the model. Then we have an array where index 1 contains the posterior/prior predictive samples, index 2 contains the posterior/prior generated quantities samples, and index 3 contains the original sampled parameters for the model. Again, we can allow misalignment of hospital and wastewater data's observed times. For this tutorial, we use the same observed points.
data_hosp = df.hosp
data_wastewater = df.log_ww_conc
obstimes_hosp = df.obstimes
obstimes_wastewater = df.obstimes
max_obstime = max(length(obstimes_hosp), length(obstimes_wastewater))
param_change_times = 1:7:max_obstime # Change every week
priors_only = false
n_samples = 500
forecast = false
forecast_days = 0
model_params = create_uciwweihr_model_params2()
samples = uciwweihr_fit(
data_hosp,
data_wastewater,
obstimes_hosp,
obstimes_wastewater,
param_change_times,
model_params;
priors_only,
n_samples,
n_discard_initial = 200
)
model_output = uciwweihr_gq_pp(
samples,
data_hosp,
data_wastewater,
obstimes_hosp,
obstimes_wastewater,
param_change_times,
model_params;
forecast=forecast, forecast_days=forecast_days
)
first(model_output[1][:,1:5], 5)
1 | 1 | 1 | 0.215332 | 0.447321 | 0.470766 |
2 | 2 | 1 | 0.109267 | 0.433796 | 0.512615 |
3 | 3 | 1 | 0.275915 | 0.350607 | 0.53743 |
4 | 4 | 1 | 0.295538 | 0.328554 | 0.552504 |
5 | 5 | 1 | 0.191069 | 0.30882 | 0.522431 |
first(model_output[2][:,1:5], 5)
1 | 1 | 1 | 197.921 | 78.6483 | 19.7918 |
2 | 2 | 1 | 197.921 | 78.6483 | 19.7918 |
3 | 3 | 1 | 197.921 | 78.6483 | 19.7918 |
4 | 4 | 1 | 164.084 | 77.9884 | 16.6532 |
5 | 5 | 1 | 177.662 | 79.2565 | 21.1659 |
first(model_output[3][:,1:5], 5)
1 | 201 | 1 | -0.0522502 | -1.20092 | -0.0523232 |
2 | 202 | 1 | -0.0522502 | -1.20092 | -0.0523232 |
3 | 203 | 1 | -0.0522502 | -1.20092 | -0.0523232 |
4 | 204 | 1 | -0.989687 | -1.24305 | -0.915657 |
5 | 205 | 1 | -0.59217 | -1.16241 | 0.283301 |
We also provide a very basic way to visualize some MCMC diagnostics along with effective sample sizes of desired generated quantities(does not include functionality for time-varying quantities). Along with this, we can also visualize the posterior predictive distribution with actual observed values, which can be used to examine forecasts generated by the model. We can also add certain parameters to ensure priors will be plotted alongside their corresponding posteriors.
uciwweihr_visualizer(
data_hosp,
data_wastewater,
forecast_days,
obstimes_hosp,
obstimes_wastewater,
param_change_times,
2024,
forecast,
model_params;
pp_samples = model_output[1],
gq_samples = model_output[2],
obs_data_hosp = data_hosp,
obs_data_wastewater = data_wastewater,
actual_rt_vals = df.rt,
actual_w_t = df.wt,
actual_E_ode_sol = df.E_ode_comp_sol,
actual_I_ode_sol = df.I_ode_comp_sol,
actual_H_ode_sol = df.H_ode_comp_sol,
actual_non_time_varying_vals = params,
bayes_dist_type = "Posterior",
save_plots = true
)
Effective Sample Size for E_init for Chain 1: 350.0
Effective Sample Size for I_init for Chain 1: 439.0
Effective Sample Size for H_init for Chain 1: 769.0
Effective Sample Size for gamma for Chain 1: 555.0
Effective Sample Size for nu for Chain 1: 797.0
Effective Sample Size for epsilon for Chain 1: 519.0
Effective Sample Size for rt_init for Chain 1: 299.0
Effective Sample Size for w_init for Chain 1: 797.0
Effective Sample Size for sigma_w for Chain 1: 754.0
Effective Sample Size for sigma_Rt for Chain 1: 553.0
Effective Sample Size for sigma_ww for Chain 1: 824.0
Effective Sample Size for sigma_hosp for Chain 1: 674.0
Effective Sample Size for rho_gene for Chain 1: 867.0
Plot saved to plots/mcmc_diagnosis_plots.png
Generating time varying parameter plots (with priors and with wastewater)...
Using uciwweihr_model with wastewater. Priors on sigma_ww and sigma_hosp!!!
┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/FSyVk/src/sample.jl:382
Sampling (1 threads) 0%| | ETA: N/A
Sampling (1 threads) 100%|██████████████████████████████| Time: 0:00:00
Sampling (1 threads) 100%|██████████████████████████████| Time: 0:00:00
Using uciwweihr_model with wastewater. Priors on sigma_ww and sigma_hosp!!!
Generating quantities...
Plot saved to plots/mcmc_time_varying_parameter_plots.png
Generating non-time varying parameter plots (with priors and with wastewater)...
Using uciwweihr_model with wastewater. Priors on sigma_ww and sigma_hosp!!!
┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/FSyVk/src/sample.jl:382
Sampling (1 threads) 0%| | ETA: N/A
Sampling (1 threads) 100%|██████████████████████████████| Time: 0:00:00
Sampling (1 threads) 100%|██████████████████████████████| Time: 0:00:00
Using uciwweihr_model with wastewater. Priors on sigma_ww and sigma_hosp!!!
Generating quantities...
Plot saved to plots/mcmc_nontime_varying_parameter_plots.png
Plot saved to plots/mcmc_ode_solution_plots.png
Plot saved to plots/mcmc_pred_parameter_plots.png




