+
+
+
+
+
+
+
+
+
+
+Getting started with EpiAware
This tutorial introduces the basic functionality of EpiAware. EpiAware is a package for making inferences on epidemiological case/determined infection data using a model-based approach.
It is common to conceptualise the generative process of public health data, e.g a time series of reported cases of an infectious pathogen, in a modular way. For example, it is common to abstract the underlying latent infection process away from downstream issues of observation, or to treat quanitites such as the time-varying reproduction number as being itself generated as a random process.
EpiAware is built using the DynamicPPL probabilistic programming domain-specific language, which is part of the Turing PPL. The structural concept behind EpiAware is that each module of an epidemiological model is a self-contained TuringModel; that is each module is an object that can be conditioned on observable data and sampled from. A complete EpiAware model is the composition of these objects using the @submodel macro.
+
To demonstrate EpiAware we largely recreate an epidemiological model presented in On the derivation of the renewal equation from an age-dependent branching process: an epidemic modelling perspective, Mishra et al (2020). Mishra et al consider test-confirmed cases of COVID-19 in South Korea between January to July 2020. The components of the epidemilogical model they consider are:
$$I_t = \mu_t + R_t \sum_{s\geq 1} I_{t-s} g_s.$$
Where \(g_t\) is a daily discretisation of the probability mass function of an estimated serial interval distribution:
$$G \sim \text{Gamma}(6.5,0.62).$$
And \(\mu_t\) is an external importation of infection process.
$$C_t \sim \text{NegBin}(\text{mean} = I_t,~ \text{overdispersion} = \phi).$$
In the examples below we are going to largely recreate the Mishra et al model, whilst emphasing that each component of the overall epidemiological model is, itself, a stand alone model that can be sampled from.
TaskLocalRNG()
+
+
+Time-varying reproduction number as a LatentModel type
EpiAware exposes a LatentModel type system; the purpose of which is to define stochastic processes which can be interpreted as generating time-varying parameters/quantities of interest.
In the Mishra et al model the log-time varying reproductive number is a priori assumed to evolve as an auto-regressive process, AR(2):
$$\begin{align}
+Z_t &= \log R_t, \\
+Z_t &= \rho_1 Z_{t-1} + \rho_2 Z_{t-2} + \epsilon_t, \\
+\epsilon_t &\sim \text{Normal}(0, \sigma).
+\end{align}$$
EpiAware gives a concrete subtype AR <: AbstractLatentModel which defines this behaviour of the latent model. The user can supply the priors for \(\rho_1,\rho_2\), wich we call damp_priors, as well as for \(\sigma\) (std_prior) and the initial values \(Z_1, Z_2\) (init_priors).
ar = AR(
+ damp_priors = [truncated(Normal(0.8, 0.05), 0, 1),
+ truncated(Normal(0.05, 0.05), 0, 1)],
+ std_prior = HalfNormal(1.0),
+ init_priors = [Normal(-1.0, 0.1), Normal(-1.0, 0.1)]
+)
+AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2)
+
+
+The priors here are based on Mishra et al, note that we have decreased the a priori belief in the correlation parameter \(\rho_1\).
Turing model interface
As mentioned above, we can use this instance of the AR latent model to construct a TuringModel which implements the probabilistic behaviour determined by ar.
We do this with the constructor function generate_latent which combines ar with a number of time steps to generate for (in this case we choose 30).
ar_mdl = generate_latent(ar, 30)
+Model{typeof(generate_latent), (:latent_model, :n), (), (), Tuple{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64}, Tuple{}, DefaultContext}(EpiAware.EpiAwareBase.generate_latent, (latent_model = AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2), n = 30), NamedTuple(), DefaultContext())
+
+
+We can sample from this model, which is useful for model diagnostic and prior predictive checking.
+
+
+And we can sample from this model with some parameters conditioned, for example with \(\sigma = 0\). In this case the AR process is an initial perturbation model with return to baseline.
cond_ar_mdl = ar_mdl | (σ_AR = 0.0,)
+Model{typeof(generate_latent), (:latent_model, :n), (), (), Tuple{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64}, Tuple{}, ConditionContext{@NamedTuple{σ_AR::Float64}, DefaultContext}}(EpiAware.EpiAwareBase.generate_latent, (latent_model = AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2), n = 30), NamedTuple(), ConditionContext((σ_AR = 0.0,), DynamicPPL.DefaultContext()))
+
+
+
+
+
+In this note, we are going to treat \(R_t\) as varying every two days. The reason for this is to 1) reduce the effective number of parameters, and 2) showcase the BroadcastLatentModel wrapper.
In EpiAware we set this behaviour by wrapping a LatentModel in a BroadcastLatentModel. This allows us to set the broadcasting period and type. In this case we broadcast each latent process value over \(2\) days in a RepeatBlock.
twod_ar = BroadcastLatentModel(ar, 2, RepeatBlock())
+BroadcastLatentModel{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64, RepeatBlock}(AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2), 2, RepeatBlock())
+
+
+
+
+
+The Renewal model as an EpiModel type
EpiAware has an EpiModel type system which we use to set the behaviour of the latent infection model. In this case we want to implement a renewal model.
To construct an EpiModel we need to supply some fixed data for the model contained in an EpiData object. The EpiData constructor performs double interval censoring to convert our continuous estimate of the generation interval into a discretized version \(g_t\). We also implement right truncation, the default is rounding the 99th percentile of the generation interval distribution, but this can be controlled using the keyword D_gen.
truth_GI = Gamma(6.5, 0.62)
+Distributions.Gamma{Float64}(α=6.5, θ=0.62)
+
+model_data = EpiData(gen_distribution = truth_GI)
+EpiData{Float64, typeof(exp)}([0.026663134095601056, 0.14059778064943784, 0.2502660305615846, 0.24789569560506844, 0.1731751163417783, 0.09635404000022223, 0.04573437575216367, 0.019313826994143808], 8, exp)
+
+
+
+
+
+The user also needs to specify a prior for the log incidence at time zero, \(\log I_0\). The initial history of latent infections \(I_{-1}, I_{-2},\dots\) is constructed as
$$I_t = e^{rt} I_0,\qquad t = 0, -1, -2,...$$
Where the exponential growth rate \(r\) is determined by the initial reproductive number \(R_1\) via the solution to the implicit equation,
$$R_1 = 1 \Big{/} \sum_{t\geq 1} e^{-rt} g_t$$
log_I0_prior = Normal(log(1.0), 1.0)
+Distributions.Normal{Float64}(μ=0.0, σ=1.0)
+
+epi = RenewalWithPopulation(model_data, log_I0_prior, 1e8)
+RenewalWithPopulation{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.026663134095601056, 0.14059778064943784, 0.2502660305615846, 0.24789569560506844, 0.1731751163417783, 0.09635404000022223, 0.04573437575216367, 0.019313826994143808], 8, exp), Distributions.Normal{Float64}(μ=0.0, σ=1.0), 1.0e8)
+
+
+NB: We don't implement a background infection rate in this model.
Turing model interface
As mentioned above, we can use this instance of the Renewal latent infection model to construct a TuringModel which implements the probabilistic behaviour determined by epi.
We do this with the constructor function generate_latent_infs which combines epi with a provided \(\log R_t\) time series.
Here we choose an example where \(R_t\) decreases from \(R_t = 3\) to \(R_t = 0.5\) over the course of 30 days.
R_t_fixed = [0.5 + 2.5 / (1 + exp(t - 15)) for t in 1:30]
+30-element Vector{Float64}:
+ 2.9999979211799306
+ 2.9999943491892553
+ 2.9999846395634946
+ 2.99995824644538
+ 2.9998865053282437
+ 2.9996915135600344
+ 2.999161624673834
+ ⋮
+ 0.500113494671756
+ 0.5000417535546202
+ 0.5000153604365055
+ 0.5000056508107448
+ 0.5000020788200692
+ 0.5000007647555673
+
+latent_inf_mdl = generate_latent_infs(epi, log.(R_t_fixed))
+Model{typeof(generate_latent_infs), (:epi_model, :_Rt), (), (), Tuple{RenewalWithPopulation{Normal{Float64}}, Vector{Float64}}, Tuple{}, DefaultContext}(EpiAware.EpiAwareBase.generate_latent_infs, (epi_model = RenewalWithPopulation{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.026663134095601056, 0.14059778064943784, 0.2502660305615846, 0.24789569560506844, 0.1731751163417783, 0.09635404000022223, 0.04573437575216367, 0.019313826994143808], 8, exp), Distributions.Normal{Float64}(μ=0.0, σ=1.0), 1.0e8), _Rt = [1.0986115957278464, 1.098610405062754, 1.0986071685094998, 1.0985983707197156, 1.0985744563952262, 1.098509454567543, 1.0983327911702674, 1.097852790994088, 1.09654964358037, 1.0930193012626002 … -0.6808598639891831, -0.6886022683659516, -0.6914718340845708, -0.6925303979295329, -0.6929202169746165, -0.6930636769372295, -0.6931164601588107, -0.6931358790023187, -0.6931430229284499, -0.6931456510499804]), NamedTuple(), DefaultContext())
+
+
+
+
+
+Negative Binomial Observations as an ObservationModel type
In Mishra et al latent infections were assumed to occur on their observation day with negative binomial errors, this motivates using the serial interval (the time between onset of symptoms of a primary and secondary case) rather than generation interval distribution (the time between infection time of a primary and secondary case).
Observation models are set in EpiAware as concrete subtypes of an ObservationModel. The Negative binomial error model without observation delays is set with a NegativeBinomialError struct. In Mishra et al the overdispersion parameter \(\phi\) sets the relationship between the mean and variance of the negative binomial errors,
$$\text{var} = \text{mean} + {\text{mean}^2 \over \phi}.$$
In EpiAware, we default to a prior on \(\sqrt{1/\phi}\) because this quantity has the dimensions of a standard deviation and, therefore, is easier to reason on a priori beliefs.
obs = NegativeBinomialError(cluster_factor_prior = HalfNormal(0.15))
+NegativeBinomialError{HalfNormal{Float64}}(HalfNormal{Float64}(μ=0.15))
+
+
+Turing model interface
We can construct a NegativeBinomialError model implementation as a TuringModel using generate_observations
Turing uses missing arguments to indicate variables that are to be sampled. We use this to observe a forward model that samples observations, conditional on an underlying expected observation time series.
First, we set an artificial expected cases curve.
expected_cases = [1000 * exp(-(t - 15)^2 / (2 * 4)) for t in 1:30]
+30-element Vector{Float64}:
+ 2.289734845645553e-8
+ 6.691586091292782e-7
+ 1.5229979744712628e-5
+ 0.0002699578503363014
+ 0.003726653172078671
+ 0.04006529739295107
+ 0.33546262790251186
+ ⋮
+ 0.003726653172078671
+ 0.0002699578503363014
+ 1.5229979744712628e-5
+ 6.691586091292782e-7
+ 2.289734845645553e-8
+ 6.101936677605324e-10
+
+obs_mdl = generate_observations(obs, missing, expected_cases)
+Model{typeof(generate_observations), (:obs_model, :y_t, :Y_t), (), (:y_t,), Tuple{NegativeBinomialError{HalfNormal{Float64}}, Missing, Vector{Float64}}, Tuple{}, DefaultContext}(EpiAware.EpiAwareBase.generate_observations, (obs_model = NegativeBinomialError{HalfNormal{Float64}}(HalfNormal{Float64}(μ=0.15)), y_t = missing, Y_t = [2.289734845645553e-8, 6.691586091292782e-7, 1.5229979744712628e-5, 0.0002699578503363014, 0.003726653172078671, 0.04006529739295107, 0.33546262790251186, 2.187491118182885, 11.108996538242305, 43.93693362340742 … 11.108996538242305, 2.187491118182885, 0.33546262790251186, 0.04006529739295107, 0.003726653172078671, 0.0002699578503363014, 1.5229979744712628e-5, 6.691586091292782e-7, 2.289734845645553e-8, 6.101936677605324e-10]), NamedTuple(), DefaultContext())
+
+
+
+
+
+A reverse observation model, which samples the underlying latent infections conditional on observations would require a prior on the latent infections. This is the purpose of composing multiple models; as we'll see below the latent infection and latent \(R_t\) models are informative priors on the latent infection time series underlying the observations.
Composing models into an EpiProblem
As mentioned above, each module of the overall epidemiological model we are interested in is a TuringModel in its own right. In this section, we compose the individual models into the full epidemiological model using the EpiProblem struct.
The constructor for an EpiProblem requires:
An epi_model.
A latent_model.
An observation_model.
A tspan.
The tspan set the range of the time index for the models.
epi_prob = EpiProblem(epi_model = epi,
+ latent_model = twod_ar,
+ observation_model = obs,
+ tspan = (45, 80))
+EpiProblem{RenewalWithPopulation{Normal{Float64}}, BroadcastLatentModel{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64, RepeatBlock}, NegativeBinomialError{HalfNormal{Float64}}}(RenewalWithPopulation{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.026663134095601056, 0.14059778064943784, 0.2502660305615846, 0.24789569560506844, 0.1731751163417783, 0.09635404000022223, 0.04573437575216367, 0.019313826994143808], 8, exp), Distributions.Normal{Float64}(μ=0.0, σ=1.0), 1.0e8), BroadcastLatentModel{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64, RepeatBlock}(AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2), 2, RepeatBlock()), NegativeBinomialError{HalfNormal{Float64}}(HalfNormal{Float64}(μ=0.15)), (45, 80))
+
+
+
We make inferences on the unobserved quantities, such as \(R_t\) by sampling from the model conditioned on the observed data. We generate the posterior samples using the No U-Turns (NUTS) sampler.
To make NUTS more robust we provide manypathfinder, which is built on pathfinder variational inference from Pathfinder.jl. manypathfinder runs nruns pathfinder processes on the inference problem and returns the pathfinder run with maximum estimated ELBO.
The composition of doing variational inference as a pre-sampler step which gets passed to NUTS initialisation is defined using the EpiMethod struct, where a sequence of pre-sampler steps can be be defined.
EpiMethod also allows the specification of NUTS parameters, such as type of automatic differentiation, type of parallelism and number of parallel chains to sample.
num_threads = min(10, Threads.nthreads())
+1
+
+inference_method = EpiMethod(
+ pre_sampler_steps = [ManyPathfinder(nruns = 4, maxiters = 100)],
+ sampler = NUTSampler(adtype = AutoReverseDiff(true),
+ ndraws = 2000,
+ nchains = num_threads,
+ mcmc_parallel = MCMCThreads())
+)
+EpiMethod{ManyPathfinder, NUTSampler{AutoReverseDiff, MCMCThreads, UnionAll}}(ManyPathfinder[ManyPathfinder(10, 4, 100, 100)], NUTSampler{AutoReverseDiff, MCMCThreads, UnionAll}(0.8, AutoReverseDiff(true), MCMCThreads(), 1, 10, 1000.0, 0.0, 2000, AdvancedHMC.DiagEuclideanMetric))
+
+
+
In the background of this note (see hidden top cell and short R script in this directory), we load daily reported cases from South Korea from Jan-July 2020 which were gathered using covidregionaldata from ECDC data archives.
We supply the data as a NamedTuple with the y_t field containing the observed data, shortened to fit the chosen tspan of epi_prob.
south_korea_data = (y_t = daily_cases[epi_prob.tspan[1]:epi_prob.tspan[2]],
+ dates = dates[epi_prob.tspan[1]:epi_prob.tspan[2]])
+(y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], dates = [Date("2020-02-13"), Date("2020-02-14"), Date("2020-02-15"), Date("2020-02-16"), Date("2020-02-17"), Date("2020-02-18"), Date("2020-02-19"), Date("2020-02-20"), Date("2020-02-21"), Date("2020-02-22") … Date("2020-03-10"), Date("2020-03-11"), Date("2020-03-12"), Date("2020-03-13"), Date("2020-03-14"), Date("2020-03-15"), Date("2020-03-16"), Date("2020-03-17"), Date("2020-03-18"), Date("2020-03-19")])
+
+
+Sampling with apply_method
The apply_method function combines the elements above:
And returns a collection of results:
inference_results = apply_method(epi_prob,
+ inference_method,
+ south_korea_data
+)
+EpiAwareObservables(Model{typeof(generate_epiaware), (:y_t, :time_steps, :epi_model), (:latent_model, :observation_model), (), Tuple{Vector{Int64}, Int64, RenewalWithPopulation{Normal{Float64}}}, Tuple{BroadcastLatentModel{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64, RepeatBlock}, NegativeBinomialError{HalfNormal{Float64}}}, DefaultContext}(EpiAware.EpiAwareBase.generate_epiaware, (y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], time_steps = 36, epi_model = RenewalWithPopulation{Normal{Float64}}(EpiData{Float64, typeof(exp)}([0.026663134095601056, 0.14059778064943784, 0.2502660305615846, 0.24789569560506844, 0.1731751163417783, 0.09635404000022223, 0.04573437575216367, 0.019313826994143808], 8, exp), Distributions.Normal{Float64}(μ=0.0, σ=1.0), 1.0e8)), (latent_model = BroadcastLatentModel{AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}, Int64, RepeatBlock}(AR{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}}}, HalfNormal{Float64}, DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}, Int64}(Distributions.Product{Distributions.Continuous, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}, Vector{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64, Float64, Float64}}}(v=Truncated{Normal{Float64}, Continuous, Float64, Float64, Float64}[Truncated(Distributions.Normal{Float64}(μ=0.8, σ=0.05); lower=0.0, upper=1.0), Truncated(Distributions.Normal{Float64}(μ=0.05, σ=0.05); lower=0.0, upper=1.0)]), HalfNormal{Float64}(μ=1.0), DistributionsAD.TuringScalMvNormal{Vector{Float64}, Float64}(m=[-1.0, -1.0], σ=0.1), 2), 2, RepeatBlock()), observation_model = NegativeBinomialError{HalfNormal{Float64}}(HalfNormal{Float64}(μ=0.15))), DefaultContext()), nothing, MCMC chain (2000×35×1 Array{Float64, 3}), @NamedTuple{generated_y_t::Vector{Int64}, I_t::Vector{Float64}, Z_t::Vector{Float64}, process_aux::@NamedTuple{σ_AR::Float64, ar_init::Vector{Float64}, damp_AR::Vector{Float64}, sq_cluster_factor::Float64}}[(generated_y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], I_t = [0.6302189734369976, 0.4774713968245689, 0.3905291507562299, 0.3022359509208115, 0.9922833853995451, 0.813228204790894, 21.106447028610603, 46.10331357944557, 53.33687089383859, 152.24222960039674 … 197.5819007248144, 173.94553409919794, 136.98722629438038, 111.65633481799061, 125.31743099837803, 99.6995783777211, 78.81882944289, 68.09834544453484, 113.75830423055085, 99.6034860166183], Z_t = [-1.0832047855152138, -1.0832047855152138, -1.0146262702530706, -1.0146262702530706, 0.418711014559223, 0.418711014559223, 3.6849734508145495, 3.6849734508145495, 2.414227285537577, 2.414227285537577 … -0.8921461057175676, -0.8921461057175676, -0.9481069026492762, -0.9481069026492762, -0.6086302626701813, -0.6086302626701813, -0.6492378533259052, -0.6492378533259052, 0.0048214101954354716, 0.0048214101954354716], process_aux = (σ_AR = 1.5841436659699561, ar_init = [-1.0832047855152138, -1.0146262702530706], damp_AR = [0.7640084446616473, 0.06809055715841618], sq_cluster_factor = 0.079180832239541)); (generated_y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], I_t = [0.9054036333499285, 0.7298094574188798, 0.5567385913041931, 0.4503846504021242, 2.230213896626074, 1.9065600991663039, 17.881558039189844, 32.4744970077826, 51.89540975761122, 128.44184641355815 … 178.21457047327442, 153.92366069174773, 150.53268696146284, 118.86938255651403, 120.72629891594556, 96.25096518758762, 111.96747397950845, 95.52634405260403, 165.56194047256344, 157.2302637345635], Z_t = [-0.8898357447217325, -0.8898357447217325, -0.9467061642702784, -0.9467061642702784, 0.8678971922522976, 0.8678971922522976, 2.9833461455222965, 2.9833461455222965, 2.4248084414242914, 2.4248084414242914 … -1.108402692341453, -1.108402692341453, -0.9535543249686977, -0.9535543249686977, -0.6756561197376965, -0.6756561197376965, -0.3264624005485029, -0.3264624005485029, 0.33979129660429197, 0.33979129660429197], process_aux = (σ_AR = 2.006219571534054, ar_init = [-0.8898357447217325, -0.9467061642702784], damp_AR = [0.7477859993710714, 0.07099853573426045], sq_cluster_factor = 0.10824590568528863)); … ; (generated_y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], I_t = [1.2739313505255094, 0.9284479089848867, 0.8675168996177901, 0.6575416056195217, 1.7933463862437304, 1.4831973604222313, 18.312506724736963, 28.3160889159895, 82.39544494211103, 215.99113564106977 … 153.63601628462584, 135.00614788819811, 154.28236495292967, 124.1945925971968, 128.01665725054332, 106.31031490775568, 78.72349016967804, 70.33397592850444, 155.20469776895922, 138.04015530663935], Z_t = [-1.1712194264050033, -1.1712194264050033, -0.9365858749109013, -0.9365858749109013, 0.31385097788297855, 0.31385097788297855, 2.8525497842836396, 2.8525497842836396, 2.927661598221874, 2.927661598221874 … -1.0354808546735808, -1.0354808546735808, -0.7278311981768211, -0.7278311981768211, -0.47094674087799177, -0.47094674087799177, -0.6315540092573445, -0.6315540092573445, 0.28437585348657035, 0.28437585348657035], process_aux = (σ_AR = 1.638797767408246, ar_init = [-1.1712194264050033, -0.9365858749109013], damp_AR = [0.7615374108450998, 0.061451939399942826], sq_cluster_factor = 0.1010913711477206)); (generated_y_t = [0, 0, 0, 1, 1, 1, 15, 34, 75, 190 … 131, 242, 114, 110, 107, 76, 74, 84, 93, 152], I_t = [0.36557026373911194, 0.28537577652833657, 0.24465754806075726, 0.19332007790512695, 0.5784430227066278, 0.4831367513694791, 23.986389990707373, 74.85371802000027, 40.67406402063849, 127.27792831196787 … 200.46755591770244, 177.25464086769261, 169.4193439834972, 139.64155100487525, 117.50767176549644, 98.4990070786369, 75.57177784500391, 64.52735423543758, 200.2979491781433, 173.46136967899554], Z_t = [-0.9990885088367124, -0.9990885088367124, -0.9097637407764456, -0.9097637407764456, 0.41166688293517095, 0.41166688293517095, 4.320290495933675, 4.320290495933675, 1.963418480789977, 1.963418480789977 … -0.8273355767471939, -0.8273355767471939, -0.7012007179628654, -0.7012007179628654, -0.6810786163999226, -0.6810786163999226, -0.7819315466370461, -0.7819315466370461, 0.5240376203663384, 0.5240376203663384], process_aux = (σ_AR = 2.0368151186594696, ar_init = [-0.9990885088367124, -0.9097637407764456], damp_AR = [0.7551398290101086, 0.07370704512817201], sq_cluster_factor = 0.12321894564496723));;])
+
+
+Results and Predictive plotting
We can spaghetti plot generated case data from the version of the model which hasn't conditioned on case data using posterior parameters inferred from the version conditioned on observed data. This is known as posterior predictive checking, and is a useful diagnostic tool for Bayesian inference (see here).
Because we are using synthetic data we can also plot the model predictions for the unobserved infections and check that (at least in this example) we were able to capture some unobserved/latent variables in the process accurate.
We find that the EpiAware model recovers the main finding in Mishra et al; that the \(R_t\) in South Korea peaked at a very high value (\(R_t \sim 10\) at peak) before rapidly dropping below 1 in early March 2020.
Note that, in reality, the peak \(R_t\) found here and in _Mishra et al) is unrealistically high, this might be due to a combination of:
A mis-estimated generation interval/serial interval distribution.
An ascertainment rate that was, in reality, changing over time.
In a future note, we'll demonstrate having a time-varying ascertainment rate.
+
+
+Parameter inference
We can interrogate the sampled chains directly from the samples field of the inference_results object.
+
+