Inference in Gaussian state-space models with mixed effects for multiple epidemic dynamics.

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 for prevalence data, 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. Moreover, we consider here incidence data, which requires to develop a new version of the filtering algorithm. Its performances are investigated on SIR simulated epidemics for prevalence and incidence 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. 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 with a parsimonious statistical model.

Introduction

Estimation from available data of model parameters describing epidemic dynamics is a major challenge in epidemiology, especially contributing to better understand the mechanisms underlying these dynamics and to provide reliable predictions. Epidemics can be recurrent over time and/or occur simultaneously in different regions. For example, influenza outbreaks in France are seasonal and can unfold in several distinct regions with different intensities at the same time. This translates into a non-negligible variability between epidemic phenomena. In practice, this inter-epidemic variability is often omitted, by not explicitly considering specific components for each entity (population, period). Instead, each data series is analysed separately and this variability is estimated empirically. Integrating in a unique model these sources of variability allows to study simultaneously the observed data sets corresponding to each spatial (e.g. region) or temporal entity (e.g. season). This approach should improve the statistical power and accuracy of the estimation of epidemic parameters as well as refine knowledge about underlying inter-epidemic variability.

An appropriate framework is represented by the mixed-effects models, which allow to describe the variability between subjects belonging to a same population from repeated data (see e.g. Pinheiro and Bates 2000; Lavielle 2014). These models are largely used in pharmacokinetics with intra-population dynamics usually modeled by ordinary differential equations (ODE) and, in order to describe the differences between individuals, random effects on the parameters ruling these dynamics (see e.g. Collin et al. 2020). This framework was later extended to models defined by stochastic differential equations incorporating mixed effects in the parameters of these diffusion processes (Donnet and Samson 2008, 2013; Delattre and Lavielle 2013; Delattre et al. 2018). To our knowledge, the framework of mixed-effects models has rarely been used to analyse epidemic data, except in a very few studies. Among these, in Prague et al. (2020), 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 different approach to tackle data from multiple epidemics, Bretó et al. (2020) proposed models that incorporate unit-specific parameters and shared parameters and studied a likelihood-based inference method using particle filtering techniques for non-linear and partially observed models. Indeed, various ways can be investigated to describe multiple epidemics. They differ according to the purpose. Modeling dependence between regional epidemics within the same country or between successive epidemic waves requires other models. The mixed-effect approach proposed here is a first step in the analysis of the variability present across epidemics. It presents the advantage over other models that it allows to avoid the well-known “curse of dimensionality” because of its parsimony in the model 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. Dempster et al. 1977). 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, , can be computed. In mixed-effects models, there is generally no closed form expression for . In such cases, this quantity can be approximated using a Monte-Carlo procedure (MCEM, Wei and Tanner 1990), which is computationally very demanding. A more efficient alternative is the SAEM algorithm (Delyon et al. 1999), often used in the framework of mixed-effects models (Kuhn and Lavielle 2005), which combines at each iteration the simulation of unobserved data under the conditional distribution given the observations and a stochastic approximation procedure of [(see also Delattre and Lavielle (2013), Donnet and Samson (2014) for the study and implementation of the SAEM algorithm for mixed-effects diffusion models].

Data from epidemic dynamics are mostly 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) or 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). In this paper, our concern is to consider both types of data and, 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 Narci et al. (2021) for prevalence data of single epidemics to a model with mixed effects on the parameters describing simultaneously several epidemics and their observations. Then, following (Delattre and Lavielle 2013) and building on the Kalman filtering-based inference method proposed in Narci et al. (2021), we propose to couple the SAEM algorithm with Kalman-like filtering to estimate model parameters. Afterwards, in order to handle incidence data, we propose a new version of the filtering algorithm that is coupled with SAEM to estimate the parameters. The performances of the estimation method are investigated on simulations mimicking noisy prevalence data, and second noisy incidence data for SIR epidemics. The method is then applied to the case of influenza epidemics in France over several years: the underlying dynamics is described by a SEIR model and data consist of noisy incidence data from 1990 to 2017.

The paper is organized as follows. In Sect. 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 Sect. 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.

A mixed-effects approach for a state-space epidemic model for multiple epidemics

First, we sum up the approach developed in Narci et al. (2021) in the case of single epidemics for prevalence data and extend it to incidence data (Sect. 2.1). By extending this approach, we propose a model for simultaneously considering several epidemics, in the framework of mixed-effects models (Sect. 2.2).

The basics of the modeling framework for the case of a single epidemic

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 compartments corresponding to the different health states of the infectious process within the population. These dynamics are described by a density-dependent Markov jump process with state space and transition rates depending on a multidimensional parameter . Assuming that , the normalized process representing the respective proportions of population in each health state converges, as , to a classical and well-characterized ODE:where and is explicit and easy to derive from the Q-matrix of process (see Guy et al. 2015; Narci et al. 2021).

Two stochastic approximations of are available: a d-dimensional diffusion process with drift coefficient and diffusion matrix (which is also easily deducible from the jump functions of the density-dependent jump process, see e.g. Narci et al. 2021), and a time-dependent Gaussian process with small variance coefficient (see e.g. Britton and Pardoux 2020), having for expressionwhere 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 , then is the centered Gaussian processand is the resolvent matrix associated to (1)with denoting the matrix . In the sequel, we rely on the Gaussian process (2) to represent epidemic dynamics.

Remark 1

This large population framework is valid only in case of a major outbreak. It does not properly describe the beginning and the end of the epidemic outbreak (for this supercritical and subcritical, respectively, branching processes are more appropriate). We expect that this middle part of the epidemic is sufficiently well described by the approximating model to allow parameters estimation. The value does not represent the starting point of the epidemic but the time where the epidemic reaches O(N). Indeed, we just need a time and a value to derive the ODE or the Gaussian process. Moreover, in the inference method developed in the sequel, the value is unknown and estimated, and for multiple epidemics, a random effect is present for modeling the of the U epidemics.

The epidemic is observed at discrete times , where n is the number of observations. Let us assume that the observation times are regularly spaced, that is with the time step (but the following can be easily adapted to irregularly spaced observation times). Setting and , the model can be written under the auto-regressive AR(1) formAll the quantities in (4) have explicit expressions with respect to the parameters. Indeed, using (1) and (3), we haveExample: SIR model As an illustrative example, we use the simple SIR epidemic model described in Fig. 1, but other models can be considered (see e.g. the SEIR model, used in Sect. 5).

Fig. 1

SIR compartmental model with three blocks corresponding respectively to susceptible (S), infectious (I) and recovered (R) individuals. Transitions of individuals from one health state to another are governed by the transmission rate and the recovery rate , respectively

In the SIR model, and . The parameters involved in the transition rates are and and the initial proportions of susceptible and infectious individuals are . Denoting , the ODE satisfied by isWhen there is no ambiguity, we denote by s and i the solution of (8). Then, the functions , and areWe refer the reader to Appendix 1 for the computation of , and in the SEIR model. Another parameterization, involving the basic reproduction number and the infectious period , is more often used for SIR models. Hence, we set .

Observation model for prevalence data Following (Narci et al. 2021), we assume that observations are made at times , 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 observed coordinates, and an operator . Moreover, we assume that, conditionally on the random variables , these noises are independent but not identically distributed. We approximate their distributions by q-dimensional Gaussian distributions with covariance matrix depending on and . This yields that the observations satisfyLet 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 . Denoting by p the reporting rate, one could define the operator and the covariance error as with satisfying (8). The expression of 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 Narci et al. (2021). For some compartmental models, the observations (incidence) at times can be written as the increments of a single or more coordinates, that is where, as above, is a given operator and are emission parameters. Let us write the epidemic model in this framework. For , letFrom (11), the following holds, denoting by the identity matrix,As , (12) becomes:To model the errors that affect the data collected , we assume that, conditionally on , the observations are independent and proceed to the same approximation for their distributionsConsequently, using (13), (14) and (15), the epidemic model for incidence data is adapted as follows:Contrary to (4), 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 2 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).

Example: SIR model (continued) Here, and the number of new infectious individuals at times is given by . Observing a proportion p of the new infectious individuals would lead to the operator . Mimicking binomial draws, the covariance error could be chosen as where satisfies (8).

Modeling framework for multiple epidemics

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

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 , and . Using (10), conditionally to , the epidemic observations for unit u are described as in Sect. 2.1.

For prevalence data, ,[see (5), (6) and (7) for the expressions of , , and (9) for and ].

For incidence data,[see (14) for the expression of and (15) for and ].

2. Inter-epidemic variability

We assume that the epidemic-specific parameters are independent and identically distributed (i.i.d.) random variables with distribution defined as follows,where and . The vector contains known link functions (a classical way to obtain parameterizations easier to handle), is a vector of fixed effects and 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.

Remark 2

As far as inference is concerned, there is a compromise to look for between a parsimonious description of the variability between the U epidemics and a more detailed one. The set-up of mixed-effects SDE or Gaussian processes allows to describe simultaneously the stochasticity within and between epidemics. In this framework, epidemics are seen as independent and the presence of structural dependencies between the for regional epidemics or the for different periods cannot be described in this set-up. This would be conceivable but at the cost of many additional parameters.

Example: SIR model (continued) Let and where is the population size in unit u. The random parameter is and has to fulfill the constraintsTo meet these constraints, one could introduce the following function :where and .

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 Sect. 4.1 for instance).

Parametric inference

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 Delattre and Lavielle (2013) allowing to derive the maximum likelihood estimator (MLE).

Maximum likelihood estimation

The model introduced in Sect. 2.2 can be seen as a latent variable model with the observed data and the latent variables. Denote respectively by ), and 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 is given by:Computing the distribution of the observations for any epidemic u requires the integration of the conditional density of the data given the unknown random effects 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 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 . In latent variable models, the use of the EM algorithm (Dempster et al. 1977) 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 , denoted (E-step); (2) the update of the parameter estimates by maximization of (M-step). In our case, the E-step cannot be performed because does not have a simple analytic expression. We rather implement a Stochastic Approximation-EM (SAEM, Delyon et al. 1999) which combines at each iteration the simulation of unobserved data under the conditional distribution given the observations (S-step) and a stochastic approximation of (SA-step). In our case, an exact sampling under in the S-step is not feasible. In such intractable cases, MCMC algorithms such as Metropolis-Hastings algorithm can be used (Kuhn and Lavielle 2004). In the sequel, we combine the S-step of the SAEM algorithm with a MCMC procedure.

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

(S-step) Simulate a realization of the random parameters under the conditional distribution given the observations for a current parameter denoted .

(SA-step) Update according to where is a sequence of positive step-sizes s.t. and .

(M-step) Update the parameter estimate by maximizing

Computation of the S-step by combining the Metropolis-Hastings algorithm with Kalman filtering techniques

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 calculateLet , . In both models (17) and (18), the conditional densities are Gaussian densities. In model (17) involving prevalence data, their means and variances can be exactly computed with Kalman filtering techniques (see Narci et al. 2021). In model (18), the Kalman filter can not be used in its standard form. We therefore develop an alternative filtering algorithm.

Generate a candidate for a given proposal distribution q

Take where

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 Narci et al. (2021). Assume that and set . Then, the Kalman filter consists in recursively computing for : Incidence data Let us consider model (16). Assume that and . Let and . Then, at iterations , the filtering steps are: The equations are deduced in Appendix 2, the difficult point lying in the prediction step, i.e. the derivation of the conditional distribution .

Prediction:

Updating:

Marginal:

Prediction:

Updating:

Marginal:

Convergence of the SAEM-MCMC algorithm

Generic assumptions guaranteeing the convergence of the SAEM-MCMC algorithm were stated in Kuhn and Lavielle (2004). 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 are such that and , then the sequence 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 Sect. 2.2, our approach for multiple epidemics fulfills the exponentiality condition stated in (M1) provided that all the components of 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’).

There is no theoretical guarantee that the algorithm converges to a global maximum of the likelihood. It is a classical problem in statistics which concerns the majority of algorithms developed to optimize non convex functions. In practice, to prevent convergence of the algorithm to a local maximum of the likelihood, it is possible to consider different starting values for the parameters and to finally choose the set of estimated values associated with the highest likelihood value among the ones obtained with these different starting points. Nevertheless, depending on the complexity of the model and the number of observations to process, the computation time of the algorithm for a given set of starting values can be important. Therefore, the strategy adopted in the paper is to use a simulated annealing version of SAEM in order to have more flexibility in the first iterations and thus to escape more easily from potential local maxima of the likelihood at the beginning of the algorithm [cf. Appendix 3, 5th item and pages 249–252 in Lavielle (2014)]. This does not completely prevent from converging to a local optimum but we can reasonably hope to reach the global optimum by considering fewer different initializations than with a standard version of the algorithm. In practice, the a priori knowledge of specialists in the field of study, in this case epidemiologists, can help to initialize the algorithm close to the optimum.

Assessment of parameter estimators performances on simulated data

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 (Gillespie 1977) to simulate a two-dimensional Markov jump process . Then, choosing a sampling interval and a reporting rate p, we consider prevalence data simulated as binomial trials from a single coordinate of the system . We refer the reader to Appendix 5 for an assessment of the performances of our inference method on simulated incidence data.

Simulation setting

Model Recall that the epidemic-specific parameters are . In the sequel, for all , we assume that and are random parameters. We also set (which means that the initial number of recovered individuals is zero), with being a random parameter. Moreover, we consider that the infectious period 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 3 for implementation details.

Four fixed effects and three random effects are considered. Therefore, using (19) and (20), we assume the following model for the fixed and random parameters:In other words, random effects on 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 , and , 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 4 for details).where stands for the coefficient of variation of a random variable . Let us note that the link between and for p and does not have an explicit expression.

Setting (i): and corresponding to , ; ; , ; , ;

Setting (ii): and corresponding to , ; ; , ; , ;

Data simulation The population size is fixed to . For each , data sets, each composed of U SIR epidemic trajectories, are simulated. Independent samplings of , , , are first drawn according to model (23). Then, conditionally to each parameter set , a bidimensionnal Markov jump process is simulated. Normalizing with respect to and extracting the values of the normalized process at regular time points , , gives the ’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, is defined as the first time point at which the number of infected individuals becomes zero. Two values of are considered () corresponding to an average number of time-point observations . Only trajectories that did not exhibit early extinction were considered for inference. The theoretical proportion of these trajectories is given by (Andersson and Britton 2000). Then, given the simulated ’s and parameters ’s, the observations are generated from binomial distributions .

Point estimates and standard deviations for inferred parameters

Tables 1 and 2 show the estimates of the expectation and standard deviation of the mixed effects , computed from the estimations of and using functions h defined in (23), for settings (i) and (ii). For each parameter, the reported values are the mean of the parameter estimates , , and their standard deviations in brackets.

Table 1

Estimates for setting (i): high inter-epidemic variability

Parameters		ER_0,u	d	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values		1.500	2.500	0.739	0.119	0.250	0.226	0.079
n¯=20	U=20	1.580	2.584	0.688	0.126	0.335	0.193	0.078
		(0.135)	(0.293)	(0.117)	(0.024)	(0.151)	(0.051)	(0.020)
	U=50	1.574	2.538	0.704	0.122	0.359	0.201	0.079
		(0.111)	(0.220)	(0.089)	(0.019)	(0.149)	(0.030)	(0.014)
	U=100	1.583	2.564	0.700	0.124	0.385	0.199	0.081
		(0.105)	(0.210)	(0.083)	(0.015)	(0.134)	(0.023)	(0.011)
n¯=100	U=20	1.501	2.502	0.734	0.118	0.292	0.217	0.075
		(0.080)	(0.159)	(0.059)	(0.021)	(0.105)	(0.035)	(0.019)
	U=50	1.510	2.522	0.729	0.120	0.305	0.217	0.080
		(0.054)	(0.126)	(0.038)	(0.014)	(0.070)	(0.022)	(0.012)
	U=100	1.503	2.508	0.738	0.119	0.308	0.216	0.079
		(0.047)	(0.097)	(0.030)	(0.010)	(0.054)	(0.016)	(0.009)

For each combination of and for each model parameter (defined in the first line of the table, with the true value displayed in bold in the second line), point estimates and precision are calculated as the mean of the individual estimates and their standard deviations (in brackets)

Table 2

Estimates for setting (ii): moderate inter-epidemic variability

Parameters		ER_0,u	d	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values		1.500	2.500	0.777	0.109	0.125	0.143	0.049
n¯=20	U=20	1.619	2.764	0.666	0.127	0.190	0.117	0.053
		(0.120)	(0.256)	(0.099)	(0.022)	(0.106)	(0.034)	(0.014)
	U=50	1.638	2.789	0.653	0.128	0.213	0.122	0.056
		(0.103)	(0.233)	(0.087)	(0.018)	(0.099)	(0.018)	(0.010)
	U=100	1.623	2.769	0.658	0.128	0.209	0.122	0.056
		(0.081)	(0.194)	(0.075)	(0.013)	(0.056)	(0.017)	(0.007)
n¯=100	U=20	1.540	2.627	0.732	0.118	0.176	0.143	0.050
		(0.066)	(0.143)	(0.057)	(0.017)	(0.055)	(0.035)	(0.012)
	U=50	1.539	2.622	0.733	0.117	0.183	0.145	0.052
		(0.044)	(0.098)	(0.041)	(0.009)	(0.038)	(0.018)	(0.007)
	U=100	1.541	2.629	0.732	0.118	0.187	0.149	0.053
		(0.040)	(0.078)	(0.030)	(0.008)	(0.035)	(0.016)	(0.006)

For each combination of and for each model parameter (defined in the first line of the table, with the true value displayed in bold in the second line), point estimates and precision are calculated as the mean of the individual estimates and their standard deviations (in brackets)

Remark 3

Via the link functions h, the random parameters can be expressed as a function of the fixed effects and the random effects . When the link has a suitable form (for example, a log link function), it is possible to explicitly obtain the mean and variance of the ’s. When it is not the case (for example, with a logit link function), we can compute the empirical mean and variance based on simulations of the ’s. For more complex link functions (such as the logit link function), this is no longer true and the empirical mean and variance are computed via simulations of the ’s. This is the method used here.

The results show that all the point estimates are close to the true values (relatively small bias), whatever the inter-epidemic variability setting, even for small values of 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 (large ) lead to lower bias and standard deviations. It is particularly marked concerning both expectation and standard deviations of the random parameters and . Irrespective to the level of inter-epidemic variability, the estimations are quite satisfactory. While standard deviations of are slightly over-estimated, even for large U and , this trend in bias does not affect the standard deviations of and .

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

Fig. 2

Convergence graphs of the SAEM algorithm for estimates of and . Setting (i) with and . Parameter values at each iteration of the SAEM algorithm (plain blue line) and true values of model parameters (dotted red line) (color figure online)

Comparison with an empirical two-step approach

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 Narci et al. (2021) (referred to as KM) performed in two steps: first, we compute the estimates independently on each of the U trajectories. Second, the empirical mean and variance of the ’s are computed. We refer the reader to Appendix 3 for practical considerations on implementation of the KM method.

Let us consider and . Figure 3 displays the distribution of the bias of the parameter estimates , , , obtained with SAEM-KM and KM for simulation settings (i) and (ii).

Fig. 3

Boxplots (25th, 50th and 75th percentiles) of the bias of the estimates of each model parameter, with , obtained with SAEM-KM (blue boxes) and KM (red boxes). Two levels: and epidemics. Dark colours: high inter-epidemic variability [setting (i)]. Light colours: moderate inter-epidemic variability [setting (ii)]. The symbol represents the estimated mean bias. For sake of clarity, we removed extreme values from the graphical representation. This concerns only the parameter and the KM method: 37 values for (35 in setting (i), 2 in setting (ii)) and 50 values for (47 in setting (i), 3 in setting (ii)) (color figure online)

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 ).

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

Case study: influenza outbreaks in France

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 and 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 (Cauchemez et al. 2008), leading to 28 epidemic dynamics. Two epidemics have been discarded due to their bimodality (1991–1992 and 1998). Therefore, epidemic dynamics are considered for inference.

Compartmental model Let us consider the SEIR model (see Fig. 4). An individual is considered exposed (E) when infected but not infectious. Denote , with , the parameters involved in the transition rates, where is the transition rate from E to I. ODEs of the SEIR model are as follows:

Fig. 4

SEIR compartmental model with four blocks corresponding respectively to susceptible (S), exposed (E), infectious (I) and recovered (R) individuals. Transitions of individuals from one health state to another are governed by the transmission rate , the incubation rate and the recovery rate

Another parametrization exhibits the basic reproduction number , the incubation period and the infectious period . Thus, the epidemic parameters are . Let us describe the two-layer model used in the sequel.

Intra-epidemic variability For each epidemic u, let andwhere the population size is fixed at . Denote by the number of newly infected individuals at time for epidemic u. We haveObservations are modeled as incidence data observed with Gaussian noises. We draw our inspiration from Bretó (2018) to account for over-dispersion in data. Therefore, assuming a reporting rate for epidemic u, the mean and the variance of the observed newly infected individuals are respectively defined as and , where parameter is introduced to handle over-dispersion in the data. Denote . Therefore, we use the model defined in (18) with , , , , and deriving from (26) in Appendix 1, andwhere is the ODE solution of (24).

Inter-epidemic variability In this real data study, due to identifiability issues, we have to perform inference by fixing parameters (incubation period), (infectious period) and (initial proportion of removed individuals). Let us first comment on the two parameters , . Studies in the literature found discrepant values of these durations (see Cori et al. 2012 for a review), varying from 0.64 (Fraser et al. 2009) to 3.0 (Pourbohloul et al. 2009) days for the incubation period and from 1.27 (Fraser et al. 2009) to 8.0 (Pourbohloul et al. 2009) days for the infectious period. For example, Cori et al. (2012) estimated that and days on average using excretion profiles from experimental infections. In two other papers, these durations were fixed according to previous studies (e.g. Mills et al. 2004; Ferguson et al. 2005): days (Chowell et al. 2008); days (Baguelin et al. 2013). Performing a systematic review procedure from viral shedding and/or symptoms, Carrat et al. (2008) estimated to be between 1.7 and 2.0 on average. Therefore in what follows, we consider the latent and infectious periods and known and test three combinations of values: , (0.8, 1.8) and (1.9, 4.1).

We consider that the basic reproduction number 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 random and unknown (i.e. the proportion of initial exposed and infectious individuals is variable between epidemics). Cauchemez et al. (2008) assumed that at the start of each influenza season, a fixed average of of the population is immune, that is . To assess the robustness of the model with respect to the value, we test three values: . This leads to random and unknown. Finally, we assume that is fixed and unknown. Thus, we have to study nine candidate models with: known parameters and ; fixed and unknown parameter ; random and unknown parameters , and p.

Therefore, using (19), we consider the following model for random parameters:where fixed effects and the random effects are with a covariance matrix assumed to be diagonal.

Parameter estimates These nine candidate models correspond to different combinations of values of . They have exactly the same structure and the same complexity in terms of number of parameters to be estimated. After the inference is performed for each of these nine candidate models, we have to choose the best candidate values. 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 value, we find that the model with outperforms the two other models in terms of log-likelihood value. Moreover, for a given combination of values of , the estimated log-likelihood values are quite similar according to the three tested values.

Table 3

Estimated values of the observed log-likelihood of the model obtained by testing nine combinations of values of

(d_E,d_I)	r_0	Estimated log-likelihood
(0.8,1.8)	0.1	9011.752
	0.27	8827.870
	0.5	8499.452
(1.6,1.0)	0.1	9147.108
	0.27	8961.991
	0.5	8643.562
(1.9,4.1)	0.1	10270.000
	0.27	10216.260
	0.5	9905.436

Remark 4

Model comparison is usually performed by using information criteria like BIC which are defined by adding a penalty term, depending on the total number of model parameters, to the log-likelihood. The best model according to these criteria is the model that leads to the smallest criterion value. We could have used BIC to compare the nine candidate models, but as they have the same number of parameters, the penalty is useless and the comparison is fully based on the log-likelihood term of the criterion. That is why we only show the estimated log-likelihood values in Table 3. Instead of comparing values of log-likelihood, we directly compare log-likelihood values, so higher values are better.

Let us focus on the model with . Table 4 presents the estimation results of the model parameters obtained by testing the three values of : 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 , assuming and testing three values of : 0.1, 0.27 and 0.5

		R_0,u	p_u	i_0,u	τ^2
Estimated mean	r_0=0.1	1.810	0.069	0.010	0.025
	r_0=0.27	2.238	0.084	0.008	0.013
	r_0=0.5	3.281	0.119	0.006	0.037
Estimated [5th,95th] percentiles	r_0=0.1	[1.470,2.264]	[0.026,0.138]	[0.003,0.023]	—
	r_0=0.27	[1.787,2.825]	[0.031,0.169]	[0.002,0.019]	—
	r_0=0.5	[2.696,3.977]	[0.044,0.238]	[0.002,0.014]	—
Estimated CV	r_0=0.1	14 %	53 %	67 %	—
	r_0=0.27	14 %	52 %	72 %	—
	r_0=0.5	12 %	51 %	74 %	—

For fixed parameter, only the estimated mean is available

The average estimated value of is quite contrasted according to the value: between 1.81 and 3.28 from to . By comparison, in Cauchemez et al. (2008), is estimated to be 1.7 during school term, and 1.4 in holidays, using a population structured into households and schools. Chowell et al. (2008) estimated a different reproduction number , measuring the transmissibility at the beginning of an epidemic in a partially immune population, from mortality data. In our case, the average value of is estimated to 1.63, 1.63 and 1.64 when , 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 together with . 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 irrespective of the value). Moreover, we observe that together with p and seem to be variable from season to season, with moderate coefficient of variation close to and high coefficients of variation and around and respectively.

It is possible to perform a maximum a posteriori (MAP) estimation of the parameters corresponding to each period, by computing where is the parameter estimate obtained with the SAEM algorithm. We refer the reader to Appendix 6 for a graphical representation of the time-series behaviour of and , which could be interesting from an epidemiological point of view.

The post-predictive check is shown in Fig. 5. The difference between the average simulated curves obtained with estimated parameter values is negligible according to the value. Considering the values of , 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.

Fig. 5

Post-predictive check. Observations (number of ILI as proxy for new infectious for each of the U epidemics) (blue). Simulated trajectories obtained for (red), (magenta) and (green) in three steps: (i) generation of 1000 values based on estimated values of parameters; (ii) given , simulation of 1000 epidemics according to the model (18); (iii) computation of average trajectory (solid line) and 5th and 95th percentiles (dotted lines) of the 1000 simulated epidemics. Population size fixed to (color figure online)

Remark 5

We observe on Fig. 5 that the two epidemics lasting longer, corresponding to the seasons 1998–1999 and 2012–2013, tend to be above the 95% percentile near the end, which could be explained by the fact that they grow very slowly the first three weeks. Also, considering a different threshold defining the epidemic season (here taken equal to 160 cases per 100,000 inhabitants) could change the data-points considered for these two trajectories and hence their positioning with respect to the average or confidence bound trajectories. Finally, the predicted envelope of the 5th and 95th percentiles is ensured to contain, by construction, only of observations. Therefore, some observations can be found below the 5th percentile or above the 95th percentile.

Discussion

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 Narci et al. (2021) 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 () 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 of and U. The performances, in term of precision, are improved when increasing U, whereas the bias and standard deviations of the estimations decrease when increasing . 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 (Narci et al. 2021). 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 and , we find that leads to the best fitting model. Then, irrespective to the value, we estimated an average value of 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 () and the reporting rate (p), as a combined effect of observational uncertainties and differences between seasons. Although to a lesser extent, 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. Bretó et al. (2020) 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. Prague et al. (2020) 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.

This paper focuses of independent epidemics. Even in this apparently simple case, a non negligible number of technical and methodological difficulties arise. Given these difficulties and as this setting is rarely suitable in practice, there is a compromise to find for inference, between a parsimonious description of the U epidemics and a more detailed one. The set-up of mixed-effects SDE allows to describe simultaneously the within and between epidemic stochasticity. This study can be considered as the first investigation step of the multiple epidemics data set, that does not prevent a second investigation step with a more accurate description including for instance some shared parameters, unit specific parameters, a dependence structure, etc. Extensions of this work would imply modifications of the model but also important modifications of the algorithm, the latter being necessarily specific to the way the dependence is accounted for. A first strategy could be to incorporate a given dependence structure between the ’s directly in the time series equations by specifying the ’s about this way: , where is known up to parameter . This mechanically increases the number of parameters in , which may lead to high computation times when U is large. A second strategy should be to introduce a correlation between the ’s by defining the epidemic-specific parameters as with a non block diagonal covariance matrix of size . Here again, the number of parameters to be estimated can be increased significantly depending on the way the dependency is parametrized. In both cases, the simulation step of the algorithm has to be modified because the random effects can no longer be simulated independently. The modifications of the algorithm implied by these two ways of modeling the dependence between epidemics require an important additional work and are thus not considered in this paper.

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 Delattre et al. (2014) and Delattre and Poursat (2020)]. 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.

Table 5

Estimates for setting (i): high inter-epidemic variability

Parameters		ER_0,u	d	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values		3.000	3.000	0.739	0.119	1.000	0.226	0.079
n¯=20	U=20	3.085	2.889	0.758	0.111	1.477	0.205	0.075
		(0.460)	(0.205)	(0.060)	(0.016)	(0.666)	(0.036)	(0.018)
	U=50	3.152	2.926	0.761	0.111	1.509	0.199	0.075
		(0.360)	(0.170)	(0.049)	(0.011)	(0.457)	(0.025)	(0.012)
	U=100	3.116	2.904	0.765	0.111	1.517	0.200	0.077
		(0.307)	(0.152)	(0.046)	(0.008)	(0.366)	(0.018)	(0.009)
n¯=100	U=20	2.929	2.932	0.742	0.116	1.124	0.212	0.075
		(0.263)	(0.144)	(0.047)	(0.016)	(0.332)	(0.029)	(0.017)
	U=50	3.002	2.973	0.749	0.116	1.186	0.207	0.075
		(0.242)	(0.116)	(0.031)	(0.012)	(0.315)	(0.022)	(0.011)
	U=100	2.952	2.942	0.751	0.115	1.159	0.212	0.075
		(0.148)	(0.090)	(0.022)	(0.008)	(0.155)	(0.018)	(0.007)

For each combination of and for each model parameter (defined in the first line of the table), point estimates and precision are calculated as the mean of the individual estimates and their standard deviations (in brackets)

Table 6

Estimates for setting (ii): moderate inter-epidemic variability

Parameters		ER_0,u	d	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values		3.000	3.000	0.777	0.109	0.500	0.143	0.049
n¯=20	U=20	3.183	3.051	0.771	0.106	0.811	0.128	0.046
		(0.292)	(0.164)	(0.046)	(0.012)	(0.321)	(0.029)	(0.011)
	U=50	3.201	3.050	0.765	0.106	0.874	0.132	0.048
		(0.208)	(0.116)	(0.035)	(0.008)	(0.241)	(0.018)	(0.007)
	U=100	3.232	3.068	0.765	0.106	0.906	0.132	0.048
		(0.189)	(0.103)	(0.028)	(0.005)	(0.212)	(0.013)	(0.005)
n¯=100	U=20	3.037	3.051	0.770	0.110	0.563	0.135	0.046
		(0.169)	(0.100)	(0.037)	(0.012)	(0.206)	(0.026)	(0.011)
	U=50	3.064	3.055	0.764	0.110	0.632	0.139	0.048
		(0.117)	(0.080)	(0.023)	(0.009)	(0.142)	(0.016)	(0.007)
	U=100	3.059	3.057	0.768	0.110	0.619	0.141	0.048
		(0.088)	(0.061)	(0.019)	(0.005)	(0.094)	(0.013)	(0.004)

For each combination of and for each model parameter (defined in the first line of the table), point estimates and precision are calculated as the mean of the individual estimates and their standard deviations (in brackets)

Table 7

Estimates for the first set of parameters values

Parameters	ER_0,u	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values (setting (i))	1.500	0.739	0.119	0.250	0.226	0.079
U=20	1.516	0.725	0.120	0.349	0.180	0.079
	(0.076)	(0.045)	(0.021)	(0.169)	(0.039)	(0.023)
U=50	1.508	0.730	0.120	0.365	0.185	0.081
	(0.053)	(0.029)	(0.013)	(0.097)	(0.022)	(0.012)
U=100	1.517	0.731	0.120	0.373	0.184	0.081
	(0.038)	(0.024)	(0.010)	(0.059)	(0.017)	(0.009)

For each 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 individual estimates and their standard deviations (in brackets)

Table 8

Estimates for the second set of parameters values

Parameters	ER_0,u	Ep_u	Ei_0,u	sdR_0,u	sdp_u	sdi_0,u
True values (setting (i))	3.000	0.739	0.119	1.000	0.226	0.079
U=20	2.871	0.790	0.163	1.618	0.174	0.117
	(0.315)	(0.049)	(0.034)	(0.973)	(0.037)	(0.033)
U=50	2.856	0.804	0.165	1.466	0.166	0.121
	(0.169)	(0.030)	(0.019)	(0.397)	(0.025)	(0.017)
U=100	2.842	0.804	0.168	1.486	0.169	0.125
	(0.131)	(0.021)	(0.017)	(0.285)	(0.019)	(0.015)

For each 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 individual estimates and their standard deviations (in brackets)