key: cord-0463012-yzvml7zo
authors: Narci, Romain; Delattre, Maud; Lar'edo, Catherine; Vergu, Elisabeta
title: Inference in Gaussian state-space models with mixed effects for multiple epidemic dynamics
date: 2021-09-17
journal: nan
DOI: nan
sha: d4f286c32beb670640a80b3e41819f52d4124f1e
doc_id: 463012
cord_uid: yzvml7zo

The estimation from available data of parameters governing epidemics is a major challenge. In addition to usual issues (data often incomplete and noisy), epidemics of the same nature may be observed in several places or over different periods. The resulting possible inter-epidemic variability is rarely explicitly considered. Here, we propose to tackle multiple epidemics through a unique model incorporating a stochastic representation for each epidemic and to jointly estimate its parameters from noisy and partial observations. By building on a previous work, a Gaussian state-space model is extended to a model with mixed effects on the parameters describing simultaneously several epidemics and their observation process. An appropriate inference method is developed, by coupling the SAEM algorithm with Kalman-type filtering. Its performances are investigated on SIR simulated data. Our method outperforms an inference method separately processing each dataset. An application to SEIR influenza outbreaks in France over several years using incidence data is also carried out, by proposing a new version of the filtering algorithm. Parameter estimations highlight a non-negligible variability between influenza seasons, both in transmission and case reporting. The main contribution of our study is to rigorously and explicitly account for the inter-epidemic variability between multiple outbreaks, both from the viewpoint of modeling and inference.

of mixed-effects models has rarely been used to analyse epidemic data, except in a very few studies. Among these, in ( [33] ), the dynamics of the first epidemic wave of COVID-19 in France were analysed using an ODE system incorporating random parameters to take into account the variability of the dynamics between regions. Using a slightly different approach to tackle data from multiple epidemics, [4] proposed a likelihood-based inference method using particle filtering techniques for non-linear and partially observed models. In particular, these models incorporate unitspecific parameters and shared parameters.

In addition to the specific problem of variability reflected in multiple data sets, observations of epidemic dynamics are often incomplete in various ways: only certain health states are observed (e.g. infected individuals), data are temporally discretized or aggregated, and subject to observation errors (e.g. under-reporting, diagnosis errors). Because of this incompleteness together with the non-linear structure of the epidemic models, the computation of the maximum likelihood estimator (MLE) is often not explicit. In hidden or latent variable models which are appropriate representations of incompletely observed epidemic dynamics, estimation techniques based on Expectation-Maximization (EM) algorithm can be implemented in order to compute the MLE (see e.g. [17] ). However, the E-step of the EM algorithm requires that, for each parameter value θ, the conditional expectation of the complete log-likelihood given the observed data, Q(θ), can be computed. In mixed-effects models, there is generally no closed form expression for Q(θ). In such cases, this quantity can be approximated using a Monte-Carlo procedure (MCEM, [34] ), which is computationally very demanding. A more efficient alternative is the SAEM algorithm ( [16] ), often used in the framework of mixed-effects models ( [27] ), which combines at each iteration the simulation of unobserved data under the conditional distribution given the observations and a stochastic approximation procedure of Q(θ) (see also [13] , [20] for the study and implementation of the SAEM algorithm for mixed-effects diffusion models).

In this paper, focusing on the inference for multiple epidemic dynamics, we intend to meet two objectives. The first objective is to propose a finer modeling of multiple epidemics through a unique mixed-effects model, incorporating a stochastic representation of each epidemic. The second objective is to develop an appropriate method for jointly estimating model parameters from noisy and partial observations, able to estimate rigorously and explicitly the inter-epidemic variability. Thus, the main expected contribution is to provide accurate estimates of common and epidemic-specific parameters and to provide elements for the interpretation of the mechanisms underlying the variability between epidemics of the same nature occurring in different locations or over distinct time periods. For this purpose, we extend the Gaussian state-space model introduced in ([30]) for single epidemics to a model with mixed effects on the parameters describing simultaneously several epidemics and their observations. Then, following ( [13] ) and building on the Kalman filtering-based inference method proposed in ( [30] ), we propose to couple the SAEM algorithm with Kalman-like filtering to estimate model parameters. The performances of the estimation method are investigated on simulations mimicking noisy prevalence data (i.e. the number of cases of disease in the population at a given time or over a given period of time). The method is then applied to the case of influenza epidemics in France over several years using noisy incidence data (i.e. the number of newly detected cases of the disease at a given time or over a given period of time), by proposing a new version of the filtering algorithm to handle this type of data.

The paper is organized as follows. In Section 2 we describe the epidemic model for a single epidemic, specified for both prevalence and incidence data, and its extension to account for several epidemics through a two-level representation using the framework of mixed-effects models. Section 3 contains the maximum likelihood estimation method and convergence results of the SAEM algorithm. In Section 4, the performances of our inference method are assessed on simulated noisy prevalence data generated by SIR epidemic dynamics sampled at discrete time points. Section 5 is dedicated to the application case, the influenza outbreaks in France from 1990 to 2017. Section 6 contains a discussion and concluding remarks.

2 A mixed-effects approach for a state-space epidemic model for multiple epidemics First, we sum up the approach developed in ( [30] ) in the case of single epidemics for prevalence data and extend it to incidence data (Section 2.1). By extending this approach, we propose a model for simultaneously considering several epidemics, in the framework of mixed-effects models (Section 2.2).

The epidemic model Consider an epidemic in a closed population of size N with homogeneous mixing, whose dynamics are represented by a stochastic compartmental model with d +1 compartments corresponding to the successive health states of the infectious process within the population. These dynamics are described by a density-dependent Markov jump process Z(t) with state space {0, . . . , N} d and transition rates depending on a multidimensional parameter ζ. Assuming that Z(0)/N → x 0 (0, . . . , 0) , the normalized process Z(t)/N representing the respective proportions of population in each health state converges, as N → ∞, to a classical and well-characterized ODE:

where η = (ζ, x 0 ) and b(η, ·) is explicit and easy to derive from the Q-matrix of process Z(t) (see ( [24] ), ( [30] )).

Two stochastic approximations of Z(t)/N are available: a d-dimensional diffusion process Z(t k ) with drift coefficient b(η, ·) and diffusion matrix 1 N Σ(η, ·) (which is also easily deducible from the jump functions of the density-dependent jump process, see e.g. ( [30] )), and a time-dependent Gaussian process G N (t) with small variance coefficient (see e.g. [5] ), having for expression

where g(η, t) is a centered Gaussian process with explicit covariance matrix. There is a link between these two processes: let W(t) be a Brownian motion in R d , then g(η, t) is the centered Gaussian process

with ∇ x b(η, x) denoting the matrix ( ∂bi ∂x j (η, x)) 1≤i, j≤d . In the sequel, we rely on the Gaussian process (2) to represent epidemic dynamics.

The epidemic is observed at discrete times t 0 = 0 < t 1 , · · · , < t n = T , where n is the number of observations. Let us assume that the observation times t k are regularly spaced, that is t k = k∆ with ∆ the time step (but the following can be easily adapted to irregularly spaced observation times). Setting X k := G N (t k ) and X 0 = x 0 , the model can be written under the auto-regressive AR(1) form

All the quantities in (4) have explicit expressions with respect to the parameters. Indeed, using (1) and (3), we have

Example: SIR model. As an illustrative example, we use the simple SIR epidemic model described in Figure 1 , but other models can be considered (see e.g. the SEIR model, used in Section 5).

λI/N γ In the SIR model, d = 2 and Z(t) = (S (t), I(t)) . The parameters involved in the transition rates are λ and γ and the initial proportions of susceptible and infectious individuals are

When there is no ambiguity, we denote by s and i the solution of (8). Then, the functions b(η, ·), Σ(η, ·) and σ(η,

We refer the reader to Appendix A for the computation of b(η, ·), Σ(η, ·) and σ(η, ·) in the SEIR model. Another parameterization, involving the basic reproduction number R 0 = λ γ and the infectious period d = 1 γ , is more often used for SIR models. Hence, we set η = (R 0 , d, s 0 , i 0 ) .

Observation model for prevalence data Following ( [30] ), we assume that observations are made at times t k = k∆, k = 1, . . . , n, and that some health states are not observed. The dynamics is described by the d-dimensional AR(1) model detailed in (4) . Some coordinates are not observed and various sources of noise systematically affect the observed coordinates (measurement errors, observation noises, under-reporting, etc.). This is taken into account by introducing an additional parameter µ, governing both the levels of noise and the amount of information which is available from the q ≤ d observed coordinates, and an operator B(µ) : R d → R q . Moreover, we assume that, conditionally on the random variables (B(µ)X k , k = 1, . . . , n), these noises are independent but not identically distributed. We approximate their distributions by q-dimensional Gaussian distributions with covariance matrix P k (η, µ) depending on η and µ. This yields that the observations (Y k ) satisfy

Let us define a global parameter describing both the epidemic process and the observational process,

Finally, joining (4), (9) and (10) yields the formulation (for both epidemic dynamics and observation process) required to implement Kalman filtering methods in order to estimate the epidemic parameters:

Example: SIR model (continued). The available observations could be noisy proportions of the number of infectious individuals at discrete times t k . Denoting by p the reporting rate, one could define the operator B(µ) = B(p) = (0 p) and the covariance error as P k (φ) = 1 N p(1 − p)i(η, t k ) with i(η, t) satisfying (8) . The expression of P k (φ) mimics the variance that would arise from assuming the observations to be obtained as binomial draws of the infectious individuals.

Observation model for incidence data For this purpose, we have extended the framework developed in ( [30] ). For some compartmental models, the observations (incidence) at times t k can be written as the increments of a single or more coordinates, that isB(µ)(X k−1 − X k ) where, as above,B(µ) : R d → R q is a given operator and µ are emission parameters. Let us write the epidemic model in this framework. For k = 1, . . . , n, let

From (11), the following holds, denoting by I d the d × d identity matrix,

As X k−1 = k−1 l=1 ∆ l X + x 0 , (12) becomes:

To model the errors that affect the data collected (Y k ), we assume that, conditionally on (∆ k X, k = 1, . . . , n), the observations are independent and proceed to the same approximation for their distributions

Consequently, using (13) , (14) and (15) , the epidemic model for incidence data is adapted as follows:

Contrary to (4), (∆ k X, k = 1, . . . , n) is not Markovian since it depends on all the past observations. Therefore, it does not possess the required properties of classical Kalman filtering methods. We prove in Appendix B that we can propose an iterative procedure and define a new filter to compute recursively the conditional distributions describing the updating and prediction steps together with the marginal distributions of the observations from the model (16) . 

Consider now the situation where a same outbreak occurs in many regions or at different periods simultaneously. We use the index 1 ≤ u ≤ U to describe the quantities for each unit (e.g. region or period), where U is the total number of units. Following Section 2.1, for unit u, the epidemic dynamics are represented by the d-dimensional process (X u (t)) t≥0 corresponding to d + 1 infectious states (or compartments) with state space E = [0, 1] d . It is assumed that (X u (t)) t≥0 is observed at discrete times t k = k∆ on [0, T u ], T u = n u ∆, where ∆ is a fixed time step and n u is the number of observations, and that Y u,k are the observations at times t k . Each of these dynamics has its own epidemic and observation parameters, denoted φ u .

To account for intra-and inter-epidemic variability, a two level representation is considered, in the framework of mixed-effects models. First, using the discrete-time Gaussian state-space for prevalence (11) or for incidence data (16) , the intra-epidemic variability is described. Second, the inter-epidemic variability is characterized by specifying a set of random parameters for each epidemic.

1. Intra-epidemic variability Let us define X u,k := X u (t k ), X u,0 = x u,0 and ∆ k X u := X u (t k ) − X u (t k−1 ). Using (10), conditionally to φ u = ϕ, the epidemic observations for unit u are described as in Section 2.1.

For prevalence data, 1 ≤ k ≤ n u ,

(see (5) , (6) and (7) for the expressions of F k (·), A k−1 (·), T k (·) and (9) for B(·) and P k (·)).

For incidence data,

(see (14) for the expression of G k (·) and (15) forB(·) andP k (·)).

We assume that the epidemic-specific parameters (φ u , 1 ≤ u ≤ U) are independent and identically distributed (i.i.d) random variables with distribution defined as follows,

where c = dim (φ u ) and h(β, x) : R c ×R c → R c . The vector h(β, x) = (h 1 (β, x), . . . , h c (β, x)) contains known link functions (a classical way to obtain parameterizations easier to handle), β ∈ R c is a vector of fixed effects and ξ 1 , . . . , ξ U are random effects modeled by U i.i.d centered random variables. The fixed and random effects respectively describe the average general trend shared by all epidemics and the differences between epidemics. Note that it is sometimes possible to propose a more refined description of the inter-epidemic variability by including unit-specific covariates in (19) . This is not considered here, without loss of generality.

Example: SIR model (continued). Let s 0,u = S u (0) Nu and i 0,u = Iu(0) Nu where N u is the population size in unit u. The random parameter is φ u = (R 0,u , d u , p u , s 0,u , i 0,u ) and has to fulfill the constraints

To meet these constraints, one could introduce the following function h(β, x) :

In this example, we supposed that all the parameters have both fixed and random effects, but it is also possible to consider a combination of random-effect parameters and purely fixed-effect parameters (see Section 4.1 for instance).

To estimate the model parameters θ = (β, Γ), with β and Γ defined in (19) , containing the parameters modeling the intra-and inter-epidemic variability, we develop an algorithm in the spirit of ( [13] ) allowing to derive the maximum likelihood estimator (MLE).

The model introduced in Section 2.2 can be seen as a latent variable model with y = (y u,k , 1 ≤ u ≤ U, 0 ≤ k ≤ n u ) the observed data and Φ = (φ u , 1 ≤ u ≤ U) the latent variables. Denote respectively by p(y; θ), p(Φ; θ) and p(y|Φ; θ) the probability density of the observed data, of the random effects and of the observed data given the unobserved ones. By independence of the U epidemics, the likelihood of the observations y u = (y u,1 , . . . , y u,nu ) is given by:

Computing the distribution p(y u ; θ) of the observations for any epidemic u requires the integration of the conditional density of the data given the unknown random effects φ u with respect to the density of the random parameters:

Due to the non-linear structure of the proposed model, the integral in (21) is not explicit. Moreover, the computation of p(y u |φ u ; θ) is not straightforward due to the presence of latent states in the model. Therefore, the inference algorithm needs to account for these specific features.

Let us first deal with the integration with respect to the unobserved random variables φ u . In latent variable models, the use of the EM algorithm ( [17] ) allows to compute iteratively the MLE. Iteration k of the EM algorithm combines two steps: (1) the computation of the conditional expectation of the complete log-likelihood given the observed data and the current parameter estimate θ k , denoted Q(θ|θ k ) (E-step); (2) the update of the parameter estimates by maximization of Q(θ|θ k ) (M-step). In our case, the E-step cannot be performed because Q(θ|θ k ) does not have a simple analytic expression. We rather implement a Stochastic Approximation-EM (SAEM, [16] ) which combines at each iteration the simulation of unobserved data under the conditional distribution given the observations (S-step) and a stochastic approximation of Q(θ|θ k ) (SA-step).

a) General description of the SAEM algorithm Given some initial value θ 0 , iteration m of the SAEM algorithm consists in the three following steps:

(S-step) Simulate a realization of the random parameters Φ m under the conditional distribution given the observations for a current parameter θ m−1 denoted p(·|y; θ m−1 ).

(M-step) Update the parameter estimate by maximizing Q m (θ)

In our case, an exact sampling under p(·|y; θ m−1 ) in the S-step is not feasible. In such intractable cases, MCMC algorithms such as Metropolis-Hastings algorithm can be used ( [26] ).

In the sequel, we combine the S-step of the SAEM algorithm with a MCMC procedure.

For a given parameter value θ, a single iteration of the Metropolis-Hastings algorithm consists in:

To compute the rate of acceptation of the Metropolis-Hastings algorithm in (22), we need to calculate p(y u |φ u ; θ) = p(y u,0 |φ u ; θ) nu k=1 p(y u,k |y u,0 , . . . , y u,k−1 , φ u ; θ), 1 ≤ u ≤ U.

Let y u,k:0 := (y u,0 , . . . , y u,k ), k ≥ 1. In both models (17) and (18), the conditional densities p(y u,k |y u,k−1:0 , φ u ; θ) are Gaussian densities. In model (17) involving prevalence data, their means and variances can be exactly computed with Kalman filtering techniques (see ( [30] )). In model (18), the Kalman filter can not be used in its standard form. We therefore develop an alternative filtering algorithm.

From now on, we omit the dependence in u and Φ for sake of simplicity.

Prevalence data Let us consider model (11) and recall the successive steps of the filtering developed in ( [30] ). Assume that X 0 ∼ N d (x 0 , T 0 ) and setX 0 = x 0 ,Ξ 0 = T 0 . Then, the Kalman filter consists in recursively computing for k ≥ 1:

Incidence data Let us consider model (16) . Assume that L(

Then, at iterations k ≥ 1, the filtering steps are:

The equations are deduced in Appendix B, the difficult point lying in the prediction step, i.e. the derivation of the conditional distribution L(∆ k+1 X|Y k , · · · , Y 1 ).

Generic assumptions guaranteeing the convergence of the SAEM-MCMC algorithm were stated in ( [26] ). These assumptions mainly concern the regularity of the model (see assumptions (M1-M5)) and the properties of the MCMC procedure used in step S (SAEM3'). Under these assumptions, and providing that the step sizes (α m ) are such that ∞ m=1 α m = ∞ and ∞ m=1 α 2 m < ∞, then the sequence (θ m ) obtained through the iterations of the SAEM-MCMC algorithm converges almost surely toward a stationary point of the observed likelihood.

Let us remark that by specifying the inter-epidemic variability through the modeling framework of Section 2.2, our approach for multiple epidemics fulfills the exponentiality condition stated in (M1) provided that all the components of φ u are random. Hence the algorithm proposed above converges almost surely toward a stationary point of the observed likelihood under the standard regularity conditions stated in (M2-M5) and assumption (SAEM3').

First, the performances of our inference method are assessed on simulated stochastic SIR dynamics. Second, the estimation results are compared with those obtained by an empirical two-step approach.

For a given population of size N and given parameter values, we use the Gillespie algorithm ( [23] ) to simulate a twodimensional Markov jump process Z(t) = (S (t), I(t)) . Then, choosing a sampling interval ∆ and a reporting rate p, we consider prevalence data (O(t k ), k = 1, . . . , n) simulated as binomial trials from a single coordinate of the system I(t k ).

Model Recall that the epidemic-specific parameters are φ u = R 0,u , d u , p u , s 0,u , i 0,u . In the sequel, for all u ∈ {1, . . . , U}, we assume that R 0,u > 1 and 0 < p u < 1 are random parameters. We also set s 0,u + i 0,u = 1 (which means that the initial number of recovered individuals is zero), with 0 < i 0,u < 1 being a random parameter. Moreover, we consider that the infectious period d u = d > 0 is a fixed parameter since the duration of the infectious period can reasonably be assumed constant between different epidemics. It is important to note that the case study is outside the scope of the exponential model since a fixed parameter has been included. We refer the reader to Appendix C for implementation details.

Four fixed effects β ∈ R 4 and three random effects ξ u = (ξ 1,u , ξ 3,u , ξ 4,u ) ∼ N 3 (0, Γ) are considered. Therefore, using (19) and (20), we assume the following model for the fixed and random parameters:

In other words, random effects on (R 0 , p, i 0 ) and fixed effect on d are considered. Moreover, these random effects come from a priori independent sources, so that there is no reason to consider correlations between ξ 1,u , ξ 3,u and ξ 4,u , and we can assume in this set-up a diagonal form for the covariance matrix

Parameter values We consider two settings (denoted respectively (i) and (ii) below) corresponding to two levels of inter-epidemic variability (resp. high and moderate). The fixed effects values β are chosen such that the intrinsic stochasticity of the epidemic dynamics is significant (a second set of fixed effects values leading to a lower intrinsic stochasticity is also considered; see Appendix D for details). where CV φ stands for the coefficient of variation of a random variable φ. Let us note that the link between φ u and (β, ξ u ) for p and i 0 does not have an explicit expression.

The population size is fixed to N u = N = 10, 000. For each U ∈ {20, 50, 100}, J = 100 data sets, each composed of U SIR epidemic trajectories, are simulated. Independent samplings of φ u, j = R 0,u , d u , p u , i 0,u j , u = 1, . . . U, j = 1, . . . , J, are first drawn according to model (23) . Then, conditionally to each parameter set φ u, j , a bidimensionnal Markov jump process Z u, j (t) = (S u, j (t), I u, j (t)) is simulated. Normalizing Z u, j (t) with respect to N u and extracting the values of the normalized process at regular time points t k = k∆, k = 1, . . . , n u, j , gives the X u,k, j = S u,k, j Nu , Iu,k, j Nu 's. A fixed discretization time step is used, i.e. the same value of ∆ is used to simulate all the epidemic data. For each epidemic, T u, j is defined as the first time point at which the number of infected individuals becomes zero. Two values of ∆ are considered (∆ ∈ {0.425, 2}) corresponding to an average number of time-point observations n j = 1 U U u=1 n u, j ∈ {20, 100}. Only trajectories that did not exhibit early extinction were considered for inference. The theoretical proportion of these trajectories is given by 1 − (1/R 0 ) I0 ( [1] ). Then, given the simulated X u,k, j 's and parameters φ u, j 's, the observations Y u,k, j are generated from binomial distributions B(I u,k, j , p u, j ). The results show that all the point estimates are close to the true values (relatively small bias), whatever the interepidemic variability setting, even for small values ofn and U. When the number of epidemics U increases, the standard error of the estimates decreases, but it does not seem to have a real impact on the estimation bias. Besides, observations of higher frequency of the epidemics (largen) lead to lower bias and standard deviations. It is particularly marked concerning both expectation and standard deviations of the random parameters R 0,u and p u . Irrespective Table 1 Estimates for setting (i): high inter-epidemic variability. For each combination of (n, U) and for each model parameter (defined in the first line of the table), point estimates and precision are calculated as the mean of the J = 100 individual estimates and their standard deviations (in brackets).

Parameters Parameters to the level of inter-epidemic variability, the estimations are quite satisfactory. While standard deviations of R 0,u are slightly over-estimated, even for large U and n, this trend in bias does not affect the standard deviations of p u and i 0,u .

For a given data set, Figure 2 displays convergence graphs of the SAEM algorithm for each estimates of model parameters in setting (i) with U = 100 andn = 100. Although the model does not belong to the curved exponential family, convergence of model parameters towards their true value is obtained for all parameters. 

The inference proposed method (referred to as SAEM-KM) is compared to an empirical two-step approach not taking into account explicitly mixed effects in the model. For that purpose, let us consider the method presented in ( [30] ) (referred to as KM) performed in two steps: first, we compute the estimatesφ u independently on each of the U trajectories. Second, the empirical mean and variance of theφ u 's are computed. We refer the reader to Appendix C for practical considerations on implementation of the KM method.

Let us considern = 50 and U ∈ {20, 100}. Figure 3 displays the distribution of the bias of the parameter estimates φ u, j , j ∈ {1, . . . , J}, J = 100, obtained with SAEM-KM and KM for simulation settings (i) and (ii). We notice a clear advantage to consider the mixed-effects structure. Overall, the results show that SAEM-KM outperforms KM. This is more pronounced for standard deviation estimates in the large inter-epidemic variability setting (i) than in the moderate inter-epidemic variability setting (ii). Concerning the expectation estimates, their dispersion around the median is lower for KM than for SAEM-KM, especially in setting (ii), but the bias of KM estimates is also higher. When the inter-epidemic variability is high (setting (i)), the performances of the two inference methods are substantially different. In particular, KM sometimes fails to provide plausible estimates (especially for parameter R 0 ).

We also tested other values forn and N (not shown here), e.g.n = 20 (lower amount of information) and N = 2000 (higher intrinsic variability of epidemics). In such cases, KM also failed to provide satisfying estimations whereas the mixed-effects approach was much more robust.

Data The SAEM-KM method is evaluated on a real data set of influenza outbreaks in France provided by the Réseau Sentinelles (url: www.sentiweb.fr). We use the daily number of influenza-like illness (ILI) cases between 1990-2017, considered as a good proxy of the number of new infectious individuals. The daily incidence rate was expressed per 100, 000 inhabitants. To select epidemic periods, we chose the arbitrary threshold of weekly incidence of 160 cases per 100, 000 inhabitants ( [7] ), leading to 28 epidemic dynamics. Two epidemics have been discarded due to their bimodality (1991-1992 and 1998 ). Therefore, U = 26 epidemic dynamics are considered for inference.

Compartmental model Let us consider the SEIR model (see Figure 4 ). An individual is considered exposed (E) when infected but not infectious. Denote η = (λ, , γ, x 0 ), with x 0 = (s 0 , e 0 , i 0 , r 0 ), the parameters involved in the transition rates, where is the transition rate from E to I. ODEs of the SEIR model are as follows: where the population size is fixed at N u = N = 100, 000. Denote by Inc u (t k ) the number of newly infected individuals at time t k for epidemic u. We have

Observations are modeled as incidence data observed with Gaussian noises. We draw our inspiration from ( [3] ) to account for over-dispersion in data. Therefore, assuming a reporting rate p u for epidemic u, the mean and the variance of the observed newly infected individuals are respectively defined as p u Inc u (t k ) and p u Inc u (t k ) + τ 2 u p 2 u Inc u (t k ) 2 , where parameter τ u is introduced to handle over-dispersion in the data. Denote φ u = η u , p u , τ 2 u . Therefore, we use the model defined in (18) 

where x(·, t) is the ODE solution of (24). We consider that the basic reproduction number R 0 and the reporting rate p are random, reflecting the assumptions that the transmission rate of the pathogen varies from season to season and the reporting could change over the years. Moreover, we assume e u (0) = i u (0) random and unknown (i.e. the proportion of initial exposed and infectious individuals is variable between epidemics). [7] assumed that at the start of each influenza season, a fixed average of 27% of the population is immune, that is r 0,u = r 0 = 0.27. To assess the robustness of the model with respect to the r 0 value, we test three values: r 0 ∈ {0.1, 0.27, 0.5}. This leads to s 0,u = 1 − r 0 − 2i 0,u random and unknown. Finally, we assume that τ 2 u = τ 2 is fixed and unknown. To sum up, we have to consider in the model: known parameters (d E , d I ) ∈ {(0.8, 1.8), (1.6, 1.0), (1.9, 4.1)} and r 0 ∈ {0.1, 0.27, 0.5}; fixed and unknown parameter τ 2 ; random and unknown parameters R 0 , i 0 and p.

Therefore, using (19) , we consider the following model for random parameters:

where fixed effects β ∈ R 4 and the random effects are ξ u ∼ i.i.d. N 3 (0, Γ) with Γ a covariance matrix assumed to be diagonal.

Parameter estimates We consider nine models with different combinations of values of ((d E , d I ), r 0 ). Using importance sampling techniques, we estimate the observed log-likelihood of each model from the estimated parameters values initially obtained with the SAEM algorithm. Table 3 provides the estimated log-likelihood values of the nine models of interest. Irrespectively of the r 0 value, we find that the model with (d E , d I ) = (1.9, 4.1) outperforms the two other models in terms of log-likelihood value. Moreover, for a given combination of values of (d E , d I ), the estimated log-likelihood values are quite similar according to the three r 0 tested values. Let us focus on the model with (d E , d I ) = (1.9, 4.1). Table 4 presents the estimation results of the model parameters obtained by testing the three values of r 0 : 0.1, 0.27 and 0.5. Table 4 Estimates of the mean, 5th and 95th percentiles and coefficient of variation (CV) for model parameters R 0,u , i 0,u , p u , τ 2 , assuming (d E , d I ) = (1.9, 4.1) and testing three values of r 0 : 0.1, 0.27 and 0.5. For fixed parameter, only the estimated mean is available. The average estimated value of R 0 is quite contrasted according to the r 0 value: between 1.81 and 3.28 from r 0 = 0.1 to r 0 = 0.5. By comparison, in ( [7] ), R 0 is estimated to be 1.7 during school term, and 1.4 in holidays, using a population structured into households and schools. [8] estimated a different reproduction numberR = (1−r 0 )R 0 = 1.3, measuring the transmissibility at the beginning of an epidemic in a partially immune population, from mortality data. In our case, the average value ofR is estimated to 1.63, 1.63 and 1.64 when r 0 = 0.1, 0.27 and 0.5 respectively. Therefore, given the nature of the observations (new infected individuals) and the considered model, this appears to be difficult to correctly identify R 0 together with r 0 . Indeed, the fraction of immunized individuals at the beginning of each seasonal influenza epidemic is an important parameter for the epidemic dynamics, but its value is not well known. This has implications for the stability of the estimation of the other parameters. Interestingly, the average reporting rate is estimated particularly low (around 10% irrespective of the r 0 value). Moreover, we observe that R 0 together with p and i 0 seem to be variable from season to season, with moderate coefficient of variation CV(R 0,u ) close to 15% and high coefficients of variation CV(p u ) and CV(i 0,u ) around 50% and 70% respectively. The post-predictive check is shown in Figure 5 . The difference between the average simulated curves obtained with estimated parameter values is negligible according to the r 0 value. Considering the values ofR, very close in the three scenarios, the proximity of the predicted trajectories is not surprising. Let us emphasize that the majority of the observations are within the predicted envelope (5th and 95th percentiles). Moreover, the predicted average trajectory informs about generic trends of influenza outbreaks: on average, the epidemic peak should be reached around 25 days after the beginning of the outbreak with an incidence of 90/100, 000 inhabitants approximately.

In this paper, we propose a generic inference method taking into account simultaneously in a unique model multiple epidemic trajectories and providing estimations of key parameters from incomplete and noisy epidemic data (prevalence or incidence). The framework of the mixed-effects models was used to describe the inter-epidemic variability, whereas the intra-epidemic variability was modeled by an autoregressive Gaussian process. The Gaussian formulation of the epidemic model for prevalence data used in ( [30] ) was extended to the case where incidence data were considered. Then, the SAEM algorithm was coupled with Kalman-like filtering techniques in order to estimate model parameters.

The performances of the estimators were investigated on simulated data of SIR dynamics, under various scenarios, with respect to the parameter values of epidemic and observation processes, the number of epidemics (U), the average number of observations for each of the U epidemics (n) and the population size (N). The results show that all estimates are close to the true values (reasonable biases), whatever the inter-epidemic variability setting, even for small values ofn and U. The performances, in term of precision, are improved when increasing U, whereas the bias and standard deviations of the estimations decrease when increasingn. We also compared our method with a two-step empirical approach that processes the different data sets separately and combines the individual parameter estimates a posteriori to provide an estimate of inter-epidemic variability ( [30] ). When the number of observations is too low and/or the coefficient of variation of the random effects is high, SAEM-KM clearly outperforms KM.

The proposed inference method was also evaluated on an influenza data set provided by the Réseau Sentinelles, consisting in the daily number of new infectious individuals per 100, 000 inhabitants between 1990 and 2017 in France, using a SEIR compartmental model. Testing different combinations of values for (d E , d I ) and r 0 , we find that (d E , d I ) = (1.9, 4.1) leads to the best fitting model. Then, irrespective to the r 0 value, we estimated an average value ofR = (1 − r 0 )R 0 to be around 1.6. Moreover, we highlighted a non-negligible variability from season to season that is quantitatively assessed. This variability appears especially in the initial conditions (i 0 ) and the reporting rate (p), as a combined effect of observational uncertainties and differences between seasons. Although to a lesser extent, R 0 also appears to vary between seasons, plausibly reflecting the variability in the transmission rate (λ). Obviously, the estimations can strongly depend on the choice of the compartmental model, the nature and frequency of the observations and the distribution of the random parameters. Our contribution is to propose a finer estimation of the model parameters by taking into account simultaneously all the influenza outbreaks in France for the inference procedure. This leads to an explicit and rigorous estimation of the seasonal variability.

Other methods have been implemented to deal with multiple epidemic dynamics. [4] proposed a likelihood-based inference methods for panel data modeled by non-linear partially observed jump processes incorporating unit-specific parameters and shared parameters. Nevertheless, the framework of mixed-effects models was not really investigated.

[33] used an ODE system with mixed effects on the parameters to analyse the first epidemic wave of Covid-19 in various regions in France by inferring key parameters from the daily incidence of infectious ascertained and hospitalized infectious cases. To our knowledge, there are no published studies aiming at the estimation of key parameters simultaneously from several outbreak time series using both a stochastic modeling of epidemic processes and random effects on model parameters.

The main advantage of our method is to propose a direct access to the inter-epidemic variability between multiple outbreaks. Taking into account simultaneously several epidemics in a unique model leads to an improvement of statistical inference compared with empirical methods which consider independently epidemic trajectories. For example, we can mention two experimental settings: (1) the number of epidemics is high but the number of observations per epidemic is low; (2) the number of observations per epidemic is high but the number of epidemics is low. In such cases, mixed-effects approaches can provide more satisfying estimation results. This benefit more than compensates for the careful calibration of the tuning parameters of the SAEM algorithm.

In some practical cases in epidemiology, it might be difficult to determine whether a parameter is fixed or random. Consequently, our approach could be associated with model selection techniques to inform this choice, using a criterion based on the log-likelihood of observations (see for instance ( [14] ) and ( [15] )). This would allow to determine more precisely which parameters reflect inter-individual variability and thus help to better understand the mechanisms underlying this variability. Moreover, we presented a case study on influenza outbreaks, where the variability between epidemics is seasonal, but our approach can be also applied on epidemics spreading simultaneously in many regions. In this case, the inter-epidemic variability is spatial and it would be interesting to evaluate trends from one region to another.

This work was supported by the French Agence National de la Recherche [project CADENCE, ANR-16-CE32-0007-01] and by a grant from Région Île-de-France (DIM MathInnov).

We thank the Réseau Sentinelles (INSERM/Sorbonne Université, www.sentiweb.fr) for providing a real data set of influenza outbreaks in France.

In the SEIR model, epidemic parameters are the transition rates λ, and γ and the initial proportions of susceptible, exposed and infectious individuals s 0 = S (0) N , e 0 = E(0) N and i 0 = I(0) N . When there is no ambiguity, we denote by s, e and i respectively the solutions s(η, t), e(η, t) and i(η, t) of the system of ODEs defined in (24) . Then, the functions b(η, ·) and Σ(η, ·) are b(η, s, e, i) =

and the Cholesky decomposition of Σ(η, ·) yields σ(η, s, e, i) =

Consequently, conditionally to Y 1 , Y 2 , ∆ 1 X and ∆ 2 X are independent and

Then, the generalization to the case k ≥ 1 is direct, leading to the Kalman filter described in Section 3 for incidence data.

Let us make some remarks on practical implementation. -An extended algorithm for non-exponential models is proposed to include fixed effects (see e.g. [11] ). Let κ be a fixed parameter to be estimated. First, for m = 1, . . . , M 0 , we use the classical procedure of the SAEM algorithm, that is a mean and a variance of the parameter is estimated at each iteration as if it were a random parameter. Then, at each new iteration m+1, the current variance of the parameter, denoted ω (m+1) κ , is updated as: ω (m+1)

-Due to the small influence of the number of iterations in the Metropolis-Hastings procedure (see e.g. [27] ), a single iteration is used. Furthermore, if the proposal distribution is the marginal distribution p(Φ;θ), the expression of the acceptance probability is simplified as follows:

-A stopping criterion for the SAEM algorithm is considered. Denote by θ (m) j the j-th component of θ estimated at iteration m of the SAEM algorithm. Then, the algorithm stops either when the criterion

is satisfied several times consecutively or when a limit of M max iterations is reached. The value of µ 0 is chosen sufficiently small (e.g. of the order of 10 −3 or 10 −4 ).

-As the convergence of the SAEM algorithm can strongly depend on the initial guess, a simulated annealing version of SAEM ( [25] ) is used to escape from potential local maxima of the likelihood during the first iterations and converge to a neighborhood of the global maximum. LetΓ φ ( j) m the estimated variance of the j-th component of Φ m at iteration m of the SAEM algorithm. Then, while m ≤ M 0 ,

with 0 < τ 0 < 1. For m > M 0 , the usual SAEM algorithm is used to estimate the variances at each iteration (see e.g. [28] ).

-For the initialization of the SAEM algorithm, the starting parameter values β 0 of the fixed effects β are uniformly drawn from a hypercube encompassing the likely true values. The initial variances Γ 0 are chosen sufficiently large (1 by default).

-When the sampling intervals between observations ∆ are large, the approximation of the resolvent matrix proposed in ( [30] ), Appendix A, is used.

-Concerning the KM approach, we use the Nelder-Mead method implemented in the optim function of the R software to maximize the approximated log-likelihood given by the Kalman filter. This requires to provide some initial values for the unknown parameters. As the optimization can be very sensitive to initialisation, 10 different starting values are considered and the maximum value for the log-likelihood among them are chosen. The starting parameter values for the maximization algorithm are uniformly drawn from a hypercube encompassing the likely true values.

For simulation studies in Section 4.1, the tuning parameters values are chosen as: M 0 = 500, ν 0 = 0.6, K 0 = 0.87, µ 0 = 0.001, M max = 1000 and τ 0 = 0.98. Concerning the investigation of influenza outbreaks in Section 5, we chose: M 0 = 5000, ν 0 = 0.6, K 0 = 0.87, µ 0 = 0.0001 and τ 0 = 0.98. The algorithm stops when the criterion is checked 100 times successively.

D Estimation results for a second set of parameter values

We consider a second set of parameter values which induces a lower intrinsic variability between epidemics. As for the first set of values, we consider two settings (denoted respectively (i) and (ii)) corresponding to two levels of inter-epidemic variability (resp. high and moderate): As for the first set of parameters values, all point estimates are closed to the true values. The standard error of the estimates decreases when the number of epidemics U and the number of observationsn increases, whereas the bias is only sensitive ton (bias decreasing whenn increasing).

For a given data set, Figure 6 displays convergence graphs for model parameters in setting (i) with U = 100 and n = 100. 

Stochastic epidemic models and their statistical analysis

Assessing optimal target populations for influenza vaccination programmes: An evidence synthesis and modelling study

Modeling and inference for infectious disease dynamics: A likelihood-based approach

Panel data analysis via mechanistic models

Stochastic epidemic models with inference

Time Lines of Infection and Disease in Human Influenza: A Review of Volunteer Challenge Studies

Estimating the impact of school closure on influenza transmission from sentinel data

Seasonal influenza in the united states, france, and australia: transmission and prospects for control

Estimation for dynamical systems using a population-based kalman filterapplications to pharmacokinetics models

Estimating influenza latency and infectious period durations using viral excretion data

On the curved exponential family in the Stochatic Approximation Expectation Maximization Algorithm

Parametric inference for discrete observations of diffusion processes with mixed effects

Coupling the saem algorithm and the extended kalman filter for maximum likelihood estimation in mixed-effects diffusion models

A note on BIC in mixed-effects models

An iterative algorithm for joint covariate and random effect selection in mixed effects models

Convergence of a stochastic approximation version of the em algorithm

Maximum likelihood from incomplete data via the em algorithm

Parametric inference for mixed models defined by stochastic differential equations

A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models

Using pmcmc in em algorithm for stochastic mixed models: theoretical and practical issues

Strategies for containing an emerging influenza pandemic in southeast asia

Pandemic potential of a strain of influenza a (h1n1): early findings

Exact stochastic simulation of coupled chemical reactions

Approximation of epidemic models by diffusion processes and their statistical inference

Optimization by simulated annealing: Quantitative studies

Coupling a stochastic approximation version of em with an mcmc procedure

Maximum likelihood estimation in nonlinear mixed effects models

Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools

Transmissibility of 1918 pandemic influenza

Inference for partially observed epidemic dynamics guided by kalman filtering techniques

Mixed-Effects Models in S and S-PLUS

Initial human transmission dynamics of the pandemic (h1n1) 2009 virus in north america

Population modeling of early covid-19 epidemic dynamics in french regions and estimation of the lockdown impact on infection rate

A monte carlo implementation of the em algorithm and the poor man's data augmentation algorithms