Romain Narci
email: romain.narci@inrae.fr

Maud Delattre

Catherine Larédo

Elisabeta Vergu

Inference for partially observed epidemic dynamics guided by Kalman filtering techniques †

Keywords: Approximate maximum likelihood, Diffusion approach, Kalman filter, Measurement errors, Partially-observed Markov process, Epidemic dynamics

Despite the recent development of methods dealing with partially observed epidemic dynamics (unobserved model coordinates, discrete and noisy outbreak data), limitations remain in practice, mainly related to the quantity of augmented data and calibration of numerous tuning parameters. In particular, as coordinates of dynamic epidemic models are coupled, the presence of unobserved coordinates leads to a statistically difficult problem. The aim is to propose an easy-to-use and general inference method that is able to tackle these issues. First, using the properties of epidemics in large populations, a two-layer model is constructed. Via a diffusion-based approach, a Gaussian approximation of the epidemic density-dependent Markovian jump process is obtained, representing the state model. The observational model, consisting of noisy observations of certain model coordinates, is approximated by Gaussian distributions. Then, an inference method based on an approximate likelihood using Kalman filtering recursion is developed to estimate parameters of both the state and observational models. The performance of estimators of key model parameters is assessed on simulated data of SIR epidemic dynamics for different scenarios with respect to the population size and the number of observations. This performance is compared with that obtained using the well-known maximum iterated filtering method. Finally, the inference method is applied to a real data set on an influenza outbreak in a British boarding school in 1978.

Introduction

The interest and impact of mathematical modeling and inference methods for infectious diseases have considerably grown in recent years in a context of increasing complex models and abundant data of varying quality. Estimating the parameters governing epidemic dynamics from available data has become a major challenge, in particular from the perspective of subsequently providing reliable predictions of such dynamics. Many authors have addressed the problem of key epidemic parameter estimation based on likelihood approaches (e.g., [START_REF] Cauchemez | Likelihood-based estimation of continuous-time epidemic models from timeseries data: application to measles transmission in london[END_REF]). While estimation is quite straightforward for complete observations, this is no longer true in the incomplete observation setting which occurs in practice, regardless of the mathematical formalism used. Indeed, available data tends to be only partially observed (e.g., certain health statuses such as asymptomatic infected stages cannot be observed at all; infectious and recovery dates are not observed for all individuals during the outbreak; not all infectious individuals are reported) and may also be temporally and/or spatially aggregated. Various approaches have been developed to deal with these types of data (e.g., see [START_REF] O'neill | Introduction and snapshot review: Relating infectious disease transmission models to data[END_REF], [START_REF] Britton | Introduction to statistical inference for infectious diseases[END_REF] for reviews). In the general framework of partially-observed Markov processes, some of these methods have been implemented in the R package POMP [START_REF] King | Statistical inference for partially observed markov processes via the r package pomp[END_REF]). Among these, we cite maximum iterated filtering (MIF: [START_REF] Ionides | Inference for nonlinear dynamical systems[END_REF], [START_REF] Ionides | Inference for dynamic and latent variable models via iterated, perturbed bayes maps[END_REF]) in which the parameter space is explored by considering that parameters follow a random walk over time with variance decreasing over filtering iterations, and the likelihood being stochastically estimated. Theoretical justification for convergence to the maximum likelihood estimates in the parameter space has been provided for this method [START_REF] Ionides | Iterated filtering[END_REF]). Furthermore, likelihood-free methods, such as approximate Bayesian computation based on sequential Monte Carlo (ABC-SMC, [START_REF] Sisson | Sequential monte carlo without likelihoods[END_REF], [START_REF] Toni | Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems[END_REF]) and particle Markov chain Monte Carlo (PMCMC, [START_REF] Andrieu | Particle markov chain monte carlo methods[END_REF]), have opened some of the most promising pathways for improvement. Nevertheless, these algorithms do not provide a definitive solution to statistical inference from incomplete epidemic data. Indeed, there are real limitations in practice due to the amount of augmented data and fitting the numerous tuning parameters involved. That can lead to substantial computational overheads.

In this paper, we consider a different approach to deal with the presence of missing coordinates, discrete observations, and reporting and measurement errors. Our goal is to propose a useful and coherent latent variable model that allows key epidemic parameters to be estimated from imperfect observations from outbreaks.

A multidimensional Markov jump process describes the epidemic dynamics in a closed population of size N. Using the large population framework, i.e., with N large, we first build an approximation of epidemic dynamics using an autoregressive Gaussian process via a diffusion approach (see e.g., [START_REF] Ethier | Markov processes: characterization and convergence[END_REF], [START_REF] Guy | Approximation of epidemic models by diffusion processes and their statistical inference[END_REF]). Then we simultaneously account for a given missing coordinate value and systematic noise present in observations by applying a projection operator to the process and adding heteroscedastic Gaussian errors. This yields the theoretical framework that allows recursive computations of an approximate likelihood. This approach, based on Kalman filtering, enables the computation of the approximate log-likelihood of the available observations and, consequently, the estimation of model parameters. An initial innovative aspect of this method with respect to others is the use of a Kalman filter to recursively compute the approximate likelihood in the non-standard case of the small noise framework (i.e., with noise covariance matrix proportional to 1/N), rather than the classical recurrent case coupled with a large observation time-window (with the number of observations going to infinity). In addition, the explicit integration into the algorithm of the data sampling interval, and an alternative point of view in the prediction of successive model states-given the observations-are further innovative points.

The derivation and accuracy assessment of Gaussian process approximation for stochastic epidemic models have previously been described in [START_REF] Buckingham-Jeffery | Gaussian process approximations for fast inference from infectious disease data[END_REF], along with maximum likelihood inference for parameters underlying epidemic dynamics. However, that study does not rely on Kalman filtering, nor does it consider noise in outbreak data. Computation of the approximate likelihood of the associated statistical model, as well as parameter estimation, performed via Kalman filtering recursion was proposed in [START_REF] Favetto | Parameter estimation for a bidimensional partially observed ornstein-uhlenbeck process with biological application[END_REF], but for simpler models without nonlinear terms in the drift, and with no parameter to estimate in the diffusion term.

For the sake of simplicity, we consider here an epidemic with homogeneous mixing in a closed population whose dynamics are described by a compartmental model, with each compartment containing individuals with identical health states. We focus on the simple SIR (susceptible -infectious -recovered) epidemic model characterized by a two-dimensional jump process, partially observed at regularly-spaced discrete times, with measurement errors. The approach can be easily extended to broader epidemic models observed with various sampling intervals. The paper is organized as follows. In Section 2 we introduce the general framework and related inference issues, and propose the model approximation. Section 3 contains the main methodological developments of our paper: construction of the approximate log-likelihood, its computation based on Kalman filtering recursion, and the associated parameter estimation. In Sections 4 and 5 we assess the performance of our estimators on both simulated data and real data from an influenza outbreak in a British boarding school in 1978, and compare our results with those obtained using the MIF method. Section 6 contains a discussion and concluding remarks.

Gaussian model approximation for large population epidemics

Preliminary comments on inference in epidemic models

Epidemic dynamics can be naturally described using compartmental models, which are by essence mechanistic and include parameters in their characterization. In such models, the population is partitioned into compartments corresponding to different stages of the infection process, whose temporal evolution is described. As an illustrative example throughout the article, we will use the simple SIR epidemic model. At any time, each individual is either susceptible (S), infectious (I), or recovered (R). In this model, there are two mechanistic parameters of interest that govern the transitions of individuals between states S, I, and R: the transmission rate of the pathogen λ and the recovery rate γ. More precisely, individuals can move from state S to I according to λ, or from state I to R according to γ (Figure 1). One of the main goals of epidemic studies is to estimate such mechanistic parameters from the available data. One of the most natural probabilistic representations of compartmental epidemic models is the continuous-time Markov jump process (see Section 2.2). Inference for Markov jump processes is straightforward when sample paths are completely observed. In the context of epidemics, this is equivalent to the observation of all infection and recovery times for all individuals in the population. This rarely occurs in practice; often one or more of the coordinates (i.e., S (t), I(t)) are not observed, and available observations are only collected at discrete time points

S I R

λI/N γ

t k with 0 = t 0 < t 1 < t 2 < • • • < t n = T over a finite time interval [0, T ].
More specifically, the data often consists of counting newly infected individuals N I (t k ) on successive time intervals

[t k-1 , t k ].
Alternatively, the successive numbers of infectious individuals I(t k ) are sometimes available, especially for low population sizes. Moreover, it is common that the available data is affected by several sources of noise such as under-reporting of infection events or-when reported-imperfect diagnostic tests. Essentially, the nature of such data makes it difficult to infer key epidemic parameters: (i) observations are available at discrete time points, (ii) not all coordinates of the dynamical model are observed, and (iii) systematic reporting and measurement errors have to be taken into account.

Approximation of large population epidemic models and the autoregressive point of view

Consider an epidemic in a closed population with homogeneous mixing modeled by a ddimensional Markov jump process Z(t), where d is the number of compartments corresponding to successive health statuses within the population. If N is the population size, the state space of (Z(t), t ≥ 0) is E = {0, . . . , N} d . Let Q = (q k,l , k, l ∈ E) denote its Q-matrix; the latter satisfies ∀l k, q k,l ≥ 0, and q k,k = -l∈E,l k q k,l . There are two standard ways of describing this jump process (see e.g., [START_REF] Norris | Markov chains[END_REF]):

(i) By the underlying jump chain and holding times. Starting from Q, set π k,l = q k,l q k with q k = -q k,k , π k,k = 0 if q k 0, and π k,k = 1 if q k = 0. The process stays in state k according to an exponential distribution E(q k ) and jumps to state l with probability π k,l .

(ii) Using its infinitesimal generator: as h → 0, P(Z(t

+ h) = l|Z(t) = k) = δ k,l + q k,l h + o P (h),
where δ k,l denotes the Kronecker function (

δ k,l = 1 if l = k, δ k,l = 0 if l k).
Hence, for f a measurable function E → R, if E k denotes the expectation conditional on

Z(0) = k, [Q f ](k) = l∈E q k,l f (l) = lim t→0 1 t (E k f (Z(t)) -f (k))
. Simulations of Z(t) are usually based on (i), while (ii) relies on general properties of Markov processes.

For any vector V or matrix M, let V t or M t denote their transpose. For a jump (0, . . . , 0) t of Z(t), we define the jump function:

α (k) = q k,k+ for k, k + ∈ E.
Consider now the normalized Markov jump process (Z N (t)) t≥0 :

Z N (t) = Z(t) N ∈ E N = {k/N, k ∈ E}. (1)
The associated jump functions are, for

x ∈ E N , α N (x) = 1 N α ([N x]). Assume that the process (Z(t)) is density-dependent, i.e., H1 : ∀ , ∀x ∈ [0, 1] d , 1 N α ([N x]) → N→+∞ β (x), H2 : ∀ , β ∈ C 2 ([0, 1] d , R), where [N x] is the vector of integers [N x 1 ], . . . , [N x d ], with [N x i ] the integer part of N x i . Next, define for x ∈ [0, 1] d the function b(•) and the d × d symmetric non-negative matrix Σ(•): b(x) = ∈E - β (x) ; Σ(x) = ∈E - β (x) t . (2)
For the SIR epidemic model in a closed population, we have that S (t) + I(t) + R(t) = N for all t. Therefore, its state space is E = {0, . . . , N} 2 . Only two jumps are possible from k = (S , I) t :

• 1 = (-1, +1) t : (S , I) → (S -1, I + 1) ⇒ q k,k+ 1 = λS I/N = α 1 (k), • 2 = (0, -1) t : (S , I) → (S , I -1) ⇒ q k,k+ 2 = γI = α 2 (k).
This process is density dependent

: if s = S N , i = I N , then 1 N α 1 ([N s], [Ni]) = 1 N (λ[N s]) [Ni] N → λsi and 1 N α 2 ([N s], [Ni]) = 1 N γ[Ni] → γi as N → ∞. Moreover (2) is, for x = s i , b(x) = λsi -1 1 + γi 0 -1 , Σ(x) = λsi 1 -1 -1 1 + γi 0 0 0 1 .
We now recall the law of large numbers result stated (for instance) in [START_REF] Britton | Stochastic epidemic models with inference[END_REF].

Lemma 1. Assume that (Z(t)) satisfies (H1), (H2), and Z N (0) → x 0 as N → +∞. Then, (Z N (t)) converges almost surely uniformly on [0, T ] to the solution x(t) of the ordinary differential equa-

tion dx dt = b(x(t)); x(0) = x 0 . (3)
If x 0 = (0, . . . , 0) t , then x(t) = 0 for all t and (3) no longer adequately describes the epidemic dynamics (see e.g., [START_REF] Britton | Stochastic epidemic models with inference[END_REF] Part I). Equation (3) describes the dynamics in the case of a major outbreak corresponding to x 0 (0, . . . , 0) t . In [START_REF] Guy | Approximation of epidemic models by diffusion processes and their statistical inference[END_REF], by extending the results of [START_REF] Ethier | Markov processes: characterization and convergence[END_REF], another approximation of the epidemic model was proposed, leading to a diffusion process (Z N (t)) t≥0 with the small diffusion matrix 1 N Σ(x), where Σ is the matrix defined in (2):

       dZ N (t) = b(Z N (t)) + 1 √ N σ(Z N (t)) dB(t), Z N (0) = x 0 , (4)
where (B(t)) t≥0 is a d-dimensional Brownian motion and σ a d × d matrix such that

σ(x)σ t (x) = Σ(x). (5)
For stochastic differential equations with small noise (i.e., proportional to 1 N ), an approximation of Z N (t) can be obtained using (2)-( 5), based on the theory of perturbations of dynamical systems (see e.g., [START_REF] Azencott | Formule de taylor stochastique et développement asymptotique intégrales de feynmann[END_REF], [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF]):

             Z N (t) = x(t) + 1 √ N g(t) + 1 √ N R N (t), dg(t) = ∇ x b(x(t)) g(t) dt + σ(x(t)) dB(t) ; g(0) = 0, with sup t R N (t) → 0 in probability as N → +∞, (6)
where ∇ x b(x) denotes the matrix ( ∂b i ∂x j (x)) 1≤i, j≤d . The stochastic differential equation for g(•) defined in (6) can be solved explicitly (see e.g., [START_REF] Guy | Parametric inference for discretely observed multidimensional diffusions with small diffusion coefficient[END_REF] for details) and its solution is the time-inhomogeneous Gaussian process

g(t) = t 0 Φ(t, s)σ(x(s)) dB(s), (7)
where

Φ(t, s) satisfies ∂Φ ∂t (t, s) = ∇ x b(x(t))Φ(t, s), Φ(s, s) = I d . Hence, Φ(t, s) is the d × d matrix Φ(t, s) = exp t s ∇ x b(x(u)) du . (8)
Using ( 3) and ( 7), let us define the Gaussian process G N (t):

G N (t) = x(t) + 1 √ N g(t). (9)
Consider now the Wasserstein-1 distance on the interval [0, T ] between R d -valued processes

U t , V t on [0, T ]. W 1,T (U, V) = inf E(||U -V|| T ), where if x : [0, T ] → R d , ||x|| T = sup 0≤t≤T ||x(t)||,
and the above infimum is over all couplings of two processes. According to [START_REF] Britton | Stochastic epidemic models with inference[END_REF], Part I, Theorem 2.4.1, the following holds.

Proposition 1. For all T > 0, the Wasserstein-1 distances on [0, T ] between the three processes

(Z N (•)), (Z N (•)), and (G N (•)) defined in (1),(4),(9) satisfy, as N → ∞, √ NW 1,T (Z N , Z N ) → 0, √ NW 1,T (Z N , G N ) → 0, and
√ NW 1,T (Z N , G N ) → 0.
From a statistical point of view, this proposition has important consequences: given the fact that these distances are o(N -1/2 ), we develop our inference method by plugging the observations into the likelihood of either the diffusion process (Z N ) or the Gaussian process (G N ). This approach is often used to derive approximate likelihoods or contrasts for stochastic processes. For instance, for discretely observed diffusion processes, parametric inference is often based on the likelihood of the Euler scheme of the diffusion (see e.g., [START_REF] Kessler | Statistical Methods for Stochastic Differential Equations[END_REF]). Moreover, it was proved in [START_REF] Guy | Parametric inference for discretely observed multidimensional diffusions with small diffusion coefficient[END_REF] that parametric inference based on (G N ) leads to efficient estimators for the parameters ruling the jump process.

From here on, we will use the approximation of (Z N ) by the Gaussian process (G N ). Let us now consider a parametric model for epidemic dynamics. This yields a parametric continuoustime approximate model for epidemic dynamics, with parameter

η = (ζ, x 0 ), ( 10
)
where ζ contains the parameters found in the transition rates of the jump process, and therefore in the functions β (x) defined in H1, and x 0 is the initial point of the ordinary differential equation (ODE) defined in Lemma 1. As mentioned in Section 2.1, the process is however observed at discrete times t k , where (t k ) is an increasing sequence on [0, T ], with

t 0 = 0 < t 1 • • • < t n = T .
We therefore deduce from above a discrete-time representation of the epidemic evolution.

Let us denote by F t = σ(B(s), s ≤ t). Then the following holds.

Proposition 2. There exists a sequence of independent Gaussian random variables (U k ) such that (i) For all k, U k is F t k -measurable and independent of F t k-1 .

(ii) The process G N defined in (9) is an AR(1) process and satisfies, using (3), ( 8), G N (0) = x 0 , for k ≥ 1,

G N (t k ) = F k (η) + A k-1 (η) G N (t k-1 ) + U k , where A k-1 (η) = A(η, t k-1 ) = Φ(η, t k , t k-1 ), F k (η) = F(η, t k ) = x(η, t k ) -Φ(η, t k , t k-1 )x(η, t k-1 ),
and (U k ) are independent random variables such that

U k ∼ N d (0, T k (η)), with T k (η) = 1 N t k t k-1 Φ(η, t k , s)Σ(η, x(η, s)) Φ t (η, t k , s)ds.
The proof of Proposition 2 is given in the Appendix. Using now that sup

t ||Z N (t) -G N (t)|| = 1 √ N o P (1), Proposition 2 becomes, setting X k := X(t k ) = Z N (t k ), X 0 = x 0 , for k ≥ 1, X k = F k (η) + A k-1 (η)X k-1 + U k . (11)

Approximation of the observation model

Assume now that there are noisy observations O(t k ) of the original jump process Z(t) (with state space E = {0, . . . , N} d at discrete times t k ). As mentioned in Section 2.1, it often occurs in practice that not all epidemiological health states are observed. We account for this by introducing a projection operator B : R d → R q with q ≤ d, where BX(•) contains only the coordinates that can be observed. Therefore B is a d × q matrix whose elements are 0 and 1. For k = 0, . . . , n, define C(t k ) = (C 1 (t k ), . . . , C q (t k )) t = BZ(t k ) ∈ {0, . . . , N} q .

In an initial approach, assume that each component of C is observed with independent reporting rate p i and measurement errors. In this way, we propose a rather general model for the observations conditional on Z(t), for 1

≤ i ≤ q: O i (t k ) = O i,1 (t k ) + O i,2 (t k ), with O i,1 (t k ) ∼ Binomial(C i (t k ), p i ), O i,2 (t k ) ∼ N 0, τ 2 i C i (t k ) , (12)
where, conditional on σ(Z(s), 0 ≤ s ≤ t k ), the variables O i,1 (t k ) and O i,2 (t k ) are independent. This yields a new higher-dimensional parameter containing parameters for both the epidemic (i.e., η defined in (10)) and observation processes: θ = (η, (p 1 , . . . , p q ), (τ 2 1 , . . . , τ 2 q )).

Consider now the normalized process Z N (t). We can then define C N (t) = BZ N (t) and associated normalized observations

O N (t k ) = 1 N O(t k ).
A Gaussian approximation of the observation process has first and second moments which satisfy

E(O N,i (t k )|Z(t k )) = p i C N,i (t k ), Var(O N,i (t k )|Z(t k )) = 1 N (p i (1 -p i ) + τ 2 i )C N,i (t k ).
Using now (6) and Proposition 1, we get that

C N (t) = BZ N (t) = Bx(η, t) + 1 √ N Bg(η, t) + 1 √ N o P (1).
The Gaussian process g(η, t) is uniformly bounded in probability on [0, T ], so we have that

Var(O N,i (t k )|Z(t k )) = 1 N (p i (1 -p i ) + τ 2 i )(Bx(η, t k )) i + O P (N -3/2 ).
Let us next define the q-dimensional matrices

P(θ) = diag(p i ) 1≤i≤q , Q k (θ) = 1 N diag (p i (1 -p i ) + τ 2 i )(Bx(η, t k )) i , (13)
and the q × d matrix B(θ) = P(θ)B.

The Gaussian approximations (Y k ) of the observations O N (t k ) satisfy that conditionally on Z(t k ),

Y k = B(θ)X k + V k with V k ∼ N q (0, Q k (θ)), (14)
where (V k ) are independent random variables such that for all k, V k is independent of Z N (t k ).

Application on the SIR epidemic model

Let us now illustrate the model approximations derived in Sections 2.2 and 2.3 on the simple SIR model introduced in Section 2.1. The Markov jump process Z(t) = (S (t), I(t))), t ≥ 0 is defined in Section 2.2. The parameters controlling the dynamics of the system are η = (λ, γ, x 0 ) = (λ, γ, s 0 , i 0 ), which include the transition rates λ and γ, and the initial point x 0 = (s 0 , i 0 ) (cf Lemma 1).

Dynamical state model. Let us define the key quantities necessary to derive the appropriate Gaussian process (G N (t)) as defined in ( 9), including the dependence on η:

G N (t) = x(η, t) + 1 √ N g(η, t).
The first important element is x(η, t) = (s(η, t), i(η, t)) t , solution of the following ODEs:

             ds dt (η, t) = -λs(η, t)i(η, t), di dt (η, t) = λs(η, t)i(η, t) -γi(η, t), x 0 = (s 0 , i 0 ).
When there is no ambiguity, we denote by s and i respectively s(η, t) and i(η, t). Then, to get g(η, •), we need to derive the functions b(η, •) and

Σ(η, •) from (2) (see Section 2.2): b(η, s, i) = -λsi λsi -γi ; Σ(η, s, i) = λsi -λsi -λsi λsi + γi , (15)
and the Cholesky decomposition of Σ(η, •):

σ(η, s, i) = √ λsi 0 - √ λsi √ γi .
From (15), we deduce the gradient of b:

∇ x b(η, s, i) = -λi -λs λi λs -γ ,
and the resolvent matrix defined in (8):

Φ(η, t, s) = exp t s ∇ x b(η, x(η, u))du .
Finally, we obtain

g(η, t) = t 0 Φ(η, t, u)σ(η, x(η, u))dB(u),
where (B(u)) u≥0 is a bidimensionnal Brownian motion.

Discrete-time system. For simplicity, we assume a regular sampling:

t k = k∆, k = 0, . . . , n, T = n∆.
The dependence with respect to ∆ is explicitly given in the equations. The approximate autoregressive model, setting

X k = Z N (t k ) = (S N (k∆), I N (k∆)) t
, is given by:

             X k = F k (η, ∆) + A k-1 (η, ∆)X k-1 + U k , where F k (η, ∆) = x(η, t k ) -Φ(η, t k , t k-1 )x(η, t k-1 ), A k-1 (η, ∆) = Φ(η, t k , t k-1 ), U k ∼ N 2 (0, T k (η, ∆)) with T k (η, ∆) = 1 N t k t k-1 Φ(η, t k , s)Σ(η, x(η, s))Φ t (η, t k , s)ds.
(16) Observation model. Suppose for example that only the infected individuals are observed with reporting and measurement errors. This corresponds to considering in (12):

O 1 (t k ) ∼ Binomial(I(t k ), p), O 2 (t k ) ∼ N(0, τ 2 I(t k )). ( 17
)
Hence the full parameter vector is θ = (λ, γ, s 0 , i 0 , p, τ 2 ). To derive (14) from this example, we define the operator B(θ) = pB, where B : (x 1 , x 2 ) t → x 2 is the projection operator on the infected compartment, and

Q k (θ) = 1 N (p(1 -p) + τ 2 )i(η, t k ), with Q k is defined in (13)
. By joining ( 16) with the Gaussian approximate observation model defined above, we get the following discrete-time state-space model:

       X k = F k (η, ∆) + A k-1 (η, ∆)X k-1 + U k , with U k ∼ N 2 (0, T k (η, ∆)) , Y k = p 0 1 X k + V k , with V k ∼ N 0, 1 N (p(1 -p) + τ 2 )i(η, t k ) .

Parameter estimation using Kalman filtering techniques

Approximate likelihood inference

The parameters of interest in the general case are denoted by θ = (η, (p 1 , . . . , p q ), (τ 2 1 , . . . , τ 2 q )), where η contains the parameters controlling the dynamics and x 0 , whereas (p 1 , . . . , p q ) and (τ 2 1 , . . . , τ 2 q ) are derived from the reporting and measurements errors in the observations. Our aim is to estimate the unknown parameters θ from observations y n:0 = (y 0 , . . . , y n ) obtained at discrete time points t 0 < t 1 < • • • < t n . Joining ( 11) and ( 14), we get the following discrete-time Gaussian state-space setting that is more convenient for inference:

       X k = F k (η) + A k-1 (η)X k-1 + U k , Y k = B(θ)X k + V k , (18)
where all quantities are explicitly defined in Sections 2.2 and 2.3. Using ( 18), we propose to estimate θ by maximizing the associated likelihood

L(•; Y 0 , . . . , Y n ): θ = argmax θ L(θ; Y 0 , . . . , Y n ). ( 19
)
The log-likelihood of the observations y 0 , . . . , y n is given by:

L(θ; y 0 , . . . , y n ) = log f (θ, y 0 ) + n k=1 log f k (θ; y k |y k-1:0 ). ( 20
)
Computing L(θ; y 0 , . . . , y n ) requires the computation of the two first moments of the Gaussian conditional distributions corresponding to each log f k (θ; . . .) term. This relies on the computation of the predictive distributions

ν k|k-1:0 (θ; dx) = L(X k |y k-1:0 ), k ≥ 1, from which we derive the conditional densities f k (θ; y k |y k-1:0 ) = f (y k |x)ν k|k-1:0 (θ; dx).
Usually, these conditional distributions are obtained by means of filtering methods, based on the iterative computations of the conditional distributions:

• the predictive distribution: L(X k |y k-1 , . . . , y 0 ) = ν k|k-1:0 (dx), k ≥ 1, with the convention ν 0,0 (dx) = L(X 0 ),

• the updating distribution: L(X k |y k , . . . , y 0 ) = ν k|k:0 (dx), k ≥ 0,

• the marginal distribution: L(Y k |y k-1 , . . . , y 0 ) = µ k|k-1:0 , k ≥ 1, with the convention µ 0|0:0 (dx) = L(Y 0 ).
In the special case of the Gaussian state space model and Gaussian noise, all of these distributions are Gaussian and therefore characterized by their mean and covariance matrix. Using notation specific to Kalman filtering, let us set

L(X k |y k-1 , • • • , y 0 ) = ν k|k-1:0 (dx) = N d ( Xk , Ξk ) (predictive distribution). L(X k |y k , • • • , y 0 ) = ν k|k:0 (dx) = N d ( Xk , Tk ) (updating distribution). L(Y k |y k-1 , • • • , y 0 ) = µ k|k-1:0 = N q ( Mk , Ωk ) (marginal distribution).
The Gaussian approximations defined in ( 11), ( 14) and ( 18) allow us to use specific properties of Gaussian distributions that are recalled below.

3.1.1. Preliminary results in the general framework of Kalman filtering Let (X i , i ≥ 0) be a non-centered d-dimensional Gaussian AR(1) process and assume that only q coordinates of (X i ) are observed, with Gaussian noise. Computations of the conditional distributions rely on a Kalman filter approach, which is derived from the following lemma.

Lemma 2. Assume that X is a random variable with distribution N d (ξ, T ) which conditional on X, Y has distribution N q (BX, Q). Then, L(X|Y) is Gaussian:

N d ( ξ(y), T ), with ξ(y) = ξ + T B t (BT B t + Q) -1 (y -Bξ); T = T -T B t (BT B t + Q) -1 BT. ( 21
)
Remark 1. We stress that Lemma 2 holds even if Q is singular. In particular, the formula holds when Q = 0 and B is a projection operator, i.e., the observations are Y k = BX k , provided that T is non-singular.

Let us now go back to our general setting (X k , Y k ) defined in (18).

Proposition 3. Assume that (X k , Y k ) are defined as in (18). Then, ν k|k-1:0 (dx), ν k|k:0 (dx), and µ k|k-1:0 )(dy) satisfy, with the initialization X0 = ξ 0 , Ξ0 = T 0 , for k ≥ 0,

(i) Prediction: ν k|k-1:0 (dx) ∼ N d ( Xk , Ξk ) with Xk = F k + A k-1 Xk-1 , Ξk = A k-1 Tk-1 A t k-1 + T k . (ii) Updating: ν k|k:0 (dx) ∼ N d ( Xk , Tk ) with Xk = Xk + Ξk B t (B Ξk B t + Q k ) -1 (Y k -B Xk ), Tk = Ξk -Ξk B t (B Ξk B t + Q k ) -1 B Ξk . (iii) Marginal distribution: µ k+1|k:0 (dy) ∼ N q ( Mk+1 , Ωk+1 ) with Mk+1 = B Xk+1 , Ωk+1 = B Ξk+1 B t + Q k+1 .
Using specific notation from Kalman filtering, we recover a modified version of the Kalman algorithm. Assume that X 0 ∼ N d (ξ 0 , T 0 ) and that, for all k ≥ 0, the matrices Γ k defined below are non-singular. Then, setting X0 = ξ 0 , Ξ0 = T 0 , we have

k-1 = Y k-1 -B Xk-1 , (innovation) Γ k-1 = B Ξk-1 B t + Q k-1 , (innovation covariance) H k-1 = A k-1 Ξk-1 B t Γ -1 k-1 , (Kalman gain) Xk = F k + A k-1 Xk-1 + H k-1 k-1 , (predicted mean state estimation) Ξk = (A k-1 -H k-1 B) Ξk-1 A t k-1 + T k . (predicted error covariance)
Therefore, the marginal distributions appearing in the computation of the log-likelihood (20) are µ k+1|k:0 (dy) ∼ N q ( Mk+1 , Ωk+1 ), with

Mk+1 = B Xk+1 , Ωk+1 = B Ξk+1 B t + Q k+1 . (22)
3.1.2. Recursive computation of the approximate log-likelihood An important consequence of the previous section is that we can compute (20) based on the recursive computations of the first moments of the Gaussian distributions corresponding to each term of the log-likelihood. By explicitly accounting for the dependence on θ of moments given in ( 22), we obtain:

L(θ; y 0 , . . . , y n ) = C + log f (θ; y 0 ) - 1 2 n k=1 log | Ωk (θ)| + (y i -Mk (θ)) t Ωk (θ) -1 (y i -Mk (θ)) ,
with C a constant (independent of the parameters) and |A| denoting the determinant of the matrix A.

Note that the sampling interval ∆ plays an important role in the various key quantities involved in the Kalman recursions (see Appendix A for details).

Application on the SIR epidemic model

Let us again take the example of SIR epidemics, when only the infected individuals are observed with reporting and measurement errors, considered in Section 2.4. By assuming an initial distribution X 0 ∼ N 2 (ξ 0 , T 0 ), setting X0 = ξ 0 , Ξ0 = T 0 , and applying the algorithm given in Proposition 3, we have, for k = 0, . . . , n -1:

k-1 (θ) = Y k-1 -p Îk-1 (θ), (scalar) Γ k-1 (θ) = p 2 ( Ξk-1 (θ)) 22 + 1 N (p(1 -p) + τ 2 )i(η, t k-1 ), (scalar)
H k-1 (θ) = pA k-1 (η) Ξk-1 (θ) 0 1 Γ -1 k-1 (θ), (vector) Xk (θ) = F k (η) + A k-1 (η) Xk-1 (θ) + H k-1 (θ) k-1 (θ), (vector) Ξk (θ) = A k-1 (η) -pH k-1 (θ) 0 1 Ξk-1 (θ)A k-1 (η) t + T k (η). (2 × 2 matrix)
This yields the marginal distributions:

Mk+1 (θ) = p Îk+1 (θ), Ωk+1 (θ) = p 2 Ξk+1 (θ) 22 + 1 N (p(1 -p) + τ 2 )i(η, t k+1 )),
which are used to compute the likelihood

L(θ, y 1 , . . . , y n ) - 1 2 n k=1 log Ωk (θ) - 1 2 n k=1 (y k -Mk (θ)) 2
Ωk (θ) .

Simulation study

We assessed the performance of our method on simulated SIR epidemics in which only the infectious compartment is observed at discrete time points (see Section 2.4 where the model is fully described).

Simulation settings

Data simulation. We first simulated SIR dynamics according to the Markov jump process using the Gillespie algorithm [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]). Only trajectories that did not exhibit early extinction were considered for inference. The theoretical proportion of these trajectories is given by 1 -(γ/λ) I 0 [START_REF] Andersson | Stochastic epidemic models and their statistical analysis[END_REF]), where I 0 is the number of infectious individuals at time 0. We simulated two cases. First, for the emergent trajectories, the observations were generated by binomial draws from I(t) at n + 1 discrete time points t 0 < t 1 < . . . < t n . In (17), this amounts to considering τ = 0, with simulated observations finally obtained via

O(t k ) = O 1 (t k ) ∼ Binomial(I(t k ), p)
. Second, we considered the more general case where observations are

O(t k ) = O 1 (t k ) + O 2 (t k ), with O 1 (t k ) ∼ Binomial(I(t k ), p), O 2 (t k ) ∼ N 0, τ 2 I(t k ) ,
where the non-zero measurement error τ is an additional parameter to estimate. Figure 2 represents epidemic trajectories corresponding to the various steps of data simulation. These plots illustrate the variability in the stochastic trajectories compared to the deterministic counterpart of the SIR model, and the loss of information from the unobservable real dynamics to the observations available for inference. Moreover, the second source of error, driven by the measurement error τ, seems to have a minor impact on the global observational noise compared to the reporting error. The evolution of the number of susceptible individuals is not shown in Figure 2. From the point of view of inference, the S compartment is a latent variable, the observations being only available for the infected state. Numerical scenarios. We used the following parameter values for the simulation of the epidemics: λ = 1, γ = 1/3, and initial starting points S (0)/N = s 0 = 0.99, I(0)/N = i 0 = 0.01, R(0) = 0 (hence with s 0 + i 0 = 1). Observations were generated under two scenarios: i) high reporting rate p = 0.8 and ii) low reporting rate p = 0.3. Two experiments were considered concerning the measurement error: τ = 0 (experiment 1) and τ = 0.5 (experiment 2). Scenarios combining three population sizes (N ∈ {1000, 2000, 10000}) with different values for the number of observations (n) for each epidemic trajectory were also investigated. For each value of N, conditionally on non extinction, 500 SIR epidemic dynamics were simulated. Observations were generated at regularly-spaced time points t k = k∆ using, for a given scenario, the same value of ∆ for each of the 500 epidemics (obtained by dividing the mean epidemic duration over 500 trajectories by a target number of observations n). As the epidemic duration is stochastic, we considered slightly different observation intervals [0, T ] for each epidemic and set the value of T as the first time point when the number of infected individuals became zero. This generates slightly different numbers of observations per epidemic trajectory.

Inference: settings, performance comparison, and implementation

The unknown parameters to be estimated are either θ = (λ, γ, p, i 0 ) or θ = (λ, γ, p, i 0 , τ), according to the experiment. Here, we do not need to estimate s 0 as s 0 = 1i 0 . When τ 0, the observational model used for the two estimation methods was a Gaussian model given as the sum of the two sources of noise in the data (reporting: Gaussian approximation of a binomial model; measurement: Gaussian model). For each simulated dataset, θ is estimated with our Kalman filter-based estimation method (KM) and with the MIF algorithm [START_REF] Ionides | Inference for nonlinear dynamical systems[END_REF], [START_REF] Ionides | Iterated filtering[END_REF], [START_REF] Ionides | Inference for dynamic and latent variable models via iterated, perturbed bayes maps[END_REF]), which is widely used in practice for statistical inference of epidemics. The simulation study was performed with the R software on a Bi-pro Xeon E5-2680 processor with 2.8 Ghz, 96 Go RAM, and 20 cores. MIF estimation was performed with the mif2 function of the POMP-package [START_REF] King | Statistical inference for partially observed markov processes via the r package pomp[END_REF]). We provide user-friendly code on the RunMyCode website (see Appendix F for details).

Let us make some initial remarks on the algorithms and their practical implementations. Regardless of the estimation method used, maximisation of the log-likelihood requires considering several constraints: (i) strict positivity of λ, γ, i 0 , (ii) s 0 +i 0 = 1 (or s 0 +i 0 ≤ 1 in the general case), and (iii) 0 < p ≤ 1. To facilitate optimization, a different parameterization was implemented:

λ = exp(µ 1 ), γ = exp(µ 2 ), p = (1 + exp(µ 3 )) -1 , i 0 = (1 + exp(µ 4 )) -1 , where µ 1 , µ 2 , µ 3 , µ 4 ∈ R.
With no constraints on this new set of parameters, numerical optimization was more stable in practice. The approximated log-likelihood given by Kalman filtering techniques cannot be maximized explicitly. We instead used the Nelder-Mead method implemented in the optim function in R, which requires inputting initial values for the unknown parameters. According to the amount of information available in the observations, the result of the optimization is more or less sensitive to these initial points. The same problem can occur for the MIF algorithm. The dependence on the initialization can be circumvented by trying different starting values (10 in the present case) and choosing the maximum value for the log-likelihood among them. The starting parameter values for the maximization algorithm were uniformly drawn from a hypercube encompassing the likely true values.

When the time intervals ∆ between observations are large (which often occurs for low values of n), we computed the resolvent matrix defined in (8) as in (A.1) in order to obtain the approximated log-likelihood with Kalman filtering techniques. MIF, based on particle filtering, returns an estimate of the log-likelihood of the observations by using resampling techniques. The parameter space is investigated by randomly perturbing the parameters of interest at each iteration, the amplitude of the perturbation decreasing as the iterations progress. The MIF algorithm has a complexity of O(JM), where J and M are respectively the number of particles and the number of iterations. Running MIF requires specifying several tuning parameters. For the present study, the best results were obtained using M = 100 iterations, J = 500 particles, standard deviation rw.sd equal to 0.2 for the random walk for each parameter, and a cooling of the perturbations of cooling.fraction.50=0.05 in the POMP-package (we drew inspiration from [START_REF] Stocks | Iterated filtering methods for markov process epidemic models[END_REF] for this choice of tuning parameters).

Concerning implementation issues, in our experience, the tuning of the MIF algorithm (number of particles, number of iterations, etc.) can greatly affect the quality of the estimates. In particular, it seems that there is an important interplay between the tuning parameters and the initialization values of model parameters to be inferred. In comparison, our method has only one main calibration parameter in practice. In the filtering step, it is necessary to initialize the covariance matrix (i.e., T 0 in Section 3.2) of the state variables, conditional on the observations, but it seems that this initialization does not have a noticeable influence on the accuracy of estimates. 1 and2 respectively display the results for the high reporting scenario (p = 0.8) and the low reporting scenario (p = 0.3). Each table compares estimates obtained with KM and MIF. For each parameter and each estimation method, the reported values are the mean of the 500 parameter estimates, with their standard deviations in brackets.

These results show that, irrespective of the reporting rate p, when the population size N and the number of observations n per epidemic increase, the bias and the standard error of the estimates obtained decrease, whichever method is used for inference. For a given (N, n), the estimation bias is higher when the reporting rate is low (p * = 0.3, where the star here designates the true value). This may be partly related to the fact that the information contained in the data decreases as p * decreases. Both methods provide estimates with comparable levels of accuracy.

Table 1: First experiment (τ = 0). Estimation of θ = (λ, γ, p, i 0 ) under the constraint s 0 + i 0 = 1 in Setting 1 with true parameter values (λ * , γ * , p * , i * 0 )=(1, 1/3, 0.8, 0.01). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by our Kalman-based method (KM) and Maximum Iterated Filtering (MIF). The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. Table 2: First experiment (τ = 0). Estimation of θ = (λ, γ, p, i 0 ) under the constraint s 0 + i 0 = 1 in Setting 2 with true parameter values (λ * , γ * , p * , i * 0 )=(1, 1/3, 0.3, 0.01). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by KM and MIF. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. The estimates are less computationally demanding and require less algorithmic tuning with the Kalman filtering approach. This simulation study was also performed for a second set of parameter values (λ = 0.6, γ = 0.4, i 0 = 0.01), under the constraint s 0 + i 0 = 1 and for p = 0.8 and p = 0.3, and naturally led to greater variability between simulated trajectories. These results are provided in Appendix E for comparative purposes.

Simulation results for the second experiment (τ 0)

Here, we present the estimation results when the simulated observations are obtained with a non-zero measurement error τ, which is to be estimated. As noticed in [START_REF] Stocks | Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in germany[END_REF], the initial conditions of the system are difficult to estimate, and usually set at plausible values. Consequently, we distinguish two situations, where either (i) i 0 is unknown and estimated; or (ii) i 0 is known and fixed.

Unknown starting point i 0 . Five different target values for sample sizes were considered: n = 10, n = 30, n = 100, n = 500, and n = 1000. The unknown parameters to be estimated were θ = (λ, γ, p, i 0 , τ) under the constraint s 0 + i 0 = 1. For the sake of clarity, we do not show the results when N = 2000 and p = 0.3. Results are displayed in Table 3.

Table 3: Second experiment (τ 0). Estimation of θ = (λ, γ, p, i 0 , τ) under the constraint s 0 + i 0 = 1 in Setting 1 with true parameter values (λ * , γ * , p * , i * 0 , τ * )=(1, 1/3, 0.8, 0.01, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by our Kalman-based method and the MIF algorithm. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. As in the first experiment where τ = 0, the results show that the estimations provided by KM and MIF are of the same order of accuracy. The pattern concerning the bias and the standard error observed in the case τ = 0 also occurs when τ = 0.5 is estimated, i.e., bias decreasing and accuracy increasing when N and n increase. We remark that the estimation is more difficult, inducing larger bias, when the measurement error τ is unknown, even for a quite large number of observations n ≈ 100. Consider for example N = 1000, n = 101 and p * = 0.8. The point estimate value of p obtained by KM with τ = 0 (cf. Table 1) and τ 0 (cf. Table 3) is respectively 0.79 and 0.74. This is more marked for the second set of parameters values (λ = 0.6 and γ = 0.4), presented in Appendix E, which induces more variability between epidemics. For N = 1000, n = 99 and p * = 0.8, comparing the results in Tables E.6 and E.8 shows that p passes from 0.75 to 0.66 when τ = 0.5 unknown. Higher frequency observations of the epidemics lead to more satisfactory estimations: considering n = 998 when τ = 0.5 unknown leads to p = 0.78. The estimates obtained with MIF behave similarly. In summary, when the measurement error τ is non-zero and estimated, a greater number of observations is needed in order to obtain estimates without bias for both the Kalman-based and MIF methods.

Known starting point i 0 . The unknown parameters to be estimated are θ = (λ, γ, p, τ). Tables 4 and5 respectively display the results for the high reporting scenario (p = 0.8) and low reporting scenario (p = 0.3).

Table 4: Second experiment (τ 0). Estimation of θ = (λ, γ, p, τ) with s 0 = 0.99 and i 0 = 0.01 known in Setting 1 with true parameter values (λ * , γ * , p * , τ * )=(1, 1/3, 0.8, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by KM and MIF. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. Table 5: Second experiment (τ 0). Estimation of θ = (λ, γ, p, τ) with s 0 = 0.99 and i 0 = 0.01 known in Setting 2 with true parameter values (λ * , γ * , p * , τ * )=(1, 1/3, 0.3, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by KM and MIF. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. It appears that the influence of knowing or not knowing the initial condition i 0 is different according to the values of the parameters used to simulate the data. For the setting where λ = 1 and γ = 1/3, Tables 3 and4 does not exhibit major differences between estimates. On the contrary, the impact of knowing or not knowing the initial condition i 0 is more visible when considering λ = 0.6 and γ = 0.4 (see Appendix E). Tables E.8 and E.9 show that the quality of estimates deteriorates when i 0 is unknown, leading in particular to more significant biases. For N = 10000, n = 101 and p * = 0.8, p passes from 0.77 when i 0 is known to 0.65 when it is not. Once again, higher frequency observations of the epidemics lead to more satisfactory estimates (see Table E.8). Tables 4 and5 suggest that the estimation bias obtained for the measurement error τ increases when p * decreases.

N = 1000 N = 2000 N = 10000 n = 11 n = 31 n = 101 n = 11 n = 31 n = 102 n = 10 n = 30 n = 100 (7,
N = 1000 N = 2000 N = 10000 n = 11 n = 31 n = 101 n = 11 n = 31 n = 102 n = 10 n = 30 n = 100 (7,

Additional comments

In the simulation study, we also considered cases where only the susceptible individuals are observed (not shown here). We noticed that the estimates provided by our Kalman-based method and the MIF algorithm were more accurate when considering the S rather than the I values. As the S values are several orders of magnitude larger than the I ones, a plausible explanation is that the observation noise (due to imperfect reporting and measurement errors) has a lower impact on the S values.

As for the computation times of both methods, these are sensitive to the number of observations n per epidemic: the computation time increases linearly with n. Concerning the population size N, only the computation time for MIF-based inference increased when N increased, while our method was insensitive to it. As an example, for the scenario with N = 10000, n = 30 and p = 0.8 (which corresponds to Table 1), the average computation time for a single estimate (i.e., a single trajectory) was 31 seconds with KM and 97 seconds with the MIF algorithm. For n = 100, the average computation times were 81 and 147 seconds for the KM and MIF algorithms, respectively.

Confidence interval estimates based on profile likelihood

Following other authors (see [START_REF] Ionides | Monte carlo profile confidence intervals for dynamic systems[END_REF] for instance), we provide profile-likelihood confidence intervals of estimated parameters, for which we briefly recall the principle. Let us denote a general parameter vector ψ = (ψ 1 , ψ 2 ), where ψ 1 ∈ R is the parameter of interest and ψ 2 contains the remaining parameters. The profile log-likelihood of ψ 1 is built by maximizing the approximate log-likelihood function (proposed in Section 3) over ψ 2 , for fixed values of ψ 1 : L pro f ile (ψ 1 ) = max ψ 2 L(ψ 1 , ψ 2 ). A 95% confidence interval for ψ 1 is given by:

{ψ 1 : L pro f ile ( ψ) -L pro f ile (ψ 1 )} < 1.92, ( 23
)
where ψ is the maximum approximate likelihood estimator (see ( 19)). The threshold value of 1.92 comes from Wilks' theorem and corresponds to the quantile of order 0.95 of the χ 2 distribution with 1 degree of freedom.

As an illustrative example, 95% profile likelihood confidence intervals were constructed for the key epidemic parameters λ and γ on two particular trajectories of SIR simulated dynamics in the first experiment (τ = 0). A graphical representation is provided in Figure 3 for parameter λ and in Figure 4 for parameter γ. The first confidence interval (left panel of both figures) is obtained with a sample of n = 30 observation of an SIR epidemic for a population of size N = 2000 with reporting rate p = 0.3. The second confidence interval (right panel of each Figures) is obtained with a sample of n = 100 observation of an SIR epidemic for a population of size N = 10000 with reporting rate p = 0.8. For each of the two parameters (playing the role of ψ 1 in (23)), 20 values were considered in a relevant interval containing the point estimate. For each of the 20 values of the parameter of interest, the remaining parameters (playing the role of ψ 2 in (23)), on which the likelihood is optimized (corresponding to L pro f ile (ψ 1 ) in ( 23))), were randomly initialized, with 10 different initialization values, the best being stored. The 20 values of maximum log-likelihood were reported on a graph, linked up by a smoothing curve. The two vertical lines, going through the intersection of this curve with the horizontal line at the y-value equal to the maximum log-likelihood for all parameters minus 1.92 (cf. equation ( 23)), determine the x-value for the CI95%. Based on Figures 3 and4, we see that the widths of the confidence intervals CI95%(λ) = [0.96, 1.10] and CI95%(γ) = [0.31, 0.48] are naturally greater in the case where N = 2000, n = 30 and p = 0.3 (which is a more difficult case for performing estimates, due to an increased stochasticity of epidemic trajectories and significant noise in the observations) than for N = 10000, n = 100 and p = 0.8 (a much more tractable case with low variability amongst trajectories and low levels of noise in observations): CI95%(λ) = [0.95, 1.00] and CI95%(γ) = [0.33, 0.36].

Application on real data

We applied our inference method on the data from an influenza outbreak that occurred in January 1978 in a boarding school in the north of England (Anonymous (1978)), with N = 763. The observations correspond to the daily number of infectious boys across 14 days (n = 14). It is known that the epidemic started from a single infectious student. Here we also assumed that the epidemic dynamics followed an SIR model. Hence, S (0) = 762 and I(0) = 1, and the parameters to be estimated are the epidemic parameters (λ, γ), the reporting rate p, and the parameter τ related to observational noise.

Estimates were performed with both KM and MIF. For the MIF method, we used the same tuning parameters values as those chosen in the simulation study. Both series of results were graphically assessed by post-predictive checks. For this, the Markov jump processes of the SIR model were simulated using each set of parameter estimates. We kept 1000 trajectories that did not exhibit early extinction, according to the theoretical criterion used in Section 4.1. From these 1000 trajectories, we then generated equally-spaced observations with n = 14. Empirical mean, 5th, 50th and 95th percentiles were extracted at each time point and superimposed on the real data (Figure 5).

The following estimates were obtained, with the profile likelihood-based confidence intervals (CI95%) provided in brackets:

• λKM = 1. 72 [1.61, 1.83]; γKM = 0.48 [0.43, 0.52]; pKM = 1.00 [0.92, 1.00]; τKM = 0.91 [0.42, 1.62] with KM,

• λMIF = 1. 85 [1.62, 2.15]; γMIF = 0.47 [0.39, 0.54]; pMIF = 0.97 [0.84, 1.00]; τMIF = 1.58 [0.80, 2.80] with MIF.

The estimated values for λ, γ and p are similar in both methods, but the estimated values for τ are rather different. The confidence intervals provided by the MIF method are larger than those obtained by our Kalman-based method, but this could be due to non-optimal tuning in the MIF case. Moreover, we see that the confidence interval for τ is particularly wide for both methods, which is in agreement with the fact that a moderate number of observations is needed in order to properly estimate τ (as showed in the simulation analyses). A post-predictive check (Figure 5) indicates that both methods provide estimates and hence predictions that are consistent with the data. Estimation took 22.7 seconds with our method, versus 46.5 seconds using MIF.

Discussion

In this paper we have proposed a general and practical inference method for continuoustime epidemics involving discrete, partially and noisily observed time-series data. We derived a Gaussian approximation of an epidemic's density-dependent Markovian jump process underlying its dynamics using a diffusion based approach and a Gaussian approximation of observations model. This two-level Gaussian approximation allowed us to develop an inference method based on Kalman filtering for the calculation of the likelihood, to estimate key epidemic parameters (such as transmission and recovery rates), the initial state of the system (number of susceptible and infectious individuals), and parameters of the observation model (such as the reporting rate) from incomplete and noisy data (proportion of infectious individuals over time).

The performance of the estimators obtained with the Kalman-based method was investigated on simulated data under various scenarios with respect to the parameter values of epidemic and observation processes, the population size (N), the number of observations (n), and the nature of the data (number of susceptible S or infectious I individuals over time). Performance, in terms of bias and in particular accuracy, improved when increasing N and (especially) n, and was satisfactory for a realistic observation design (e.g., n = 30, which corresponds in our case to one observation per day or every two days) and moderate community size (N = 2000).

The influence of N and n is less pronounced when data are more complete, here in the case where p, the proportion of available data-corresponding to the reporting rate-was equal to 0.8, and τ, corresponding to the measurement error, was zero. Estimation was more challenging when the measurement error τ was unknown. In the latter case, higher frequency observations were needed in order to obtain more accurate estimates. When, in addition to a non-zero measurement error, the initial point i 0 is unknown, the quality of the estimates could deteriorate in some cases.

A similar performance was observed irrespective of data type (when observations were sampled from S instead of I; results not shown). In addition, our method seemed to be little-impacted by tuning aspects. Indeed, the only obvious tuning parameter, concerning the initialization of the covariance matrix of the state variables conditionally upon the observations-in the filtering step-did not seem to influence estimation accuracy. Besides simulated data, our method provided quite plausible estimates when applied to real data from an influenza outbreak in a British boarding school in 1978, supported by the fact that the post-predictive check showed consistency with data. The good performance seen here is all the more noteworthy given that the data came with certain difficulties (low N and n).

Estimates obtained with KM were compared to those using MIF [START_REF] Ionides | Iterated filtering[END_REF], [START_REF] King | Statistical inference for partially observed markov processes via the r package pomp[END_REF]). The MIF algorithm is efficient in terms of inference quality, but computationally expensive and uses tuning parameters (number of particles, number of iterations, etc.) that are crucial for the successful functioning of the procedure. Importantly, our method does not require such specific computational calibration and its results are computed faster.

In terms of limitations of our method, we observed that the joint estimation of parameters from epidemic and observation models (λ, γ, p), along with the initial conditions of the underlying epidemic process (proportions of susceptible and infectious individuals (s 0 , i 0 )), can lead to difficulties when no constraint (e.g., s 0 + i 0 = 1) is imposed, and when only one discretized and perturbed coordinate of the system (here I) is observed. This occurred even in a "simple" scenario where N = 10000, n = 100, and p = 0.8 (low stochasticity and little loss of information in the data). This difficulty is no longer encountered if the two coordinates of the system (S and I) are observed. As well as this issue, two blocks of dependance between estimates were observed: (λ, γ, s 0 ) on the one hand, and (p, i 0 ) on the other. Therefore, an incorrect estimate of i 0 or s 0 will be reflected in the estimate of p and (λ, γ), respectively. One potential way to solve this problem could be to consider a prior for the initial conditions of the system. For more details on how to overcome this practical issue, see [START_REF] Stocks | Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in germany[END_REF], [START_REF] Stocks | Iterated filtering methods for markov process epidemic models[END_REF], who also emphasize the fact that inference algorithms are very sensitive to the initial values of the system.

Our method relies on two successive Gaussian model approximations (one for the latent state and the other for the observation model). These approximations do not seem to alter the quality of the estimates. Indeed, the small variance coefficient N -1/2 provides an advantageous framework for the approximation of the state model, for which the Kalman filter performs very well in practice (small prediction errors). The decent accuracy of Gaussian process approximations for stochastic epidemic models has previously been highlighted [START_REF] Buckingham-Jeffery | Gaussian process approximations for fast inference from infectious disease data[END_REF]). Here, we went further and examined the performance of Gaussian approximations of epidemic dynamics, not only by using a different approach based on Kalman filtering, but also by considering an even less convenient configuration where the initial conditions and observation errors had to be estimated.

Our approach can be generalized in several ways. First, although we focused in this study on the SIR model as a case study, our method is quite general since it can be extended to other mechanistic models of epidemic dynamics, including additional health states (such as an exposed state E). Second, the observations can encompass variable sampling intervals (i.e., ∆, the time step between two consecutive observations, is not necessarily constant). Third, other types of observations can be considered, both with regards to their nature (e.g., the number of new infectious individuals, which can be viewed as a function of state variables S and I) and to the error model.

Therefore, given its ease in implementation, low computation time, and satisfactory performance, we recommend the use of our Kalman filtering-based estimation method to providing an initial guess for parameters in the framework of partially observed complex epidemic dynamics. Table E.6: First experiment (τ = 0). Estimation of θ = (λ, γ, p, i 0 ) under the constraint s 0 + i 0 = 1 in Setting 1 with true parameter values (λ * , γ * , p * , i * 0 )=(0.6, 0.4, 0.8, 0.01). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by our Kalman-based method (KM) and maximum iterated filtering (MIF). The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. E.7: First experiment (τ = 0). Estimation of θ = (λ, γ, p, i 0 ) under the constraint s 0 + i 0 = 1 in Setting 2 with true parameter values (λ * , γ * , p * , i * 0 )=(0.6, 0.4, 0.3, 0.01). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by our Kalman-based method (KM) and maximum iterated filtering (MIF). The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. The results on the second set of epidemic parameters displayed in Tables E.6 and E.7 are more contrasted, since the parameter values chosen (λ * = 0.6 and γ * = 0.4) generate more stochasticity (see Figure E.6), so trajectories are less similar and further from the mean of the jump process; hence estimates are less accurate. Besides, the peak of the number of infectious individuals is clearly lower than in the λ * = 1 and γ * = 1/3 case. The estimates of p are particularly poor when n is low, which obviously impacts estimation of the other parameters.

Table E.8: Second experiment (τ 0). Estimation of θ = (λ, γ, p, i 0 , τ) under the constraint s 0 + i 0 = 1 in Setting 1 with true parameter values (λ * , γ * , p * , i * 0 , τ * )=(0.6, 0.4, 0.8, 0.01, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by our Kalman-based method and the MIF algorithm. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. Known starting point i 0 . Tables E.9 and E.10 respectively display the results obtained by our Kalman-based method and the MIF algorithm for the high-reporting scenario (p = 0.8) and low-reporting scenario (p = 0.3).

Table E.9: Second experiment (τ 0). Estimation of θ = (λ, γ, p, τ) with s 0 = 0.99 and i 0 = 0.01 known in Setting 1 with true parameter values (λ * , γ * , p * , τ * )=(0.6, 0.4, 0.8, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by KM and MIF. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets. E.10: Second experiment (τ 0). Estimation of θ = (λ, γ, p, τ) with s 0 = 0.99 and i 0 = 0.01 known in Setting 2 with true parameter values (λ * , γ * , p * , τ * )=(0.6, 0.4, 0.3, 0.5). For each combination of (N, n) and for each model parameter, point estimates and standard deviations are calculated as the mean of the 500 individual estimates and their standard deviations (in brackets) obtained by KM and MIF. The reported values for the number of observations n correspond to the average over the 500 trajectories, with the min and max in brackets.

Figure 1 :

 1 Figure 1: SIR compartmental model with three blocks corresponding respectively to susceptible (S), infectious (I), and recovered (R) individuals. Transitions of individuals from state S to I are governed by the transmission rate λ, and transitions of individuals from state I to R are governed by the recovery rate γ of the epidemic.

Figure 2 :

 2 Figure 2: Left panel: ODE solution for the number of infected individuals I (plain black line) and 20 trajectories of the Markov jump process for I (purple lines) when N = 2000. Right panel: n = 30 observations obtained from a particular trajectory of the jump process (in bold purple in the left panel) as a function of time. The points and triangles stand for observations generated with measurement error terms τ = 0 and τ = 0.5 respectively, and the blue and red symbols represent observations generated with p = 0.8 and p = 0.3 respectively.

4. 3 .

 3 Point estimates and standard deviations for key model parameters θ 4.3.1. Simulation results for the first experiment (τ = 0) Three different target values for sample sizes were considered: n = 10, n = 30 and n = 100. Tables

  n = 11 n = 31 n = 101 n = 11 n = 31 n = 102 n = 10 n = 30 n = 100 (7, 18) (21, 51) (

  n = 11 n = 31 n = 101 n = 11 n = 31 n = 102 n = 10 n = 30 n = 100 (7

  n = 11 n = 31 n = 101 n = 501 n = 1001 n = 10 n = 30 n = 100 n = 499 n =

Figure 3 :

 3 Figure 3: Profile likelihood and confidence intervals (CI95%) for λ. Left panel: data simulated with N = 2000, n = 30 and p = 0.3; the true value λ * = 1, the point estimate λ = 1.02, and CI95% = [0.96, 1.10]. Right panel: data simulated with N = 10000, n = 100 and p = 0.8; the true value λ * = 1, the point estimate λ = 1.00, and CI95% = [0.95, 1.00].

Figure 4 :

 4 Figure 4: Profile likelihood and confidence intervals (CI95%) for γ. Left panel: data simulated with N = 2000, n = 30 and p = 0.3; the true value γ * = 1/3, the point estimate γ = 0.32, and CI95% = [0.31, 0.48]. Right panel: data simulated with N = 10000, n = 100 and p = 0.8; the true value γ * = 1/3, the point estimate γ = 0.34, and CI95% = [0.33, 0.36].

Figure 5 :

 5 Figure 5: Post-predictive checks for the Kalman-based (KM, left panel) and the maximum iterated filtering (MIF, right panel) estimates. In blue: observations (number of infectious boys). Solid red line: average trajectory over 1000 Markov jump processes from the estimated model. Dotted red lines: 5th, 50th, and 95th percentiles.

Figure E. 6 :

 6 Figure E.6: Left panel: ODE solution for the number of infected individuals I (plain black line) and 20 trajectories of the Markov jump process for I (purple lines) when N = 2000. Right panel: n = 30 observations obtained from a particular trajectory of the jump process (in bold purple in the left panel) as a function of time. The points and triangles stand for observations generated with measurement error terms τ = 0 and τ = 0.5 respectively, and the blue and red symbols represent observations generated with p = 0.8 and p = 0.3 respectively.

  n = 10 n = 30 n = 99 n = 11 n = 31 n = 102 n = 11 n = 31 n = 101 (3

  n = 10 n = 30 n = 99 n = 11 n = 31 n = 102 n = 11 n = 31 n = 101 (3

  n = 10 n = 30 n = 99 n = 499 n = 998 n = 11 n = 31 n = 101 n = 500 n = 1001 (3

  n = 10 n = 30 n = 99 n = 11 n = 31 n = 102 n = 11 n = 31 n = 101 (3

  n = 10 n = 30 n = 99 n = 11 n = 31 n = 102 n = 11 n = 31 n = 101 (3

  confidence intervals Figures E.7 and E.8 represent the profile likelihoods and the subsequent confidence intervals (CI95%) for the parameters λ and γ obtained for our Kalman filtering-based method in two settings (first case: N = 2000, n = 30, and p = 0.3; second case: N = 10000, n = 100, and p = 0.8).

Figure E. 7 :

 7 Figure E.7: Profile likelihood and confidence intervals (CI95%) for λ. Left panel: N = 2000, n = 30, and p = 0.3. The true value λ * = 0.6, the point estimate λ = 0.47, and CI95% = [0.54, 0.76]. Right panel: N = 10000, n = 100, and p = 0.8. The true value λ * = 0.6, the point estimate λ = 0.60, and CI95% = [0.56, 0.64].

Figure E. 8 :

 8 Figure E.8: Profile likelihood and confidence intervals (CI95%) for γ. Left panel: N = 2000, n = 30, and p = 0.3. The true value γ * = 0.4, the point estimate γ = 0.21, and CI95% = [0.26, 0.48]. Right panel: N = 10000, n = 100, and p = 0.8. The true value γ * = 0.4, the point estimate γ = 0.40, and CI95% = [0.35, 0.42].

  , 19) (9, 56) (30, 182) (7, 19) (20, 56) (66, 182) (8, 17) (24, 49) (78, 160) λ

	* = 0.6
	KM	0.47	0.57	0.57	0.47	0.54	0.58	0.56	0.55	0.56
		(0.16) (0.13) (0.16)	(0.12) (0.09)	(0.09)	(0.06) (0.06)	(0.05)
	MIF	0.51	0.55	0.55	0.52	0.55	0.56	0.57	0.57	0.59
		(0.12) (0.13) (0.11)	(0.09) (0.07)	(0.07)	(0.05) (0.04)	(0.03)
	γ * = 0.4
	KM	0.20	0.29	0.35	0.26	0.30	0.37	0.37	0.35	0.35
		(0.12) (0.13) (0.14)	(0.10) (0.09)	(0.09)	(0.05) (0.06)	(0.05)
	MIF	0.22	0.26	0.32	0.25	0.30	0.34	0.38	0.37	0.39
		(0.08) (0.10) (0.09)	(0.10) (0.08)	(0.06)	(0.06) (0.04)	(0.03)
	p * = 0.3
	KM	0.10	0.19	0.24	0.15	0.20	0.27	0.28	0.25	0.25
		(0.09) (0.15) (0.12)	(0.08) (0.12)	(0.13)	(0.06) (0.08)	(0.10)
	MIF	0.11	0.14	0.19	0.14	0.18	0.21	0.27	0.27	0.28
		(0.06) (0.09) (0.06)	(0.08) (0.07)	(0.05)	(0.07) (0.04)	(0.03)
	τ * = 0.5
	KM	0.23	0.42	0.52	0.13	0.24	0.51	0.05	0.19	0.36
		(0.20) (0.29) (0.21)	(0.19) (0.21)	(0.18)	(0.14) (0.15)	(0.12)
	MIF	0.38	0.37	0.44	0.39	0.29	0.43	0.61	0.33	0.43
		(0.14) (0.15) (0.09)	(0.15) (0.15)	(0.07)	(0.17

Acknowledgments

We thank two anonymous referees for their constructive and helpful comments.

Funding

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

Appendix A. Remarks on the sampling interval

The sampling interval ∆ is important in our method and we distinguish between two cases: "Small ∆" and "Moderate ∆". We give below the dependencies on quantities of interest with respect to ∆.

(1) Small sampling interval ∆ Taylor expansions with respect to t at point t k-1 yield

The following additional approximations, which simplify the analytic expressions, can be used in the state space equation:

(2) Moderate ∆ Computing the approximate log-likelihood (20) with Kalman filtering techniques requires computing the resolvent matrix Φ of the ODE system (8). When the time intervals between observations are too large (i.e., ∆ is too large), we use the following approximation for matrix exponentials:

where t k = a 1 < a 2 < . . . < a J = t k+1 . This can however significantly increase computation times.

Appendix B. Proof of Proposition 2

By the semigroup property of Φ, we have that g, defined in (7), satisfies for s ≤ t,

Substituting g(s) with √ N(G N (s)x(s)) using (9) yields:

By the independent increments property of Brownian motion, we get moreover that U k is independent of F t k-1 . This achieves the proof of Proposition 2.

Appendix C. Proof of Lemma 2

Assume first that Q and T are non-singular. The joint distribution of (Y, X) is Gaussian:

Hence,

we get:

We then obtain, using the matrix relation:

the following results:

Appendix D. Proof of Proposition 3

For k = 0, we have that X 0 ∼ N(ξ 0 , Ξ0 ). The induction assumption is:

To get (i), we apply Lemma 2, noting that the distribution L(X k |Y k-1,0 ) = N d ( Xk , Ξk ) and that the distribution Y k conditional on X k is N(BX k , Q k ). Therefore, setting ξ = Xk , T = Ξk and Q = Q k , we get that the distribution of (X k |Y k:0 ) is N d ( Xk , Tk ), with Xk = ξ(Y k ), where ξ(Y k ) and Tk are given by ( 21). These are precisely the expressions for Tk and Ξk given in (i).

For (ii), we use that

For (iii), we use that

. The induction assumption is fulfilled and therefore this achieves the proof of Proposition 3.

Appendix E. Additional simulation study

Appendix E.1. Description

We reproduced the simulation study described in Section 4 with other parameter values: λ = 0.6, γ = 0.4, s 0 = 0.99, i 0 = 0.01. An extract of the simulated data is shown in Tables E.6 and E.7 respectively display the results for the high-reporting scenario (p = 0.8) and low reporting scenario (p = 0.3) when τ = 0 and is not estimated. Each table compares the Kalman-based method (KM) to the maximum iterated filtering algorithm (MIF). The first column display the true parameter values. Columns 2 to 10 display the results for different combinations of (N, n). For each parameter and each estimation method, the reported values are the mean of the 500 parameter estimates and their standard deviations (in brackets).

Appendix E.2.2. Numerical results for the second experiment (τ 0) Unknown starting point i 0 . Table E.8 displays the results obtained by our Kalman-based method and the MIF algorithm for the high-reporting scenario (p = 0.8).

Appendix F. User-friendly code

We propose user-friendly code composed of four distinct programs in the R language, available at the RunMyCode website: http://www.runmycode.org/companion/view/4074.

• KalmanFunctions.R includes general functions implementing the Kalman filter and computing the likelihood of the observations, given a specified compartmental model, with a fixed sampling interval. These functions are easily generalizable to the case where the sampling interval is variable. Moreover, this script includes a function computing the resolvent matrix defined in (8) for large time intervals between observations ∆.

• ModelFunctions.R implements the SIR and SEIR models and defines the key quantities (described in the manuscript for the SIR model) necessary to apply the Kalman filterbased method. More precisely, given a compartmental model (SIR or SEIR), the following functions are implemented: the ode system, the drift function, the gradient of the drift function, the diffusion matrix, the projection operator linking the observations to the states of the epidemic model and the variance of the observations.

• SIRexample.R and SEIRexample.R simulate respectively SIR and SEIR Markovian jump processes for a set of parameters values, using the GillespieSSA package. The observations of infectious individuals are obtained by: O 1 (t k ) ∼ Binomial(I(t k ), p), O 2 (t k ) ∼ N(0, τ 2 I(t k )), k = 1, . . . , n, at regularly-spaced time points. Finally, an estimation of key parameters λ, γ, p and τ with known starting points and, in the SEIR model, with a known transition rate from E to I, is proposed.