Mapping Incidence and Prevalence Peak Data for SIR Forecasting Applications
Alexander C. Murph G. Casey Gibson Lauren J. Beesley Nishant Panda Lauren A. Castro Sara Y. Del Valle Dave Osthus





Incorporating Incidence Peak and Time Data for SIR Forecasting Applications












A]Alexander C. Murph[label=e1]murph@lanl.gov0000-0001-7170-867X ,
A]G. Casey Gibson[label=e2]gcgibson@lanl.gov0000-0002-0370-9846,
A]Lauren J. Beesley[label=e3]lvandervort@lanl.gov0000-0002-3788-5944,
B,C]Nishant Panda[label=e4]nishpan@lanl.gov0000-0001-9754-2794 ,
B]Lauren A. Castro[label=e5]lcastro@lanl.gov0000-0002-9778-570X ,
B]Sara Del Valle[label=e6]sdelvall@lanl.gov0000-0002-0159-1952 ,
A]Dave Osthus[label=e7]dosthus@lanl.gov0000-0002-4681-091X 




[A]Statistical Sciences (CCS-6), Computer, Computational, and Statistical Sciences
Division, Los Alamos National Laboratory, Los Alamos, NM, 87545 [presep=, ]e1,e2,e3,e7
[B]Data Analytics & Forecasting, A-1: Information Systems & Modeling, Los Alamos National Laboratory, Los Alamos, NM, 87545 [presep=, ]e5,e6
[C]Information Sciences (CCS-3), Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545 [presep=, ]e4





Infectious disease modeling and forecasting have played a key role in helping assess and respond to epidemics and pandemics. Recent work has leveraged data on disease peak infection and peak hospital incidence to fit compartmental models for the purpose of forecasting and describing the dynamics of a disease outbreak.  Incorporating these data can greatly stabilize a compartmental model fit on early observations, where slight perturbations in the data may lead to model fits that project wildly unrealistic peak infection.  We introduce a new method for incorporating historic data on the value and time of peak incidence of hospitalization into the fit for a Susceptible-Infectious-Recovered (SIR) model by formulating the relationship between an SIR model's starting parameters and peak incidence as a system of two equations that can be solved computationally.  This approach is assessed for practicality in terms of accuracy and speed of computation via simulation.  To exhibit the modeling potential, we update the Dirichlet-Beta State Space modeling framework to use hospital incidence data, as this framework was previously formulated to incorporate only data on total infections.



Compartmental Models
Disease Forecasting
Hospital Incidence
Prevalence











§ INTRODUCTION

Compartmental models have seen broad usage at the onset of several disease outbreaks in the last century as a means to project expected numbers of infected and to drive healthcare response.  Broadly speaking, a compartmental model describes the dynamics of a disease spread by breaking a population into set categories and modeling the process by which individuals transition through these categories.  Perhaps the most basic is the Kermack–McKendrick model, often called the SIR or Susceptible-Infectious-Recovered model, which models the movement of subjects from being Susceptible, to Infected (and contagious), and then to Removed <cit.>.  Since the development of the SIR model, further research has extended the idea to include additional compartments, such as the Exposed category in the SEIR model – where a subject has been exposed to the disease but is not yet infectious – and the ability to move back to the Susceptible category in the SEIS model (see <cit.> for an overview of each).  Compartmental models have been applied to projection tasks such as the 2014-15 Ebola epidemic <cit.>, the 2009 A/H1N1 influenza pandemic <cit.>, several HIV outbreaks <cit.>, the recent COVID-19 pandemic <cit.>, and many other epidemiological projection tasks in the last century (see, for instance, <cit.> and references therein).

The elegance of compartmental models is in their succinct ability to describe the state of a population in terms of how subjects transfer to and from the different categories, and thus fitting these models requires estimation of interpretable quantities such as rates of infection and recovery.  For instance, denote the proportion of the population in the susceptible, infectious, and removed compartments by S_t, I_t, and R_t, respectively, such that S_t + I_t + R_t = 1 for all t.  Then the SIR model is determined by the equations


    dS_t/dt   = - β S_t I_t , 
    dI_t/dt   = β S_t I_t - γ I_t , 
    dR_t/dt   = γ I_t  ,


where β > 0 is the disease transmission rate and γ >0 is the rate of recovery.  If one knew these two rates, and the initial number of individuals in each category – S_0, I_0 and R_0 – the numbers S_t, I_t, and  R_t could be numerically simulated for any time-point t.  While modern computational tools do make simulating these quantities feasible, the need to simulate the entire system numerically to get the exact counts in each compartment is challenging. That is, fitting these models to data and can be computationally expensive for more intricate compartmental structures.  Researchers have determined analytic solutions for the entire SIR model (explicit forms/approximations to the number of susceptible, infected, and removed at a specific time) that do not require numeric simulation <cit.>.  However, these solutions involve reparameterizing the time axis, and explicit calculations back onto the original time axis require numeric integration or approximation methods.  Similar approaches are also used to obtain exact and approximate solutions to the more complicated compartmental models, such as the SEIR model <cit.> and the SIRS model <cit.>.

Analytic maps from the initial starting parameters of an SIR curve to quantities of interest (QoIs) – such as the value of the peak of the infected curve and the limiting number of susceptible individuals – were developed as early as the mid-1900s <cit.>.  The focus on these quantities has generally been motivated by their usefulness to making public health decisions.  For instance, in papers such as <cit.> and <cit.>, the maximum value of the infected compartment (not the time of maximum infection) was studied for the purpose of informing public health officials of the maximum number of infections after the initial estimation of the disease transmission and recovery rates, since knowing this maximum quantity informs how many hospital beds might be needed.  The authors in <cit.> point out that the time of peak infection is also informative for healthcare officials, and they develop an analytic form for the peak time of infection using an approximation of the SIR curve.  An exact form for the peak infection time is available in <cit.> and <cit.>, albeit in terms of an integral without a closed form.  Several modern papers study fast approximations to the integral form for peak infection time <cit.>, using the previously-mentioned analytic solutions to the SIR model. These different approximations for peak infection time are evaluated and compared in <cit.>.

Some recent papers study SIR curve QoIs for inferential tasks related to modeling an epidemic <cit.>.  In <cit.>, the Gumbel distribution is suggested as a good approximation of the infection curve in an SIR model, and maps between the SIR curve parameters and peak value/time are used to develop Method of Moments estimators of the Gumbel parameters.  The Gumbel distribution is then used to approximate the exact solution of the SIR model.  In <cit.>, the relationship between peak infection value and the SIR parameters is used to incorporate historical data on epidemic peaks into the inferential task of fitting an SIR curve during the early stages of an epidemic. The authors point out that an SIR curve is sensitive to small perturbations in the transmission and recovery rates, and that incorporating these data discourages models that well-represent early data but drastically over-predict the peak infection value.  In <cit.>, both historic and surveyed QoIs are used to constrain an SEIRH model (where the added “H" category refers to people hospitalized at a given time) in much the same way as <cit.>.



Both <cit.> and <cit.> use information on peak infection value and time to inform the modeling of a pandemic and epidemic, respectively. Since most modern monitoring systems approximate the daily number of new cases of a disease, called incidence, this is not the most practically useful development.  The number of infected individuals at any given time, called prevalence, is typically an unobserved quantity <cit.>.  For examples of all quantities in an SIR curve on a fixed time axis, see Figure <ref>.  Note that peak prevalence and peak incidence need not occur at the same time, and that it is mathematically possible for the incidence to be greater than prevalence (as visualized on the right side of Figure <ref>).  This being said, in most disease outbreaks of note, prevalence is typically greater than incidence.

The model in <cit.> incorporates incidence QoI data by extending the parameter space to include peak incidence value (PIV) and peak incidence time (PIT) and by defining a prior on these values using either historic or survey data.  While this model achieves a similar goal to the ultimate goal of this paper, its main distinction is that that it requires an importance sampling scheme to fit, which we outline here for completeness: priors on the rate parameters and initial values for the compartmental model are used to sample proposal values, the entire system is numerically simulated to determine the (PIV, PIT), then the likelihood of this system is determined via the prior probability of (PIV, PIT).  Note that while this need not be true for the SEIRH model considered in <cit.> (since it is not an SIR model), defining both a prior on the rate parameters and on PIV and PIT for an SIR model is redundant given the initial values S_0, I_0, and R_0, since a given set of starting parameters should immediately determine these QoIs.  

Incorporating historical peak values and times of an outbreak into a Bayesian compartmental forecasting model via prior specifications is highly useful to constrain forecasts[For a review of the differences between disease forecasts and disease projections, see <cit.>.], especially early in the outbreak (see <cit.>, <cit.>, and <cit.>).  Relating these historical QoIs requires determining a relationship between historical peak values and times as well as the parameters and initial conditions of the compartmental model (e.g., the SIR model).  The challenge lies in this relationship. Analytically, this relationship is known for observed peak prevalence value (PPV) and peak prevalence time (PPT). However, almost all infectious disease data is of incidence. Despite being principally unsound, the data application in <cit.> treated incidence data as if it were prevalence in order to use the known analytic relationships.  In this paper, we make three contributions.  First, we develop the methods to map peak values and times of incidence to the parameters of the SIR model, given the initial conditions.  Second, we demonstrate how to incorporate these new mappings into the model of <cit.> and provide the code used to do so.   Third, we compare the forecasts of the modeling framework of <cit.> using both incidence data as incidence, and using the misspecified prevalence data.  In addition to comparing these forecasts, we show that if SIR parameter inference is desired, then mistaking incidence data for prevalence data will result in biased estimation.

This paper is outlined as follows.  In Section <ref>, we review existing methods for mapping SIR parameters to peak/time of prevalence.  In Section <ref>, we develop identical maps for incidence, then provide methods for inverting these maps for both incidence and prevalence.  In Section <ref>, we develop and improve upon the modeling framework in <cit.> to incorporate historic incidence peak/time data when fitting an SIR curve.  We apply this model to influenza data in Section <ref>.



§ ANALYTIC SOLUTIONS TO THE SIR EQUATIONS FOR INCORPORATING PREVALENCE DATA
 
The system (<ref>)-(<ref>) can be solved analytically for all time by reparameterizing the time axis to be instead in terms of the number of individuals removed from the system beyond the initial amount in the Removed category R_0,

    
    S_t    = S_0 e^-βτ_t, 
    
            I_t    = S_0 + I_0 - S_0 e^-βτ_t - γτ_t , 
    
            R_t    = R_0 + γτ_t,


where τ_t is the inverse of the map

    t_τ := ∫_0^τ[ dτ'/S_0 + I_0 - S_0 e^βτ' - γτ']^+,

where [x]^+ := max{x,0}. The above form for the SIR dynamics is particularly useful since it allows one to calculate the number of individuals in each category without numerically simulating the entire system.  The major drawback of this form is that mapping values back to an interpretable time requires one to approximate the integral in (<ref>), since this integral is nonelementary.  There exist other methods that derive an analytic solution to the SIR system <cit.>.  We use this formulation since it is well-known (having been originally developed in <cit.>), and because it simplifies the calculation of SIR curve QoIs, such as PPV and PPT.  For instance, the maximum of (<ref>) with respect to τ occurs at

    τ^* = 1/βlog( β S_0/γ),

which can be derived via straightforward calculus.  PPT can then be approximated from (<ref>), and PPV calculated from (<ref>).

The process above of calculating PPV and PPT, denoted (I_t_τ^*, t_τ^*), from the initial starting parameters can be reversed for fixed S_0, I_0, R_0.  That is, given (I_t_τ^*, t_τ^*), we propose a method to calculate (β, γ).  Define ρ as 

    ρ = β/γ.

The quantity S_0 ρ is more widely known as the basic reproduction number, which measures the average number of additional infections generated by a single new infection <cit.>.  Plugging (<ref>) into (<ref>) gives

    I_t = S_0 + I_0 - 1/ρ - log (ρ S_0)/ρ.

From here, setting I_t:= I_t_τ^* and solving for ρ is a (piecewise) convex optimization problem, which is typically fast and accurate computationally.  Note that a change of variables for (<ref>) gives

    β t_τ = ∫_0^βτ[ dτ̂/S_0 + I_0 - S_0 e^-τ̂ - τ̂ / ρ]^+ ,

and that βτ^* = log(ρ S_0).  Thus, once the value of ρ is approximated, (<ref>) can be evaluated to get β t_τ^*.  Since t_τ^* is assumed known a priori, the values of β and γ can be obtained via arithmetic.  

We introduce the above derivations not only because they will be used in the next Section, but also to motivate the type of calculation we aim to develop in this paper.  While the above approach still requires computational methods to map from PPV and PPT to the SIR model parameters, this approach is much more direct than the brute force method of simulating several (β, γ) combinations until a prevalence curve with a peak sufficiently near (I_t_τ^*, t_τ^*) is discovered <cit.>. 



§ MAPPING INCIDENCE DATA TO SIR PARAMETERS


As mentioned above, a primary focus of this work will be to develop maps between PIV and PIT and the initial SIR curve parameters.  In this direction, one can derive an equivalent formulation of the SIR dynamics in (<ref>) - (<ref>) by replacing (<ref>) with

    I_t    = I_t-1 + β S_t-1 I_t-1 - γ I_t-1.

Putting (<ref>) into words, prevalence is equal to the prevalence at the last time step, plus those in the infected category that infect those in the susceptible category, minus the number in the infected category that are removed from the system.  The form for prevalence in (<ref>) is particularly useful for this application since it contains a term that explicitly models incidence at time t: β S_t-1 I_t-1.



 §.§ Mapping SIR Parameters to Peak Incidence Value and Time
 
Reparameterizing the time axis for the term β S_t I_t using (<ref>) and (<ref>) gives the following form

    β (S_0 e^-βτ_t) (S_0 + I_0 - S_0 e^-βτ_t - γτ_t ).

The value for τ that maximizes (<ref>) also satisfies the equation,

    -(S_0 + I_0) + γτ + 2S_0 e^-βτ - 1/ρ = 0.

Solving (<ref>) is also a convex optimization problem on a single parameter, and is thus feasible and accurate to do numerically.  A solution to (<ref>) gives the τ for the timepoint directly before the time of max incidence, τ_t^*-1, where t^* denotes the timepoint of max incidence.  Of course, this is all that is required to calculate the PIV, β S_t-1 I_t-1, by a direct application of (<ref>).  The calculation of PIT is similarly straightforward.  Using (<ref>), the prevalence at τ_t^*-1 can be directly calculated.  The prevalence value for τ_t^*-1 and (<ref>) then gives τ_t^*, and (<ref>) can then be used to calculate t^*,

    t^* = ∫_0^τ_t^*[ dτ'/S_0 + I_0 - S_0 e^βτ' - γτ']^+.




 §.§ Mapping Incidence Peak Value and Time to SIR Parameters
 

The proposed method to map PIV and PIT back to the parameters of an SIR curve solves the following system of equations implied by (<ref>) and (<ref>):

    β (S_0 e^-βτ^⋆) (S_0 + I_0 - S_0 e^-βτ^⋆ - γτ^⋆ )    = PIV
    
        -(S_0 + I_0) + γτ^⋆ + 2S_0 e^-βτ - 1/ρ   = 0 ,

where τ^⋆ solves the equation

    ∫_0^τ^⋆[ dτ'/S_0 + I_0 - S_0 e^βτ' - γτ']^+ = PIT.


While solving a system of two equations with two unknowns (for β and γ) is generally feasible computationally, the major bottleneck for this problem is the need to invert (<ref>).  A brute-force computational approach to solving this system of equations would require one to both invert and solve (<ref>) for every value of (β, γ) investigated.  While experiments in this direction have proven to be surprisingly fast and accurate, the confounding computational approximations encourage a more analytic solution, or alternative computational strategies.  We consider three possible alternatives to computationally estimating the integral (referred to the “Compute Integral" method in the subsequent).



  §.§.§ Taylor Approximation
 
Given a candidate (β, γ) in any numerical solver, we approximate τ^⋆ by taking the second degree Taylor expansion of the integral in (<ref>) and solving for τ^⋆ algebraically <cit.>.  This leads to the following closed-form approximation,

    τ^⋆ = β^2/S_0[ ( ρ S_0 - 1 ) + κtanh( γκ (PIT)/2 - ϕ)] + R_0,

where

    κ   = √(( S_0 ρ - 1 )^2 + 2 S_0 I_0 ρ^2 ), 
    ϕ   = 1/κarctanh[ S_0 ρ - 1 ].

Using this approximation for τ^⋆, we numerically solve the system of equations expressed by (<ref>) and (<ref>).



  §.§.§ Single ODE Approximation
 
Instead of using a Taylor approximation to estimate τ^⋆ it is possible to do so by numerically solving an Ordinary Differential Equation (ODE).  Using the definition of τ, we combine equations (<ref>) and (<ref>) to get

    S_t = S_0 exp( -ρ ( R_t - R_0 ) ).

From here, we use (<ref>) and the assumption that I_t = 1 - S_t - R_t to get the ODE,

    dR_t/dt = γ[ 1 - S_0 e^-ρ(R_t - R_0) - R_t ].

The authors in <cit.> point out that while (<ref>) is a transcendental equation, it can be solved numerically, and has a single, unique solution by the general existence and uniqueness theorem for the solutions of ODEs.  Using this approximation for τ^⋆, we numerically solve the system of equations expressed by (<ref>) and (<ref>).



  §.§.§ Full ODE Approximation
 
As a final computational method for mapping PIV and PIT to β and γ, we numerically solve the system of ODEs described in (<ref>) - (<ref>).  In an optimization algorithm, this would require the ODE to be solved for every possible (β, γ) pair queried.  Much like the brute force computational method, we expect this method to be accurate, but to come at a high computational cost.



  §.§.§ Comparision via Simulation
  We compare all the above approaches via simulation, repeating the following steps 1000 times:

    
  * Sample a (PPV, PPT) from the bivariate normal 𝒩(μ, Σ) where μ = (0.0144, 17.9),
    
    Σ = [ 0.000036   0.0187;  -0.0187    16.09 ],

    truncated so that PPV ∈ (θ_I_0 , 1) and PPT ∈ (1,35).  This corresponds to the set of feasible values and sampling distribution described in <cit.>;
    
  * Use the method from Section <ref> to map this back to (β, γ), then numerically simulate the system in (<ref>)-(<ref>) to get the “true" values for PIT and PIV;
    
  * With PIT and PIV, use each method described above to find approximate values (β̂, γ̂);
    
  * For all approximation methods, compare the estimated PIT and PIV (gotten by numerically simulating the system from the appropriate (β̂, γ̂)) against the true PIT and PIV.

We outline the results from the above simulation in Table <ref>.  While the fastest method is the Taylor approximation, this method is also the least accurate. This is as expected, since this approximation is only accurate for sufficiently small values of ρ (R_t - R_0) <cit.>.  The Single ODE approximation method is slightly more accurate than the Taylor approximation, but at a higher computational cost.  The Compute Integral and the Full ODE approximation are comparable, with the Compute Integral method being more accurate and the Full ODE method being faster.

While the Compute Integral method is the slowest of these approaches, it is still somewhat fast (taking around half a second on average), and it is the most accurate overall (since the Taylor approximation and Single ODE methods are greatly off for PIV).  Since the data sizes for PIT and PIV data are not exorbitantly large, this computational burden would be acceptable in converting a data set of (PIV, PIT) values to (β̂, γ̂) values.  For the applications in this paper, we will use the Compute Integral method. Future work might examine alternative approximation methods, especially for approximating the inverse of the integral in (<ref>).  





§ A BAYESIAN STATE-SPACE SIR MODEL
 
In this section, we introduce the Dirichlet-Beta state-space model (DBSSM) from <cit.> and update it to incorporate historic PIV and PIT data.  The original formulation of the DBSSM was to answer an observed issue associated with using the SIR model for early-pandemic forecasting tasks.  Namely, that two SIR curves that reasonably fit early count data may lead to drastically different PPV predictions.  In a simulation example, <cit.> show two such SIR curves that have peaks that differ by 30% of the entire population, even though they have a nearly-identical fit to the early-pandemic data observations (see Figure 3 in the cited paper).  To address this stability issue, the DBSSM incorporates historic PPV and PPT data into the prior specifications to discourage SIR curve fits with peak values that are greatly above reasonable expectations.  As mentioned previously, this incorporation of prevalence data is not the most practical approach, since incidence data is generally the observed quantity.  A further shortcoming in the original DBSSM formulation is that it learns a map between PPT and the SIR curve parameters, rather than using an analytic map.  After introducing this model, we will identify ways that the methodology developed in this paper will improve these issues for the DBSSM.

Let y_t be the observed proportion in a population that tested positive for some disease at some timepoint t, and let θ_t = (S_t, I_t, R_t)'. Then the DBSSM is defined as,

    y_t | θ_t, ϕ∼Beta(λIn_t, λ (1 - In_t) )


    θ_t | θ_t-1, ϕ∼Dirichlet(ι f(θ_t-1, β, γ) ),

where ϕ = {S_0, I_0, R_0, β, γ, λ, ι}, λ, ι are variance control parameters, and f is a map that propagates the SIR system determined by (θ_t-1, β, γ) forward one step according to (<ref>)-(<ref>).  Note that, by this set up, the set of parameters θ_0:t' = {θ_0, θ_1, …, θ_t'} is a first-order Markov chain, and that for all s≠ t, the data observations y_s and y_t are independent given θ_t.  The variable In_t denotes the incidence of the system at time t; in the original formulation of this model, the prevalence – I_t – was used here.  The incidence at time t is directly calculated using θ_t and (<ref>).

The conditional expectations of the model described by (<ref>) and (<ref>) are unbiased,

    E(y_t | θ_t, ϕ) = In_t


    E(θ_t | θ_t-1, ϕ) = f(θ_t-1, β, γ)

while their respective variances reduce to zero as λ, ι→∞.  Of course, the conditional mean in (<ref>) is dependant upon the accuracy of f in propogating the latent space θ_t-1 forward one time step.  The authors in <cit.> used a fourth-order Runge-Kutta approximation and observed reasonable accuracy. 

We review the full Bayesian framework of the DBSSM and provide the prior specification in Appendix <ref>.  The main innovation of this model is that the parameter space is expanded to include the latent variable z = (PPT, PPV), and this latent variable is given a prior that incorporates historic PPT and PPV data.  In the following, we review this prior, π(z|θ_0), and the mechanism by which this prior informs the SIR parameters β and γ.  These priors are then each updated according to the theory developed in this paper.



 §.§ Specification of TEXT
 The prior on π(z | θ_0) in the DBSSM is a minimal-assumption distribution on historic data on PPT and PPV.  Note that this is the mechanism by which the authors in <cit.> directly address the aforementioned stability issues with fitting an SIR curve with early pandemic data.  They do so by fitting a normal distribution to historic influenza QoI data, truncated to enforce that an epidemic will occur (the lower bound on PPV was set to I_0) and so that the peak happens within the influenza forecasting season (PPT was required to be between the 1^st and 35^th weeks).  While somewhat loose, this prior gives very small (or zero) probability to values of (PPV, PPT) that are drastically outside of historically observed pandemics.  

In the formulation of the DBSSM developed in this paper, the same prior used on PPT and PPV is now used on PIT and PIV.  To connect this constraint into the model, we must next define how the assumption on this latent space affects the learning of the SIR parameters β and γ.



 §.§ Specification of TEXT
 







With the addition of the latent variable z, the prior needed for the SIR parameters is π(β, γ | z, θ_0).  In the original formulation of the DBSSM, this prior is reparameterized according to (ρ, γ), then factorized.  Thus, priors are instead given to π(ρ | z, θ_0) and π(γ | ρ, z, θ_0).  This additional formulation is done to utilize the following analytic relationship from <cit.>,

    PPV = g_1(S_0, I_0, ρ) = I_0 + S_0 - S_0 ρ[ log(S_0) + 1 - log(S_0 ρ) ].

By inverting this relationship, samples from the latent quantity z immediately determine the corresponding value of ρ.  The appropriate prior on this quantity would then be

    π(ρ | z, θ_0)  ∝ δ (S_0 ρ - g_1^-1(PPV, S_0, I_0)),

where δ is the Dirac delta function.  For the prior on γ, the map between γ and PPT is estimated using a simulated data set of 5250 SIR curves.  This map,

    PPT = g_2(S_0, I_0, ρ, γ),

was then used in lieu of any analytic form, and the prior on γ was set to

    π(ρ | z, θ_0)  ∝ δ (S_0 ρ - g_2^-1(PPV, S_0, I_0, ρ)),


We reiterate that there are two major shortcomings to the above prior specifications.  First, the above priors assume that there is access to historical PPV and PPT data, which is typically not the case, as public health data are generally on incidence, not prevalence.  In the original paper, incidence data were used instead of prevalence data without explicit justification.  Second, a map between PPT and the SIR parameters is estimated even though an analytic map between these quantities exists (see (<ref>) and (<ref>)), unnecessarily introducing a source of uncertainty.

The methods developed in this paper correct the limitations found in the original formulation of the DBSSM, since they provide maps to replace g_1, g_2 above with

    (PIV, PIT) = h(S_0, I_0, β, γ),

where h^-1 denotes the algorithm discussed in Section <ref>.  Thus, the joint prior used for β, γ in this updated version of the DBSSM is

    π(β, γ | z, θ_0)  ∝ δ( || (β, γ) - h^-1(PPV, S_0, I_0)||_1 ),

where ||· ||_1 is the 1-norm.

The naive treatment of incidence data as prevalence data (as was done in <cit.>), need not necessarily lead to a loss of forecasting accuracy in the final model.  An incidence curve can well approximate, or even be equivalent to, a prevalence curve.  For instance, consider the SIR model constrained so that γ = 1, which corresponds to the Reed-Frost model <cit.>.  In this case, incidence is precisely equal to prevalence, and thus either method for incorporating historic QoI data should yield an equivalent model.  This insight leads to the following Remark:

 
Using incidence QoI data in place of prevalence QoI data naively leads to an SIR curve where the I compartment – which normally corresponds to prevalence – now models the progression of incidence.  Using historic incidence QoI data in the way outlined in this manuscript uses the incidence curve to model incidence as is desired.  While either can be viable for the purposes of forecasting, only the method described in this paper leads to rate parameter estimates that can be interpreted as infection and recovery rates. 




§ APPLICATION TO SEASONAL INFLUENZA DATA
 
We recreate the data application from <cit.>, using the updated model from Section <ref>.  The aim of this application is not to improve the forecasting in the original formulation of the DBSSM (see Remark <ref>).  Rather, we will demonstrate that the forecasting capabilities of this model remain the same, while we also observe different estimates for the infection rate β, the recovery rate γ and the basic reproduction number ρ.

The source data modified and then used for this application are counts of patients seen in the US with an influenza-like illness (ILI), where ILI is defined as having a temperature of at least 100 degrees Fahrenheit, a cough and/or a sore throat, and no known cause for those symptoms other than influenza <cit.>.  These data are collected weekly, where more than 3400 outpatient healthcare providers report to the CDC the number of patients with ILI they treated <cit.>.

The number of patients reported as having ILI will naturally also include cases of respiratory illnesses other than influenza.  Following the approach of <cit.>, we use virologic surveillance data (where patients are actually tested for influenza) to estimate the proportion of ILI patients with influenza, then multiply ILI data by this proportion.  This corrected data is referred to as ILI+.  For more details on this adjustment, see <cit.>. Note that the ILI+ data estimates the weekly incidence of influenza cases – not prevalence.



To fit both versions of the DBSSM, we use ten influenza seasons: the seasons that started in the years 2002-2007, and the seasons the started during 2010-2013.  The years 2008 and 2009 were omitted to be consistent with <cit.>; these two years correspond to a pandemic and the focus of that work was to forecast seasonal influenza.  Each season is defined as 35 consecutive weeks starting on roughly the first week of October (epidemiology week 40, treated as t = 1).

For a estimated proportion of individuals in a population infected with influenza at timepoints 𝕋 = { 1, …, T}, suppose only the ILI+ data up through t' ∈𝕋 are observed.  Given this, we simulate 62500 from the posterior π(θ_1:t', ϕ | y_1:t' ) for four separate chains, discarding the first 12500 as burn-in and thinning out all but every tenth observation in the remaining samples.  Given these draws from the posterior distribution, the posterior predictive density, π(y_(t'+1):T | y_1:t'), is used to estimate “future" observations of ILI+ data.  

We perform two separate fits of this posterior model, on the first thirteen days (t'=13) and on the first twenty two days (t' = 22), for the considered ILI+ data for the 2010 influenza season in the United States.  These fits are performed both using the original formulation of the DBSSM, which naively uses incidence data directly in place of prevalence data, and using the new formulation developed in this paper, which uses the new maps developed in Section <ref> for a more principled treatment of incidence data.  These fits and forecasts are outlined in Figure <ref>.  The dark shaded grey regions prior to t' mark the 95 percentiles of the posterior density, while the lighter grey shaded regions after t' make the 95% prediction intervals.  The forecast using incidence data up until t'=13 has a narrower prediction interval than the one using prevalence data; each of the forecasts that use data up through t'=22 are laregly comparable.



In addition to the slight improvements on forecasting we observe in Figure <ref>, this method also has strong implications for the interpretability of β, γ for the fitted model.  Indeed, only the updated version of the DBSSM developed in this paper leads to realizations of these parameters that can accurately be interpreted as the infection rate (β), recovery rate (γ), and the basic reproduction number (S_0 ρ).  



§ DISCUSSION

The main contribution of this paper is the development of methods to map the time and value of peak incidence to the SIR curve parameters, and vice versa, for the purpose of forecasting tasks and inference on disease rate parameters.  We do this by computationally solving a system of equations ((<ref>) and (<ref>)).  There are several impactful uses of these maps in the context of previous literature.  First, much like how the peak prevalence value (PPV) and time (PPT) are useful for public health response to an epidemic <cit.>, the analogous quantities for incidence are also useful, since they describe the influx new patients entering the hospital system on a given day.  Second, this work improves upon existing work that uses historical prevalence data to model epidemics by creating a map from PIT and PIV to the SIR parameters, since incidence is typically the data that is available for ongoing epidemics <cit.>.  In the case of the application in <cit.>, where incidence data were used in place of prevalence data without justification, we have shown that this leads to biased SIR parameter estimates (see Figure <ref>).  Furthermore, our results indicate that forecasts performed using the erroneous data specification leads to larger prediction intervals earlier on in the outbreak, although this forecast is largely comparable for the correct specification later on in the outbreak (see Figure <ref> and Remark <ref>).  Of course, it remains more correct in principle to use incidence data appropriately when fitting compartment models for forecasting with ongoing incidence data <cit.>.  We have provided a modeling framework that incorporates this data appropriately (see Section <ref>).   Lastly, while the methods discussed in <cit.> do incorporate incidence data to fit forecasting models, this paper uses a Bayesian framework and importance sampling.  Since the maps developed here are deterministic, they can be used in both a Bayesian and a Frequentist framework.

As a direction for future work, it would be useful to investigate better approximations for the solution to (<ref>).  The Taylor Approximation in Section <ref> was by far the fastest computationally, but it came with the highest error on PIT.  Finding a fast and accurate approximation to this equation would greatly increase the runtime for applications where the map between the SIR parameters and PIT/PIV must be evaluated several hundreds of thousands of times.  However, for most applications (including the one in this paper), the Compute Integral approximation is sufficiently fast.

As a second direction for future work, it would be interesting to investigate the analogous maps for more complicated compartmental models, such as the SEIR model, the SIRS model, and the SEIRH model.



§ OPEN RESEARCH SECTION

All software used to perform the simulations and studies in this paper are publicly available at <https://github.com/lanl/precog>.  



§ ACKNOWLEDGMENTS

Research presented in this article was partially supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20240066DR. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).


This research was partially funded by NIH/NIGMS under grant R01GM130668-01 awarded to Sara Y. Del Valle.



§ CONFLICT OF INTEREST


The authors have no non-financial or other financial competing interests to declare that are relevant to the content of this article other than the aforementioned declared funding sources.

imsart-nameyear 
       





§ SPECIFICATION OF THE DIRICHLET-BETA STATE-SPACE MODEL


The DBSSM is fit using a fully Bayesian framework, where the interest lies in estimating the joint posterior density of the latent space θ_1:t' and ϕ given the observed data y_1:t'.  Using the conditional independence assumption described above, this density can be written as follows:

    π(θ_1:t', ϕ | y_1:t' ) ∝ π(ϕ) π(y_1:t', θ_1:t' | ϕ) = π(ϕ) ∏_t = 1^t'ℒ (y_t | θ_t, ϕ) π(θ_t | θ_t-1, ϕ),

where π(ϕ) is some prior on ϕ, ℒ (y_t | θ_t, ϕ) is the data likelihood determined by (<ref>), and the distribution π(θ_t | θ_t-1, ϕ) is determined by (<ref>).  To perform forecasts on observations y_t' : T, where T is the final timepoint of the outbreak, one uses the posterior predictive distribution, where the model and latent-space parameters are integrated out:

    π(y_ (t'+1):T | y_1 : t' ) = ∫∫π(y_ (t'+1):T, θ_1:T, ϕ | y_1 : t' ) dθ_1:T, dϕ.


To complete the specification of the DBSSM, one must determine what priors to put on the model parameters in ϕ.  It is here that the authors in <cit.> directly address the aforementioned stability issues with fitting an SIR curve with early pandemic data.  Define the latent variable z = (PPT, PPV), and let ϕ be the expanded model parameter vector that includes z.  We factorize the new prior distribution on ϕ to get

    π(ϕ) = π(ι) π(λ | ι) π(θ_0 | λ, ι) π(z | θ_0, λ, ι) π(β, γ | z, θ_0, λ, ι).


Several modeling assumptions on this conditional distribution give the abbreviated form,

    π(ϕ) = π(ι) π(λ ) π(θ_0) π(z | θ_0) π(β, γ | z, θ_0).

Specification of the individual distributions in (<ref>) is what remains to fully define the DBSSM.  The priors π(ι), π(λ ), and π(θ_0) remain unchanged and can be found in the original paper.  In Section <ref>, the priors on π(z | θ_0) and π(β, γ | z, θ_0) are described and updated, when necessary, with the theory developed in this paper.



