Here we extend the previous tutorial to include forecasting capabilities. We start with generating out data using generate_simulation_data_uciwweihr
's alternate parameterization where we do not prespecify the effective reproduction number and hospitalization probability but instead preform a log-normal random walk and a logit-normal random walk respectively. We then sample from the posterior distribution using the uciwweihr_fit.jl
function. We then generate desired quantities and forecast for a given time period with the posterior predictive distribution, using uciwweihr_gq_pp.jl
.
Here we generate two datasets, one with 150 time points and one with 178 time points. We will use the 150 time point dataset for fitting and the 178 time point dataset for forecast evaluation.
using UCIWWEIHR
# Running simulation function with presets
params = create_uciwweihr_sim_params(
time_points = 70
)
df = generate_simulation_data_uciwweihr(params)
params_ext = create_uciwweihr_sim_params(
time_points = 84
)
df_ext = generate_simulation_data_uciwweihr(params_ext)
first(df, 5)
1 | 1 | 0.0882519 | 26 | 0.99777 | 0.349492 | 200.0 | 100.0 | 20.0 |
2 | 2 | 0.407202 | 20 | 1.01468 | 0.353323 | 170.442 | 129.53 | 21.6357 |
3 | 3 | 0.416729 | 26 | 0.989119 | 0.347515 | 150.77 | 149.497 | 24.1411 |
4 | 4 | 0.829472 | 23 | 1.06904 | 0.365335 | 137.051 | 162.981 | 26.8269 |
5 | 5 | 0.75552 | 31 | 1.10915 | 0.373918 | 129.444 | 172.241 | 29.9119 |
first(df_ext, 5)
1 | 1 | 0.0882519 | 17 | 0.99777 | 0.349492 | 200.0 | 100.0 | 20.0 |
2 | 2 | 0.407202 | 32 | 1.01468 | 0.353323 | 170.442 | 129.53 | 21.6357 |
3 | 3 | 0.416729 | 18 | 0.989119 | 0.347515 | 150.77 | 149.497 | 24.1411 |
4 | 4 | 0.829472 | 28 | 1.06904 | 0.365335 | 137.051 | 162.981 | 26.8269 |
5 | 5 | 0.75552 | 26 | 1.10915 | 0.373918 | 129.444 | 172.241 | 29.9119 |
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. The difference here is that we set forecast = true
and forecast_weeks = 4
to forecast 4 weeks into the future. One other thing to note, is that we 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 = true
forecast_days = 14
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.212128 | 0.525665 | 0.438046 |
2 | 2 | 1 | -0.0806217 | 0.424289 | 0.48069 |
3 | 3 | 1 | 0.0746061 | 0.216213 | 0.411932 |
4 | 4 | 1 | 0.257695 | 0.291986 | 0.525137 |
5 | 5 | 1 | 0.243721 | 0.331795 | 0.495088 |
first(model_output[2][:,1:5], 5)
1 | 1 | 1 | 195.628 | 74.0774 | 19.987 |
2 | 2 | 1 | 195.628 | 74.0774 | 19.987 |
3 | 3 | 1 | 195.628 | 74.0774 | 19.987 |
4 | 4 | 1 | 222.705 | 70.0517 | 22.2225 |
5 | 5 | 1 | 182.93 | 91.2814 | 17.9419 |
first(model_output[3][:,1:5], 5)
1 | 201 | 1 | -0.110502 | -1.5003 | -0.00326038 |
2 | 202 | 1 | -0.110502 | -1.5003 | -0.00326038 |
3 | 203 | 1 | -0.110502 | -1.5003 | -0.00326038 |
4 | 204 | 1 | 0.537644 | -1.77968 | 0.526865 |
5 | 205 | 1 | -0.446067 | -0.456116 | -0.542966 |
We can again look at model diagnostics, posterior distribution of time or non-time varying parameters, and the posterior predictive distribution extended for forecasting. 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 = df_ext.hosp,
obs_data_wastewater = df_ext.log_ww_conc,
actual_rt_vals = df_ext.rt,
actual_w_t = df_ext.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,
plot_name_to_save_mcmcdiag = "mcmc_diagnosis_plots1",
plot_name_to_save_time_varying = "mcmc_time_varying_parameter_plots1",
plot_name_to_save_non_time_varying = "mcmc_nontime_varying_parameter_plots1",
plot_name_to_save_ode_sol = "mcmc_ode_solution_plots1",
plot_name_to_save_pred_param = "mcmc_pred_parameter_plots1"
)
Effective Sample Size for E_init for Chain 1: 181.0
Effective Sample Size for I_init for Chain 1: 223.0
Effective Sample Size for H_init for Chain 1: 708.0
Effective Sample Size for gamma for Chain 1: 602.0
Effective Sample Size for nu for Chain 1: 649.0
Effective Sample Size for epsilon for Chain 1: 453.0
Effective Sample Size for rt_init for Chain 1: 188.0
Effective Sample Size for w_init for Chain 1: 444.0
Effective Sample Size for sigma_w for Chain 1: 564.0
Effective Sample Size for sigma_Rt for Chain 1: 417.0
Effective Sample Size for sigma_ww for Chain 1: 528.0
Effective Sample Size for sigma_hosp for Chain 1: 734.0
Effective Sample Size for rho_gene for Chain 1: 490.0
Plot saved to plots/mcmc_diagnosis_plots1.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:04
Sampling (1 threads) 100%|██████████████████████████████| Time: 0:00:04
Using uciwweihr_model with wastewater. Priors on sigma_ww and sigma_hosp!!!
Generating quantities...
Plot saved to plots/mcmc_time_varying_parameter_plots1.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_plots1.png
Plot saved to plots/mcmc_ode_solution_plots1.png
Plot saved to plots/mcmc_pred_parameter_plots1.png




