Generating Posterior Distribution Samples with UCIWWEIHR ODE Compartmental Based Model with Forecasting.

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.

1. Data Generation.

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)
5×8 DataFrame
Rowobstimeslog_ww_conchosprtwtE_ode_comp_solI_ode_comp_solH_ode_comp_sol
Int64Float64Int64Float64Float64Float64Float64Float64
110.0882519260.997770.349492200.0100.020.0
220.407202201.014680.353323170.442129.5321.6357
330.416729260.9891190.347515150.77149.49724.1411
440.829472231.069040.365335137.051162.98126.8269
550.75552311.109150.373918129.444172.24129.9119
first(df_ext, 5)
5×8 DataFrame
Rowobstimeslog_ww_conchosprtwtE_ode_comp_solI_ode_comp_solH_ode_comp_sol
Int64Float64Int64Float64Float64Float64Float64Float64
110.0882519170.997770.349492200.0100.020.0
220.407202321.014680.353323170.442129.5321.6357
330.416729180.9891190.347515150.77149.49724.1411
440.829472281.069040.365335137.051162.98126.8269
550.75552261.109150.373918129.444172.24129.9119

2. Sampling from the Posterior Distribution and Posterior Predictive Distribution.

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)
5×5 DataFrame
Rowiterationchaindata_wastewater[1]data_wastewater[2]data_wastewater[3]
Int64Int64Float64Float64Float64
1110.2121280.5256650.438046
221-0.08062170.4242890.48069
3310.07460610.2162130.411932
4410.2576950.2919860.525137
5510.2437210.3317950.495088
first(model_output[2][:,1:5], 5)
5×5 DataFrame
RowiterationchainE_initI_initH_init
Int64Int64Float64Float64Float64
111195.62874.077419.987
221195.62874.077419.987
331195.62874.077419.987
441222.70570.051722.2225
551182.9391.281417.9419
first(model_output[3][:,1:5], 5)
5×5 DataFrame
RowiterationchainE_init_non_centeredI_init_non_centeredH_init_non_centered
Int64Int64Float64Float64Float64
12011-0.110502-1.5003-0.00326038
22021-0.110502-1.5003-0.00326038
32031-0.110502-1.5003-0.00326038
420410.537644-1.779680.526865
52051-0.446067-0.456116-0.542966

3. MCMC Diagnostic Plots/Results Along with Posterior Predictive Distribution.

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

3.1. MCMC Diagnostic Plots.

Plot 1

3.2. Time Varying Parameter Results Plot.

Plot 2

3.3. Non-Time Varying Parameter Results Plot.

Plot 3

3.4. ODE Solution Plot.

Plot 4

3.5. Posterior Predictive Distribution Plot.

Plot 4

Tutorial Contents