key: cord-0951770-isqqhki4
authors: Mummert, Anna; Otunuga, Olusegun M.
title: Parameter identification for a stochastic SEIRS epidemic model: case study influenza
date: 2019-05-06
journal: J Math Biol
DOI: 10.1007/s00285-019-01374-z
sha: 36c129aaa051ff049a4d92fbf4f388e340a865c8
doc_id: 951770
cord_uid: isqqhki4

A recent parameter identification technique, the local lagged adapted generalized method of moments, is used to identify the time-dependent disease transmission rate and time-dependent noise for the stochastic susceptible, exposed, infectious, temporarily immune, susceptible disease model (SEIRS) with vital rates. The stochasticity appears in the model due to fluctuations in the time-dependent transmission rate of the disease. All other parameter values are assumed to be fixed, known constants. The method is demonstrated with US influenza data from the 2004–2005 through 2016–2017 influenza seasons. The transmission rate and noise intensity stochastically work together to generate the yearly peaks in infections. The local lagged adapted generalized method of moments is tested for forecasting ability. Forecasts are made for the 2016–2017 influenza season and for infection data in year 2017. The forecast method qualitatively matches a single influenza season. Confidence intervals are given for possible future infectious levels.

weather, social behavior, and strain-specific factors. While some disease parameters, such as the latency period, can be estimated from laboratory experiments or surveys of infected populations, the transmission rate is difficult to estimate. In this work, we describe a new parameter identification technique that can be used to determine the transmission rate given infection data.

Parameter identification studies the inverse problem 'under what conditions can observations of a modeled system be used to identify the value of model parameters?'. Focusing on an epidemic model for influenza, 'how can infection data be used to identify an unknown transmission rate?'. In epidemic models, the observed infected data appears (almost) directly in the model, as a model state. However, the theory of parameter identification applies more broadly to observations that are instead some (possibly unknown) function of the model states. When parameters are unidentifiable, using identifiability techniques may uncover combinations of dependent parameters. Similarly, in the ideal scenario of noise-free data that is modeled with the stated model, even if the parameter values can be theoretically uniquely determined, computational limitation may prevent identifiability and may reveal dependent parameter combinations.

Epidemic S I R-type models are typically non-linear and, as such, model parameters are difficult to identify. (This is in contrast to linear models for which many analytical techniques exist (Godfrey and DiStephano 1987) ). The best known identification techniques for non-linear deterministic models and constant parameter values use Taylor series (Gunn et al. 1997; Pohjanpalo 1978) , similarity transformations (Chappell et al. 1990; Evans et al. 2002; Vajda et al. 1989) , or differential algebra (Audoly et al. 2001; Eisenberg 2013; Eisenberg et al. 2013; Ljung and Glad 1994) . For time-dependent parameters, differential algebra approaches can also be used (Hadeler 2011; Mummert 2013; Pollicott et al. 2012) , as well as modulating functions methods (method of moment functionals) (Ungarala et al. 2013) .

Parameter estimation for stochastic models, also called model calibration, is a type of statistical inference. Well-known estimation procedures for stochastic models are least squares estimation (Banks et al. 2014; Escobar 2012) , maximum likelihood estimation (Heijmans and Magnus 1986; Bishwal 2008) , Kalman filtering (linear quadratic estimation) (Cazelles and Chau 1995; Julier and Uhlmann 2004; Kalman 1960) and the generalized method of moments (Hansen 1982; Jeisman 2005; Hurn et al. 2007 ). Some of these techniques can be extended to estimate time-varying parameters, for example using recursive least squares (with a forgetting factor) (Escobar 2012) or an expectation maximization algorithm paired with Kalman filtering (Olama et al. 2009 ). We direct interested readers to Bishwal (2008) for a compilation of well-known stochastic differential equations and parameter estimation using the maximum likelihood estimate and the Bayes estimate for both continuous and discrete observations.

Parameter estimation for stochastic epidemic models has a rich history. Some early works include Bailey (1953) (and the references within), Becker (1976) using smallpox data, and Longini et al. (1982 Longini et al. ( , 1988 focusing on influenza. Becker (1989) provides a summary of the current techniques for statistical analysis of epidemic models at the time it was published.

More recently, increased computational power has allowed investigation of more complicated and more timely models. In particular, Bayesian inference aided typi-cally by Markov Chain Monte Carlo (MCMC) algorithms allows estimation of model parameter (posterior) distributions. See O'Neill (2002) for a methodological overview. (Here we focus on epidemic models, but Bayesian techniques can be applied more broadly, such as in health applications including designing clinical trials (Berry and Stangl 2018), ecological modeling (Hobbs and Hooten 2015) , and economics (Greenberg 2013)). These techniques are applicable to estimate fixed (see O'Neill and Becker 2001) or time varying parameters estimates (see Arnold and Lloyd 2018) . The Bayesian approach naturally allows for data imputation via inclusion of missing data as an additional parameter (O'Neill and Roberts 1999; Cauchemez and Ferguson 2008) . Data imputation can also be done with other techniques, such as expectation maximization (Becker 1997; Meng and van Dyk 1997) .

Modern parameter estimation techniques and increased computational power allow exploration of epidemics as they occur. See for example the dynamics and analysis of bovine spongiform encephalopathy (BSE) epidemic in cattle in Great Britain (Anderson et al. 1996) , the dynamics of the Foot and Mouth Disease (FMD) (Keeling et al. 2001) and Ferguson et al. (2001) , pH1N1 (Fraser et al. 2009; Yang et al. 2009 ), estimation analysis of pandemic risk of Middle East Respiratory Syndrome coronavirus (MERS) (Breban et al. 2013) , and Ebola (Rivers et al. 2014) . (Deterministic models can also be used to study epidemics in real-time. Examples include SARS (Lipsitch et al. 2003) , cholera (Tuite et al. 2011) , and Ebola (Fisman et al. 2014) ). See Cauchemez et al. (2006a) for proof of concept of real-time estimation techniques.

We direct readers interested in statistical epidemic models to the papers of Britton (2010) and Allen (2017) , and for modeling and parameter estimation to the books of Andersson and Britton (2000) and Becker (2015) .

In addition to identifying model parameters, it is equally important to understand and quantify the uncertainty in these parameter estimates. For simple linear regression with normally distributed errors, the standard error can be calculated directly, leading to confidence intervals for the parameter estimates. For more complicated stochastic models, several techniques exist to estimate the standard error of the resulting parameter estimates, including bootstrapping and asymptotic theory (Efron 1979; Banks et al. 2010 ). We direct interested readers to Banks (Banks et al. 2014 ) for a thorough review of uncertainty in modeling, including uncertainty propagation through time due to uncertainty in the model formulation and uncertainty in measurement error during data collection (the observation process); standard error calculations are included.

Here we describe a recent parameter identification technique for time-dependent parameters in stochastic dynamic models (Otunuga 2014; Otunuga et al. 2017 Otunuga et al. , 2019 , the local lagged adapted generalized method of moments (LLGMM), an extension of the generalized method of moments (Hansen 1982 ) (see also Sect. 3). The LLGMM was developed in the context of energy commodity spot prices, which are subject to response time delay and random environmental perturbations. The main advance of the LLGMM is to use some number of the past state values to construct local moment equations which are used estimate model parameters. In real world dynamic modeling problems, current and future states of continuous time dynamic processes can be influenced by the past state history. Therefore, the LLGMM technique is applicable to dynamic processes in the biological, financial, physical, chemical, and social sciences.

The LLGMM method is composed of the following components: (1) development of a stochastic mathematical model of continuous time dynamic process, (2) construction of an equivalent time series model, (3) development of generalized method of moment/observation equations, (4) introduction of a conceptual and computational parameter estimation scheme, (5) introduction of a conceptual and computational state estimation scheme, and (6) derivation of -best sub-optimal state and parameter estimates. In this work, we use the Monte-Carlo method and implicit Euler scheme (Kloeden and Platen 1995) for stochastic differential equation (SDE) to construct local moments equations and describe theoretical parameter estimation procedure for the SDE.

We demonstrate the local lagged adapted generalized method of moments parameter identification procedure (LLGMM) for the time-dependent transmission rate in a stochastic susceptible-exposed-infectious-recovered-susceptible (S E I RS) epidemic model, using US influenza data from 2004-2017. We assume the data is noise-free. The deterministic and corresponding stochastic S E I RS models are described in Sect. 2. The LLGMM is described in Section 3 and its application to the S E I RS model in Sect. 3.1. The identification procedure is applied to US influenza data in Sect. 4. The technique is evaluated with goodness-of-fit measures described in Sect. 4.1. The LLGMM is evaluated for forecasting ability in Sect. 4.2.

The deterministic susceptible-exposed-infectious-recovered-susceptible (S E I RS) model with vital rates (Fig. 1) is given by the equations

where S(t), E(t), I (t), R(t) are the fractions of the total population that are susceptible, infected but not infectious (exposed), infectious, and temporarily immune (temporarily recovered), respectively, at time t ≥ t 0 = 0. We assume every individual is in one of these compartments at t 0 with 1 = S 0 + E 0 + I 0 + R 0 . Adding the equations shows 1 = S(t) + E(t) + I (t) + R(t) for all t. The transmission rate is β(t), assumed to be a positive, bounded and continuous function of time t. We also assume that all susceptible individuals are exposed at the same transmission rate, β(t). The parameters α, ν, γ , and δ are each positive constants. Individuals are exposed for time 1/α, infectious for 1/ν, and temporarily immune for 1/γ . The demographic birth and death rate is δ; there are no disease related deaths. We direct readers interested in a review of epidemiological compartment modeling with systems of differential equations to Edelstein-Keshet (2005) , Kong et al. (2015) and Vynnycky and White (2010) . We assume that the transmission rate is a function of time and is influenced by a Gaussian white noise process causing the rate to fluctuate around the function β 0 (t). The function β 0 (t) corresponds to known influences such as periodic weather patterns or school terms; β 0 (t) could be constant. According to Méndez et al. (2012) , external fluctuations may be caused by variability in the number of contacts between infected and susceptible individuals and such random variations can be modeled by a white noise. Thus, we assume that external noise appears multiplicatively in the SEIRS model as follows

and is able to modify the dynamical behavior of the model. The transmission rate fluctuates rapidly so that we can model ζ(t) as a Gaussian white noise term with mean 0 and ζ(t)dt = dW (t), where W (t) is the standard Wiener process defined on a complete filtered probability space (Ω, F, (F t ) t≥0 , P), (F t 0 ) is measurable and (F t ) t≥0 is right continuous. The noise intensity, σ (t) > 0, of the white noise due to variations in infectivity is a function of time, measuring the amplitude of fluctuations in the transmission rate. We also assume that (S 0 , E 0 , I 0 , R 0 ) is F t 0 -measurable and independent of W (t) − W (t 0 ). By substituting (2) into (1), we get a stochastic differential equation, and in particular a Langevin equation. By construction, the resulting stochastic model will be an adapted (non-anticipating) process. The Itô approach on stochastic differential equation depends on Markovian and Martingale properties. These properties do not obey the traditional chain rule. Whereas, the Stratonovich approach obeys the traditional chain rule and allows white noise to be treated as a regular derivative of a Brownian or Wiener process. West et al. (1979) recommend that physical scientists avoid Itô calculus and claims that Stratonovich calculus is appropriate for Langevin equations with both internal and external noise. In their work, Moon and Wettlaufer (2014) suggests that in most areas of physical science where white noise is defined in terms of a δ-function autocorrelation, Stratonovich calculus is preferred due to the fluctuation-dissipation theorem. Also, Méndez et al. (2012) claim that any dynamics describing systems with real noise should be interpreted using Stratonovich equations. For these reasons, by substituting (2) into (1), we extend the resulting S E I RS model with vital rates and stochasticity to the Stratonovich equations and obtain

where • is the Stratonovich integral. The formation of β as in equation (2) and interpretation of the stochastic dynamic system of equations as a Stratonovich system of equations follow the work of Méndez et al. (2012) , who investigated a stochastic susceptible-infectious-susceptible (S I S) epidemic model. We direct readers interested in a review of stochastic modeling with systems of differential equations to Ladde and Ladde (2013) .

The S E I RS stochastic differential equation (3) has a unique (up to equivalence) global solution and the solution will remain within [0, 1] 4 whenever it starts from there. This follows from Witbooi (2017) , since the model presented here is a special case of the one presented there. (See also Parra et al. 2017 and Zhang and Teng 2007) .

We assume that all model parameter values are known constants, except β 0 (t) and σ (t). Also that the I (t) data is known and noise-free at discrete time points t j and can be modeled with the S E I RS model. The goal is to identify the time-dependent parameters β 0 (t) and σ (t), specifically β 0 (t j ) and σ (t j ) at each time t j when I is known.

The local lagged adapted generalized method of moments (LLGMM) builds on the generalized method of moments (GMM) (Hansen 1982) and we begin with a brief review of GMM. (See DeGroot and Schervish (DeGroot and Schervish 2011) for a review of the method of moments and other parameter estimation techniques.) Using the classical method of moments, one parameter is estimated by equating the model estimated expectation of y and the data estimated expectation of y. Assuming that y is generated by the model, the model estimated expectation is the population mean, while the data estimated expectation is the sample mean. The generalized method of moments constructs the moment conditions, also called orthogonality conditions, more generally as some function of the model parameters and data. The parameters are estimated by minimizing a certain norm of the sample averages of the moment conditions. If the model contains more than one parameter to be identified, the second moment, and higher if needed, would be used to construct the required moment conditions. Note that the minimization process in the GMM method uses the entire data set to estimate the constant parameters.

The LLGMM identifies time-dependent parameters using a limited number of past data points to form the moment conditions, not the entire data set. Thus the method is lagged. Also, the number of points used varies for each time estimate. Thus the method is local and adaptive.

Let

be an Itô stochastic model for y, with time-dependent parameter α j . The LLGMM parameter identification procedure estimates the expectation of Δy using an average of more than one past data point. The number of past data points used at a particular time point in the LLGMM is allowed to vary between 2 and some fixed maximum M. By defining m k as the local admissible sampled data observation size at time t k , the parameters and states at time t k are estimated using past data on the interval [t k−m k +1 , t k ]. The expected value of Δy is estimated from the data as

and from the stochastic model as

Here, F j−1 is the filtration process up to time t j−1 and for simplicity we assume Δt = t j − t j−1 is constant. We call M the maximum delay constant and m j the delay constant; m j may also be thought of as the observation size. We derive appropriate moment equations, for example by equating the two expectations above, then solve for the parameter estimate of α j−1 at time t j−1 using the well known integral mean value theorem (Khalili and Vasiliu 2010) . The computed time-dependent parameter depends on the number of past data points used in the average, so the parameter α j is more accurately expressed as α j,m j . In describing the LLGMM technique so far, we have at each time point, t j , an average of m j past data points. The value of the delay constant m j is determined using the following procedure. Fix m between 2 and M. At each time point, estimate α j,m in the stochastic model and use the value to generate a time series for y, say, y j,m . Compare the model output y j,m at time t j to the known value of y j , for example using sum of square distances. Since the model is stochastic, repeat the comparison some large number of times and compute the mean sum of square distances. Repeat for all m. Select the m value with mean model output closest to the known value. Record the value asm. Evaluate the state and parameter, y j,m and α j,m , atm as y j,m and α j,m , respectively.

In this section, we detail how the LLGMM technique, explained in general terms above, is used to estimate the state, parameters, and observation size for the S E I RS model. The method description is organized following the six components of the LLGMM scheme outlined in the Introduction (Sect. 1).

We use the Stratonovich-Itô conversion theorem given in Bernardi et al. (2001) and Kloeden and Platen (1995) to convert the Stratonovich dynamical model (3) to its Itô equivalent. We give the theorem below without proof.

Theorem 1 (Bernardi et al. 2001 and Platen 1995 

defined componentwise as

To apply the LLGMM to the S E I RS model, we derive the corresponding Itô equations from system (3) using Theorem 1 as

where the initial conditions S(t 0 ) = S 0 , E(t 0 ) = E 0 , I (t 0 ) = I 0 , and R(t 0 ) = R 0 are each between 0 and 1. Together they satisfy 1 = S 0 + E 0 + I 0 + R 0 . As it is well known, the Itô model (7) rests upon the Markovian and Martingale properties.

As part of the components of the LLGMM scheme, we construct the time series model from the Itô equations (7) and develop generalized method of moment equations (orthogonality conditions) used to estimate parameters β 0 (t) and σ (t).

From the Itô Lemma,

where ., . is the quadratic variation symbol. Multiplying S and the d E equation from (7) gives,

To estimate the parameters β 0 (t) and σ (t) at time t = t j , we first integrate (8) on the interval [t j−m j +1 , t j ] to give σ (t)β 0 (t)S 2 (t)I (t) dW t and generate moment equations from the integral equation. From the continuity of β 0 (t) and positivity of S(t) and I (t), it follows from the integral mean value theorem (Khalili and Vasiliu 2010) 

Discretizing the above equation and applying expectation gives

where S k = S(t k ), E k = E(t k ), I k = I (t k ), ΔE k = E k − E k−1 and β * j ≡ β j,m j is defined as the estimate of β 0 (t) at time t = t j using delay constant m j .

We use the moment equation in (9) to derive parameter estimate for β 0 (t). Let β j,m j be a parameter estimate of β 0 (t) that is estimated using m j past data points at time t = t j . It follows from (9) that

Using the variance of ΔE, the second parameter σ * j ≡ σ j,m j can be estimated similarly as

The parameter estimates for the transmission rate β 0 and noise intensity σ are welldefined. If {I k−1 } j k= j−m j +1 = {0}, then there are no infectious individuals during the time period t = t j−m j +1 and t = t j . From the discretized equations governing the disease spread (7), one can verify that with no infectious individuals there are also no exposed individuals. In particular, the disease spread has ended and at that time we no longer estimate the transmission rate. Also note, we do not have the case where {S k−1 } j k= j−m j +1 = {0}. If this happens, then it follows from the discretized equation governing S in (7) (11), if β * j = 0, then there is no disease transmission at time t j and therefore the noise intensity σ j in the transmission rate is also zero.

We assume that the infection data I t j is known at each point t j and that the data can be modeled with the S E I RS model. Using a simple discretization of the deterministic model equations (1) and the known initial conditions, it is possible to identify S t j+1 , E t j+1 , and R t j+1 . In particular, R can be determined from (1) as

Δt αΔt , and finally S = 1 − E − I − R. These values are used as the known data values in the LLGMM estimation.

When using a traditional explicit forward Euler-type discretization to construct solutions to the S E I RS model, we encountered negative solutions. Such unrealistic negative solutions is a common numerical artifact in solutions of systems with stochasticity. Schurz (1996) developed an implicit Euler method to ensure nonnegativity of solutions to SDE. Solutions to the S E I RS model are nonnegative (Sect. 2). Therefore, to estimate the state value at time t j , we discretize the Itô equations using the implicit Euler scheme described in Schurz (1996) . (See also Cyganowski and Grune 2001.) Namely, for a stochastic model with time-dependent parameters Θ, using an Itô interpretation,

the implicit Euler scheme gives that y j+1 satisfies

for 0 ≤ ε ≤ 1, where j = 0, 1, 2, . . . , N for sample size N , l = 1, 2, . . . , L for L simulations in the Monte-Carlo method. For the S E I RS model, Θ l j ≡ Θ l j,m j represents {β * j , σ * j }, the transmission rate and noise intensity described in (10) and (11), respectively. The simulated state for the l-th simulation at time t j using m j past data set is y l j ≡ y l j,m j . The solution y l j+1 can be determined iteratively. A value of ε = 1 is used when solving I ; a value of ε = 0.9 is used for S, E, and R.

Finally, there is a need to find, among the estimated values {y l j ≡ y l j,m j } M m j =2 at time t j , the value closest to the known real value. Specifically, the correct delay constant m j must be determined. Let S l j,m j , E l j,m j , I l j,m j and R l j,m j be the simulated susceptible, exposed, infectious and recovered estimate for the l-th simulation at time t j using m j past data set. We take the average

as the simulated value of S(t), E(t), I (t) and R(t) at time t j using m j past data set. Define

as the quadratic mean square error between the known data {S t j , E t j , I t j , R t j } and the averaged realizations {S j,m j , E j,m j , I j,m j , R j,m j } M m j =2 . For any arbitrary small positive number and for each time t j , we define the -sub-optimal admissible set of m j at t j as

If m j ∈ M j gives the minimum value for Ξ j,m j , then we record m j asm j . If condition (14) is not met at time t j , then the value of m j where the minimum min m j Ξ j,m j is attained is recorded asm j . The -best sub-optimal parameter and state estimates at time t j are now recorded as Θ j,m j = {β j,m j , σ j,m j } and {S j,m j , E j,m j , I j,m j , R j,m j }, respectively.

We Our goal is to use the CDC data in the S E I RS model as a proxy for the infectious class I , however, there is an obvious mismatch between the number of infectious and the number of positive tests. We make several adjustments to the data to account for the main differences-underreporting by the laboratories and underreporting due to people who do not seek health care. The details are described in the next paragraph. Following the general adjustment procedure in Bresee et al. (2013) , we multiply the number of positive tests by 2.74 to account for underreporting by the laboratories. Additionally, we divide by 44.1556% to account for individuals not seeking medical care. This value is the weighted average of the percentage of individuals who seek medical care by age group (weighted by percent of each age group in the U.S. population) as listed in Bresee et al. (2013) . Finally, the data is scaled so that the total attack rate for each season is within known estimates. We use the weighted average of the estimates computed by Molinari et al. (2007) , by age group. The resulting attack rate is 8.44%. (This matches well with WHO estimates of 5-10% of adults and 20-30% of children are infected with influenza each year worldwide, 2 when these estimates are weighted by US population age group, the attack rate is estimated at 8.867-15.156%.) After the underreporting and non-medical seeking adjustments, we determine the seasonal attack rates. A scale factor is determined to guarantee that the smallest attack rate is at the minimum 8.44% rate. The scale factor is determined by dividing the estimated attack rate (8.44% of the total population) by the minimum number infected over one season. This scale factor is applied to all influenza data; in particular, the adjusted data is multiplied by 231.53. We assume the total U.S. population is fixed at 311,749,110 people. The adjusted number of positive tests is divided by the total population to get the fraction of the population infectious each week (I j in the model).

The CDC influenza data is given weekly, so that Δt = 1 week. According to the CDC, 3 influenza has a latency period between 1 and 4 days, with an average of 2 days, and an infectious period of 5 to 7 days after symptoms appear and 1 day before symptoms. Therefore, with t measured in weeks, α and ν should be selected in the ranges [1.75, 7] and [0.875, 1.4] per week, respectively. The annual US birth rates for the years 2004-2016 range from 13.42 to 14.18 per 1000, while the death rates range from 8.14 to 8.39 per 1000. 4 Therefore, δ should be selected in the range [0.000157, 0.000277]. No good information is known about how long an individual remains temporarily immune from influenza after having the disease. Xu et al. (2010) , suggested that prior infection with the 1918 Spanish flu provided protection during the 2009 influenza pandemic, a span of almost 100 years. However, here we consider all influenza types as one 'disease', and given the conglomeration of virus strains, we assume it is possible for individuals to become infected each year. (In fact, it is possible to become infected more than once a year, though this is rare.) Kucharski et al. (2015) , estimate that adults 30 years or older become infected with influenza approximately twice each decade. Therefore, γ could range from around 0.000192, corresponding to 100 years temporarily immune, through to 0.04, corresponding to 1/2 year, though likely the value should be closer to 0.00385, corresponding to 5 years, than to the 100 year value. A complete list of values used in the simulation is given in Table 1 .

The scaled influenza data and simulated infectious values are shown in the top panel of Fig. 2 , assuming a maximum time delay (observation size) of M = 52 weeks (1 year).

The 'scaled influenza data' is the number of positive tests reported to the CDC for thirteen influenza season scaled to account for underreporting, individuals not seeking medical care, and to guarantee the attack rate is in known estimates, as detailed at the start of Sect. 4. The 'simulated infectious values' are computed using the data by applying the LLGMM allowing up to one year of past observations (M = 52) and assuming the data is well-described by an S E I RS model, as derived in Sect. 3.1.

The algorithm does a very good job of capturing the data. The only discrepancies lie at the peaks of the largest infections. The middle and bottom panels show the timedependent transmission rate and noise intensity, respectively. Both the transmission rate and the noise intensity show modest fluctuations each season (note the difference in scales of the two functions). Together these influence the peak and timing of the infections. The large peak in the 2009-2010 influenza season can be seen as a significant influence by the transmission rate; the noise intensity shows fluctuations only after the peak. The large peak in the 2014-2015 season can be seen as a very significant peak in the transmission rate, but no corresponding peak in the noise intensity. Similar behavior is seen with all of the delay constants M tested (Table 1) .

Time-dependent delay constant determined using the LLGMM procedure, with maximum delay M = 52 weeks (1 year). The delay constant is zero for the initial delay period of M = 52 weeks. Its minimum value after that is 2 Figure 3 shows the number of past data values used in the LLGMM procedure, as determined by the computation itself. Notice that the algorithm seems to prefer using many past data points, except directly after the two large outbreaks during the 2009-2010 and 2014-2015 seasons. This is an artifact of the high delay constant of M = 52. As the maximum allowable number of past data points is decreased, the number of large spikes in m j decreases. In particular, with a maximum delay constant M = 4 the algorithm uses the minimum value of 2 around 75% of the time, with M = 13 the value of 2 is used around 50% of the time, and with M = 26 around 25%. Each of the spikes in the determined delay constant comes around the largest peaks in infections. The larger the allowed maximum delay, the more peaks appear in the delay function.

We find the goodness-of-fit measures for the simulated infectious values following Czellar et al. (2007) . We compute the root mean square error of the simulated path R AM SE, the variability using the average median absolute deviation AM AD, and the average median bias AM B:

where y t is the data and y j t is a simulated data set for each j = 1, 2, . . . , J at the times t = 1, 2, . . . ., N . A small R AM SE indicates that the estimated parameter values can generate simulated paths similar to the actual data set. A small AM AD indicates that there is small variability among all simulated paths. The goodness-of-fit measures are computed for the LLGMM procedure with maximum delay constant M = 4, 13, 26, and 52, using J = 100 pseudo-data sets ( Table 2 ). The fit measures show that using a higher delay constant (more past data points) allows for a better fit, though visually the infectious curves show almost no differences. In order to statistically compare the estimation results using different delay constant, we estimate the statistics R AM SE, AM AD, and AM B defined in (15).

We test the efficacy of the LLGMM estimations in forecasting US influenza infection levels. We assume infection data is known up to some time, and thus S, E, and R are also known. Future data values are determined using the LLGMM technique in the following way.

Each of S, E, I , R, β 0 , and σ are estimated from the past data values iteratively using the implicit Euler method. These values depend on the random term dW t and so we simulate each 100 times and compute the averages. These averages are summed and compared with the known value 1. At each time t j , the value of m j is determined as the value with sum closest to 1.

We perform two forecasts, one assuming data through the end of the 2015-2016 influenza season and forecasting the values for the entire 2016-2017 season (Fig. 4 ) and the second assuming data through the end of 2016 and forecasting the values for weeks 1 through 29 of 2017 (Fig. 5) . The 95% confidence intervals are given around the forecasted values; the confidence intervals grow as the number of data points in the estimation increases. In both cases, the forecast shows a good qualitative fit.

The goodness-of-fit measures (Sect. 4.1) are computed for the forecasted valuesforecasting year 2017 only and forecasting the 2016-2017 season ( Table 2 ). The forecasted values begin with assuming the simulated values up to the end of 2016 or the end of the 2015-2016 influenza season, respectively, using a maximum delay constant of M = 52 weeks. The forecasted values are computed with M = 4, 13, 26, and 52, using J = 100 pseudo-data sets. The AM AD measure is 0 because the simulated infectious data is estimated using = 1. Based on the goodness-of-fit values, the forecast for the 2016-2017 influenza season is noticeably better than those for just the year 2017. However, visually no benefit is seen to a higher or lower M value in either forecast.

The local lagged adapted generalized method of moments (LLGMM), a parameter identification procedure for stochastic dynamic models, was described and demon- As noted in Sect. 2, we follow the basic formation of our stochastic S E I RS model as in Méndez et al. (2012) . In particular, the model we consider here captures only the stochastic nature of the transmission function with the Weiner process explicitly in the formation of the function β(t). This is done for two reasons. First, the transmission rate is the most difficult epidemiological parameter to estimate (Anderson and May 1991) and it is highly variable, making it of paramount importance to determine. The latency and recovery period can be measured from data and are well described by a constant value, the population mean. Second, we offer here a proof of concept that the LLGMM is a powerful tool for the inverse problem in epidemiological modeling by focusing on a simple S E I RS model applied to influenza data. Future applications of the LLGMM would include more realistic modeling of influenza outbreaks in the United States. As a case in point, vaccination uptake is highly variable and important to the modern control of influenza outbreaks, therefore influenza models should include vaccination. See for example Kong et al. (2015) who investigate the effect of vaccination on a SEIRA epidemic model with separate juvenile and adult groups. They introduce an inverse method for extracting time-dependent transmission rate from pre-and-post vaccination measles incident data.

The main technique for LLGMM involves computing the expected value of model variables by averaging the values of some number of past data points. We have considered here all reported influenza infections as infections arising from one disease, which is a gross generalization. Each year multiple influenza types, e.g. A and B, and strains, e.g. influenza A H1N1 or H3N2, circulate and are constantly mutating. By averaging the transmission patterns of several past weeks to determine a current transmission rate, we allow explicit computation of β assuming previous transmission patterns are influential on the transmission of current week. From a biological view, averaging assumes mutated strains are likely to have similar transmission rates, among other properties.

The last recorded influenza pandemic (pH1N1) occurred during the 2009-2010 season and the data (Fig. 2) clearly show the summer peak followed by the higher peak in November. In the United States the 2014-2015 influenza season was moderately severe, perhaps accounted for by a mutated (drifted) H3N2 strain that did not match the vaccine Appiah et al. (2015) . Again, this peak is noticeable in the data. These two high peaks in the fraction of infectious individuals correspond to the largest spikes in the transmission rate β, with rapid oscillations only seen for the 2009-2010 season. Extreme oscillations in the noise intensity σ occur around these two influenza seasons, though only the 2009-2010 season shows a larger than normal noise intensity. The delay constant m shows a sharp decrease around both seasons. Thus, the LLGMM procedure is able to identify that these two seasons are different than the other 'normal' influenza seasons. It is also able to capture some fundamental difference between these two extreme seasons.

Our main goal here was to describe the LLGMM technique. We use the S E I RS model and influenza data as a demonstration of the technique. Therefore, we considered only a basic epidemiological model and did only modest adjustments to the CDC influenza data. The technique shows promise for future study given the differences it is able to capture in and between the two 'non-normal' influenza seasons, even with the simple model and data adjustments. Researchers with a goal of analyzing trends in the identified transmission rate should consider a more realistic model for influenza, accounting for influential aspects such as vaccination, asymptomatic infection, effects based on age group, and changes in the total population size. More accurate initial conditions and modifications to the CDC data should also be considered.

A primer on stochastic epidemic models: formulation, numerical simulation, and analysis

Transmission dynamics and epidemiology of BSE in British cattle

Influenza activity-United States, 2014-15 season and composition of the 2015-16 influenza vaccine

An approach to periodic, time-varying parameter estimation using nonlinear filtering

Global identifiability of nonlinear models of biological systems

Standard error computations for uncertainty quantification in inverse problems: asymptotic theory vs. bootstrapping

Modeling and inverse problems in the presence of uncertainty

Estimation for an epidemic model

Analysis of infectious disease data. monographs on statistics and applied probability

Uses of the EM algorithm in the analysis of data on HIV/AIDS and other infectious diseases

Modeling to onform infectious disease control

Interhuman transmissibility of Middle East respiratory syndrome coronavirus: estimation of pandemic risk

Estimated influenza illnesses and hospitalizations averted by influenza vaccination-United States, 2012-13 influenza season

Parameter estimation in stochastic differential equations

Real-time estimates in early detection of SARS

Likelihood-based estimation of continuous-time epidemic models from time-series data: application to measles transmission in London

Adaptive dynamic modeling of HIV/AIDS epidemic using extended Kalman filter

Global identificaiton of the parameters of a nonlinear systems with specified input: a comparison of methods

Indirect robust estimation of the short-term interest rate process

Probability and statistics, 4th edn. Pearson, London Edelstein-Keshet L (2005) Mathematical models in biology

Generalizing the differential algebra approach to input-output equations in structural identifiability

Identifiability and estimation of multiple transmission pathways in cholera and waterborne disease

Time-varying parameter estimation under stochastic perturbations using LSM

Identifiability of uncontrolled nonlinear rational systems

Estimating individual and household reproduction numbers in an emerging epidemic

Pandemic potential of a strain of influenza A (H1N1): early findings

Transmission intensity and impact of control policies on the foot and mouth epidemic in Great Britain

Early epidemic dynamics of the West African 2014 Ebola outbreak: estimates derived with a simple two-parameter model

Introduction to bayesian econometrics

Identifiability of model parameters

Reparameterization of unidentifiable systems using the Taylor series approach

Parameter identification in epidemic models

Large sample properties of generalized method of moments estimators

Consistent maximum-likelihood estimation with dependent observations. The general (non-normal) case and the normal case

Bayesian models: a statistical primer for ecologists

Seeing the wood for the trees: a critical evaluation method to estimate the parameters of stochastic differential equations

Estimation of the parameters of stochastic differential equations

A new approach to linear filtering and prediction problems

Dynamics of the 2001 UK foot and mouth epidemic: stochastic dispersal in a heterogeneous landscape

An extension of the mean value theorem for integrals

The inverse method for a childhood infectious disease model with its application to pre-vaccination and post-vaccination measles data

Estimating the life course of influenza A(H3N2) antibody responses from cross-sectional data

Transmission dynamics and control of severe acute respiratory syndrome

On global identifiability for arbitrary model parameterizations

Estimating household and community transmission parameters for influenza

Statistical inference for infectious diseases

Stochastic fluctuations of the transmission rate in the susceptible-infected-susceptible epidemic model

The EM Algorithm: an old folk-song sung to a fast new tune

The annual impact of seasonal influenza in the US: measure disease burden and costs

On the interpretation of Stratonovich calculus

Studying the recovery procedure for the time-dependent transmission rate(s) in epidemic models

Stochastic differential equations for modeling, estimation and identification of mobile-to-mobile communication channels

Bayesian inference for partially observed stochastic epidemics

Inference for an epidemic when susceptibility varies

A tutorial introduction to Bayesian inference for stochastic epidemic models using Markov chain Monte Carlo methods

Ladde NS (2017) Local lagged adapted generalized method of moments and applications

Local lagged adapted generalized method of moments and applications: an innovative estimation and forecasting approach and its applications

Positivity and boundedness of solutions for a stochastic seasonal epidemiological model for respiratory syncytial virus (RSV)

System identifiability based on the power series expansion of the solution

Extracting the time-dependent transmission rate from infection data via solution of an inverse ODE problem

Modeling the impact of interventions on an epidemic of Ebola in Sierra Leone and Liberia

Numerical regularization for SDEs: construction of nonnegative solutions

Cholera epidemic in Haiti, 2010: using a transmission model to explain spatial spread of disease and identify optimal control interventions

On the estimation of time-varying parameters in continuous-time nonlinear systems

Similarity transformation approach to identifiability analysis of nonlinear compartment models

Stochastic process with non-additive fluctuations: I. Itô versus Stratonovich calculus and the effects of correlations

An SEIRS epidemic model with stochastic transmission

Structural basis of preexisting immunity to the 2009 H1N1 pandemic influenza virus

The transmissibility and control of pandemic influenza A (H1N1) virus

On a nonautonomous SEIRS model in epidemiology

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations