## 4.2 Modèle et propriété

### 4.2.1 Définition du modèle

Pour décrire notre nouveau modèle exposé-infecté, nous avons repris les hypothèses du modèle précédent (3.2.1), tout en introduisant un nouveau compartiment d'individus exposés E. Ces individus ne sont pas contagieux, mais peuvent devenir infectés à un taux individuel $\alpha$ (taux d'incubation). Les individus infectés et contagieux I sont considérés comme isolés à partir de la date d'hospitalisation (Fig.4.1).

La dynamique de la maladie est caractérisée par une phase d'incubation, de contagiosité et d'hospitalisation. Pendant la phase d'incubation, l'individu exposé (il a été en contact avec la maladie mais n'est pas encore infectieux) ne présente pas encore de symptômes de la maladie, et il devient infectieux après une durée moyenne $1/\alpha$, où $\alpha > 0$ est le taux d'incubation individuel. Ensuite, il y a une phase contagieuse où l'individu infecté devient contagieux et peut transmettre la maladie à un taux $\lambda > 0$ à une autre personne. La phase d'hospitalisation est caractérisée par l'isolement de l'individu à un taux $\mu > 0$, où il ne peut plus contaminer d'autres personnes. Pendant cette phase, l'individu reçoit un traitement à base d'antibiotiques et peut se rétablir en l'espace de 1 à 5 jours, en fonction de la gravité de l'infection et de la réponse du patient au traitement [71]. La contamination via l'environnement génère également de nouveaux infectés via un processus de Poisson de taux $\nu > 0$ (Fig.4.2).

<figure><img src="image_34.png" /><figcaption>FIGURE 4.2 – Dynamique de la fièvre typhoïde à Mayotte dans le modèle à deux compartiments.</figcaption></figure>

Le modèle peut être caractérisé par un processus de Markov bivarié en temps continu $(E_t, I_t)_{t \ge 0}$, à valeurs dans $S = \mathbb{N} \times \mathbb{N}$, où $E_t$ et $I_t$ représentent respectivement le nombre des individus exposés et infectieux à l'instant $t$ au sein de la population.

Étant donné l'état courant $(e, i)$, les transitions possibles du notre processus sont :

* $(e, i) \to (e + 1, i)$ avec un taux de $\lambda i + \nu$.

* $(e, i) \to (e - 1, i + 1)$ avec un taux de $\alpha e$.

* $(e, i) \to (e, i - 1)$ avec un taux de $\mu i$.

Ainsi, le générateur infinitésimal $Q$ du processus $(E_t, I_t)_{t \ge 0}$ est donné par :

$$ 
\begin{aligned}
Q_{(e,i),(e+1,i)} &= \lambda i + \nu, && \text{pour } e, i \ge 0, \\
Q_{(e,i),(e-1,i+1)} &= \alpha e, && \text{pour } e \ge 1, i \ge 0, \\
Q_{(e,i),(e,i-1)} &= \mu i, && \text{pour } e \ge 0, i \ge 1, \\
Q_{(e,i),(e,i)} &= -((\lambda + \mu)i + \alpha e + \nu), && \text{pour } e, i \ge 0,
\end{aligned} 
\quad (4.1) 
$$
où les paramètres $\lambda, \nu, \alpha$ et $\mu$ sont positifs.

Nous supposons que $(E_0, I_0) = (e, i) \in \mathbb{S}$ et nous désignons par

$$p_{(e,i)(e',i')}(t) = \mathbb{P}((E_t, I_t) = (e', i') | (E_0, I_0) = (e, i)) \quad (e', i') \in \mathbb{S}, t \ge 0,$$

les probabilités de transition de $(E_t, I_t)_{t \ge 0}$. Elles satisfont l'équation de Kolmogorov progressive suivante :

$$\frac{d}{dt} p_{(e,i),(e',i')}(t) = (\lambda i' + \nu) p_{(e,i),(e'-1,i')}(t) + \alpha(e' + 1)p_{(e,i),(e'+1,i'-1)}(t) + \mu(i' + 1)p_{(e,i),(e',i'+1)}(t) \\ - ((\lambda + \mu)i' + \alpha e' + \nu) p_{(e,i),(e',i')}(t), \quad (4.2)$$

avec la condition initiale $p_{(e,i),(e,i)}(0) = 1$ et $p_{(e,i),(e',i')}(0) = 0$ pour $(e', i') \neq (e, i)$.

Dans la sous-section suivante, notre intérêt se porte sur l'étude théorique du comportement moyen et asymptotique du processus $(E_t, I_t)_{t \ge 0}$. Nous allons explorer les propriétés mathématiques de ce modèle et comprendre comment il se comporte à long terme. Nous allons également illustrer l'évolution de ce modèle au cours du temps grâce à des simulations numériques. Cette analyse approfondie nous permettra de prédire l'évolution de l'épidémie dans le temps en fonction des paramètres.

### 4.2.2 Propriétés mathématiques du modèle

Nous étudions maintenant quelques propriétés mathématiques du processus $(E_t, I_t)_{t \ge 0}$.

**Théorème 4.1.** (Comportement moyen de $(E_t, I_t)_{t \ge 0}$). Pour tout $e_0, i_0 \ge 0$,

*si* $\lambda = \mu$, alors

$$\left\{ \begin{aligned} \mathbb{E}_{(e_0,i_0)}[E_t] &= \frac{\lambda\nu}{\lambda+\alpha}t - c\left(e^{-(\lambda+\alpha)t} + \frac{\lambda}{\alpha}\right) + i_0\frac{\lambda}{\alpha} + \frac{\nu}{\lambda+\alpha}, \\ \mathbb{E}_{(e_0,i_0)}[I_t] &= \frac{\alpha\nu}{\lambda+\alpha}t + c\left(e^{-(\lambda+\alpha)t} - 1\right) + i_0, \end{aligned} \right.$$

et *si* $\lambda \neq \mu$, alors

$$\left\{ \begin{aligned} \mathbb{E}_{(e_0,i_0)}[E_t] &= c_1 e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_2 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\mu\nu}{\alpha(\mu-\lambda)}, \\ \mathbb{E}_{(e_0,i_0)}[I_t] &= c_3 e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_4 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\nu}{\mu-\lambda}, \end{aligned} \right.$$

où

$$c = \left(e_0 + i_0 - \frac{\nu}{\lambda + \alpha}\right) \frac{\alpha}{\lambda + \alpha},$$

$$c_1 = \frac{1}{a} \left( i_0 + \frac{\mu - \alpha + a}{2} e_0 - \frac{\nu}{\mu - \lambda} \left( 1 + \frac{\mu(\mu - \alpha + a)}{2\alpha} \right) \right),$$

$$c_2 = e_0 - c_1 - \frac{\mu\nu}{\alpha(\mu - \lambda)},$$

$$c_3 = \frac{-\mu + \alpha + a}{2a} \left( i_0 + \frac{\mu - \alpha + a}{2} e_0 - \frac{\nu}{\mu - \lambda} \left( 1 + \frac{\mu(\mu - \alpha + a)}{2\alpha} \right) \right),$$

$$c_4 = -\frac{\mu - \alpha + a}{2} c_2,$$

$$a = \sqrt{(\mu - \alpha)^2 + 4\alpha\lambda}.$$
Preuve. Soit $(e_0, i_0) \in \mathbb{S}$, pour toute fonction $f$ positive ou bornée sur $\mathbb{S}$ et pour tout $e, i, \ge 0$, nous avons :

$$ \mathbb{E}_{(e_0, i_0)}[f(E_t, I_t)] = \sum_{(e,i) \in \mathbb{S}} p_{(e_0, i_0)(e,i)}(t)f(e, i). $$

En utilisant l'équation de Kolmogorov progressive, on obtient :

$$
\begin{align*}
\frac{d}{dt} \mathbb{E}_{(e_0, i_0)}[f(E_t, I_t)] &= \sum_{(e,i) \in S} p'_{(e_0, i_0)(e,i)}(t) f(e, i) \\
&= \sum_{(e,i) \in \mathbb{S}} \left( \sum_{(e',i') \in \mathbb{S}} p_{(e_0, i_0)(e',i')} Q_{(e',i')(e,i)} \right) f(e, i) \\
&= \sum_{(e',i') \in \mathbb{S}} p_{(e_0, i_0)(e',i')} \left( \sum_{(e,i) \in \mathbb{S}} Q_{(e',i')(e,i)} f(e, i) \right) \\
&= \sum_{(e',i') \in S} p_{(e_0, i_0)(e',i')} \left( (\lambda I_t + \nu) (f(e' + 1, i') - f(e', i')) \right. \\
&\qquad \left. + \alpha e' (f(e' - 1, i' + 1) - f(e', i')) + \mu i' (f(e', i' - 1) - f(e', i')) \right) \\
&= \mathbb{E}_{(e_0, i_0)} \left[ (\lambda I_t + \nu) (f(E_t + 1, I_t)) - f(E_t, I_t) \right) \\
&\qquad + \alpha E_t (f(E_t - 1, I_t + 1) - f(E_t, I_t)) + \mu I_t (f(E_t, I_t - 1) - f(E_t, I_t)).
\end{align*}
$$

Pour $f(E_t, I_t) = E_t$, nous trouvons :

$$
\begin{align}
\frac{d}{dt} \mathbb{E}_{(e_0, i_0)}[E_t] &= \mathbb{E}_{(e_0, i_0)}[\lambda I_t + \nu - \alpha E_t] \nonumber \\
&= -\alpha \mathbb{E}_{(e_0, i_0)}[E_t] + \lambda \mathbb{E}_{(e_0, i_0)}[I_t] + \nu. \tag{4.3}
\end{align}
$$

Pour $f(E_t, I_t) = I_t$, nous avons :

$$
\begin{align}
\frac{d}{dt} \mathbb{E}_{(e_0, i_0)}[I_t] &= \mathbb{E}_{(e_0, i_0)}[\alpha E_t - \mu I_t] \nonumber \\
&= \alpha \mathbb{E}_{(e_0, i_0)}[E_t] - \mu \mathbb{E}_{(e_0, i_0)}[I_t]. \tag{4.4}
\end{align}
$$

En combinant les deux équations **4.3** et **4.4**, nous obtenons le système d'équations différentielles suivant :

$$
\left\{
\begin{aligned}
\frac{d}{dt} \mathbb{E}_{(e_0, i_0)}[E_t] &= -\alpha \mathbb{E}_{(e_0, i_0)}[E_t] + \lambda \mathbb{E}_{(e_0, i_0)}[I_t] + \nu, \\
\frac{d}{dt} \mathbb{E}_{(e_0, i_0)}[I_t] &= \alpha \mathbb{E}_{(e_0, i_0)}[E_t] - \mu \mathbb{E}_{(e_0, i_0)}[I_t].
\end{aligned}
\right.
\qquad (4.5)
$$

En résolvant le système (4.5), nous pouvons obtenir les expressions des espérances suivantes :

Si $\lambda = \mu$, alors

$$
\left\{
\begin{aligned}
\mathbb{E}_{(e_0, i_0)}[E_t] &= \frac{\lambda \nu}{\lambda + \alpha} t - c \left( e^{-(\lambda+\alpha)t} + \frac{\lambda}{\alpha} \right) + i_0 \frac{\lambda}{\alpha} + \frac{\nu}{\lambda + \alpha}, \\
\mathbb{E}_{(e_0, i_0)}[I_t] &= \frac{\alpha \nu}{\lambda + \alpha} t + c \left( e^{-(\lambda+\alpha)t} - 1 \right) + i_0.
\end{aligned}
\right.
$$
Et si $\lambda \neq \mu$, alors

$$ egin{cases} 
\mathbb{E}_{(e_0,i_0)}[E_t] = c_1 e^{\frac{-\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_2 e^{\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\mu\nu}{\alpha(\mu-\lambda)}, \\ 
\mathbb{E}_{(e_0,i_0)}[I_t] = c_3 e^{\frac{-\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_4 e^{\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\nu}{\mu-\lambda}. 
\end{cases} $$

Les constantes réelles $c$ et $c_i$ pour $i = 1, \dots, 4$ sont déterminées à partir de la condition initiale $(E_0, I_0) = (e_0, i_0)$.

□

**Corollaire 4.1.** Asymptotiquement et conditionnellement à $(E_0, I_0) = (e_0, i_0)$, dans le cas où $\lambda \ge \mu$, le comportement moyen du processus $(E_t, I_t)_{t \ge 0}$ est explosif. Lorsque $\lambda > \mu$, les nombres moyens des individus exposés et infectés augmentent de manière exponentielle (Fig.4.3). Cela indique une croissance non contrôlée de l'épidémie, où le nombre d'infections continue d'augmenter rapidement au fil du temps. Si $\lambda = \mu$, les deux espérances augmentent de manière linéaire (Fig.4.4). c'est-à-dire que le nombre moyen d'individus exposés et infectés augmente proportionnellement au temps. Cependant, le processus moyen se stabilise dans le cas où $\lambda < \mu$ (Fig.4.5), ce qui conduit à une stabilisation du nombre moyen d'individus exposés et infectés à un niveau constant donné par :

$$ egin{cases} 
\lim_{t \to +\infty} \mathbb{E}_{(e_0,i_0)}[E_t] = \frac{\mu\nu}{\alpha(\mu-\lambda)}, \\ 
\lim_{t \to +\infty} \mathbb{E}_{(e_0,i_0)}[I_t] = \frac{\nu}{\mu-\lambda}. 
\end{cases} $$

Cette situation correspond à un équilibre endémique des moyennes.

<figure><img src="image_35.png" /><figcaption>FIGURE 4.3 – Évolution du processus $(E_t, I_t)_{t \ge 0}$ et sa moyenne $(\mathbb{E}[E_t], \mathbb{E}[I_t])_{t \ge 0}$ pour les paramètres $\lambda = 0.27, \mu = 0.25, \alpha = 0.26$ et $\nu = 0.02$.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>FIGURE 4.4 – Évolution du processus $(E_t, I_t)_{t \ge 0}$ et sa moyenne $(\mathbb{E}[E_t], \mathbb{E}[I_t])_{t \ge 0}$ pour les paramètres $\lambda = \mu = \alpha = 0.2$ et $\nu = 0.02$.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 4.5 – Évolution du processus $(E_t, I_t)_{t \ge 0}$ et sa moyenne $(\mathbb{E}[E_t], \mathbb{E}[I_t])_{t \ge 0}$ pour les paramètres $\lambda = 0.048$, $\mu = 0.2$, $\alpha = 0.1$ et $\nu = 0.02$.</figcaption></figure>

Dans la suite, nous allons identifier le régime récurrent positif du modèle et étudier l’existence d’une loi invariante pour le processus exposé-infecté. Pour étudier cette propriété, nous allons utiliser l’approche de Lyapunov, qui a été développée pour étudier la stabilité des systèmes dynamiques [72].

Le théorème suivant présente une condition suffisante pour garantir l’existence d’une loi invariante pour le processus $(E_t, I_t)_{t \ge 0}$.
**Théorème 4.2.** (Comportement asymptotique de $(E_t, I_t)_{t \ge 0}$). Si $\lambda < \alpha$ et $\frac{\lambda + \alpha}{2} < \mu$, alors $(E_t, I_t)_{t \ge 0}$ est récurrent positif et admet une unique loi invariante.

*Preuve.* Le processus $(E_t, I_t)_{t \ge 0}$ est clairement irréductible, son générateur $\mathcal{L}$ est donné par : pour toute fonction $f : \mathbb{N}^2 \to \mathbb{R}$

$$ \mathcal{L}f(e, i) = (\lambda i + \nu) [f(e+1, i) - f(e, i)] + \alpha e[f(e-1, i+1) - f(e, i)] \\ + \mu i[f(e, i-1) - f(e, i)]. $$

On applique ce générateur $\mathcal{L}$ à la fonction $f(e, i) = e^2 + i^2$, on obtient :

$$ \mathcal{L}f(e, i) = (\lambda i + \nu) [(e+1)^2 - e^2] + \alpha e[(e-1)^2 - e^2] + \alpha e[(i+1)^2 - i^2] \\ + \mu i[(i-1)^2 - i^2]. $$

Par conséquent, on peut facilement voir qu'il existe un $(e_0, i_0)$, tel que lorsque $(e, i) \ge (e_0, i_0)$,

$$ \begin{aligned} \mathcal{L}f(e, i) &\le 2\lambda e i - 2\alpha e^2 + 2\alpha e i - 2\mu i^2 \\ &\le (\lambda - \alpha)e^2 + (\lambda + \alpha - 2\mu)i^2 \\ &\le \max\{\lambda - \alpha, \lambda + \alpha - 2\mu\} (e^2 + i^2), \end{aligned} $$

où $\max\{\lambda - \alpha, \lambda + \alpha - 2\mu\} < 0$, puisque $\lambda < \alpha$ et $\frac{\lambda + \alpha}{2} < \mu$.

Ainsi,

$$ \mathcal{L}f(e, i) \le \max\{\lambda - \alpha, \lambda + \alpha - 2\mu\} (e^2 + i^2) + C1_{\{(e,i) \ge (e_0, i_0)\}}, $$

où $C > 0$.

Par conséquence, le processus $(E_t, I_t)_{t \ge 0}$ satisfait la condition de Foster-Lyapunov [72, Théorème 7.1], donc, la chaîne de Markov $(E_t, I_t)$ admet une unique loi invariante et converge en $\mathbb{L}^2$ lorsque $t$ tend vers l'infini. $\square$

Après avoir présenté notre modèle exposé-infecté ainsi que ses propriétés mathématiques et identifié et illustré les différents régimes possibles du processus à l'aide de simulations, notre prochaine étape consiste à intégrer le schéma d'observation dans ce modèle exposé-infecté afin de construire une procédure d'estimation.

## 4.3 Estimation des paramètres

Dans cette section, nous présentons notre méthode d'estimation des paramètres du modèle exposé-infecté en adaptant la démarche réalisée tout au long du chapitre précédent **3**, qui ne s'applique pas directement cette fois-ci. Tout d'abord, nous décrivons les caractéristiques de nos observations liées au processus des exposés et infectés **4.3.1**. Nous expliquons ensuite notre stratégie d'estimation ainsi que les différentes étapes impliquées dans cette méthode **4.3.2**. Nous présentons également les expressions analytiques des estimateurs obtenus à partir des moments du processus $(E_t, I_t)_{t \ge 0}$ **4.3.3**. Nous adaptons l'algorithme de Baum-Welch à notre schéma d'observations spécifique, ce qui nous permet d'estimer les probabilités de transition de la chaîne cachée couplée à partir des données cumulées des cas isolés et de déduire les différentes estimations des moments de la chaîne couplée **4.3.4**. Enfin, nous obtenons des estimations des paramètres originaux par "plug-in".
### 4.3.1 Observations

Les nombres d'individus exposés $E_t$ et infectés $I_t$ ne sont pas directement observables. Nous avons à notre disposition uniquement les cumuls de nouveaux cas déclarés par jour, correspondant aux sauts d'amplitude $-1$ du processus $I_t$. Nous utilisons les mêmes notations que dans le chapitre précédent pour les observations (voir Partie.3.2.2). $O_{[s,t]}$ représente le nombre de sauts d'amplitude -1 du processus $(I_t)_{t \ge 0}$ sur l'intervalle $[s, t]$. Le processus joint $(E_t, I_t, O_{[0,t]})_{t \ge 0}$ est caractérisé par les probabilités de transition suivantes : pour $t > 0$

$$p_{(e,i),(e',i',o)}(t) := \mathbb{P}((E_t, I_t) = (e', i'), O_{[0,t]} = o|(E_0, I_0) = (e, i)).$$

Ces probabilités satisfont une équation de Kolmogorov présentée dans le lemme suivant.

**Lemme 4.1.** Pour $t > 0$ et $e, i, e', i', o \in \mathbb{N}$, on a

$$\frac{d}{dt}p_{(e,i),(e',i',o)}(t) = (\lambda i + \nu)p_{(e+1,i),(e',i',o)}(t) + \alpha e p_{(e-1,i+1),(e',i',o)}(t) + \mu i p_{(e,i-1),(e',i',o-1)}(t) \\ -((\lambda + \mu)i + \alpha e + \nu)p_{(e,i),(e',i',o)}(t), \qquad (4.6)$$

pour $i \le i' + o$, avec $p_{(e,i),(e',i',o)}(t) = 0$ quand $i > i' + o$ et avec la condition initiale $p_{(e,i),(e',i',o)}(0) = \delta_{i=i'}\delta_{o=0}$.

*Preuve.* Nous utilisons la caractérisation infinitésimale de la distribution du processus joint.

Pour $h > 0$ et $i \le i' + o$, on a

$$\begin{align*}
p_{(e,i),(e',i',o)}(t+h) \\
&= \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[0,t+h]} = o|(E_0, I_0) = (e, i)) \\
&= \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o, (E_h, I_h) = (e+1, i), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) \\
&\quad + \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o, (E_h, I_h) = (e-1, i+1), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) \\
&\quad + \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o-1, (E_h, I_h) = (e, i-1), O_{[0,h]} = 1|(E_0, I_0) = (e, i)) \\
&\quad + \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o, (E_h, I_h) = (e, i), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) + 0(h) \\
&= \mathbb{P}((E_h, I_h) = (e+1, i), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o | \\
&\quad (E_h, I_h) = (e+1, i)) + \mathbb{P}((E_h, I_h) = (e-1, i+1), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) \mathbb{P}((E_{t+h}, \\
&\quad I_{t+h}) = (e', i'), O_{[h,t+h]} = o | (E_h, I_h) = (e-1, i+1)) + \mathbb{P}((E_h, I_h) = (e, i-1), O_{[0,h]} = 1 | \\
&\quad (E_0, I_0) = (e, i)) \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o-1 | (E_h, I_h) = (e, i-1)) + \mathbb{P}((E_h, \\
&\quad I_h) = (e, i), O_{[0,h]} = 0|(E_0, I_0) = (e, i)) \mathbb{P}((E_{t+h}, I_{t+h}) = (e', i'), O_{[h,t+h]} = o | (E_h, I_h) = (e, i)) \\
&\quad + 0(h)
\end{align*}$$

$(E_t, I_t)_{t \ge 0}$ est un processus de Markov homogène, ce qui donne :

$$p_{(e,i),(e',i',o)}(t+h) = (\lambda i + \nu)hp_{(e+1,i),(e',i',o)}(t) + \alpha e hp_{(e-1,i+1),(e',i',o)}(t) + \mu i hp_{(e,i-1),(e',i',o-1)}(t) \\ + (1 - ((\lambda + \mu)i + \alpha e + \nu)h)p_{(e,i),(e',i',o)}(t) + 0(h).$$

Le résultat souhaité est donc obtenu. □
Les équations 4.2 et 4.6 ne possèdent pas de solution analytique, ce qui nous empêche, d'une part, d'exprimer la distribution marginale des observations en fonction des paramètres λ, μ, α et ν, et d'autre part, d'obtenir une expression explicite des probabilités de transition par la méthode des fonctions génératrices des moments comme dans le chapitre précédent. En conséquence, une méthode d'estimation directe par maximum de vraisemblance n'est pas réalisable avec ce schéma d'observations. Bien que la fonction génératrice des moments du processus $(E_t, I_t)_{t \ge 0}$ puisse être calculée (voir Annexe D), il n'est pas possible d'identifier la loi de $(E_t, I_t)_{t \ge 0}$. Pour contourner ces difficultés, nous allons adopter une stratégie d'estimation présentée dans la sous-section suivante.

### 4.3.2 Stratégie d'estimation

Comme mentionné dans le chapitre précédent, notre problème d'estimation relève toujours de la classe des problèmes d'information cachée. La suite d'observations $(O_n := O_{[n\Delta t, (n+1)\Delta t]})_{n \in \mathbb{N}}$ est la même, mais la suite d'états cachés est cette fois-ci bivariée et couplée $(E_{n\Delta t}, I_{n\Delta t})_{n \in \mathbb{N}}$. Nous avons donc affaire à un modèle de type multi-chaîne cachée. Nous définissons pour la suite le processus $N_t := \sum_{k=0}^{[t]} O_k$ décrivant le nombre cumulé de nouveaux cas isolés de 0 à l'instant $t$. Notre objectif est d'estimer le quadruplet $(\lambda, \mu, \alpha, \nu)$. La première étape de notre méthode d'estimation consiste à aborder le cas où les variables $E_t$ et $I_t$ sont observées à chaque instant $t$. Nous étudions le comportement moyen du processus joint $(E_t, I_t, N_t)$, nous sommes en mesure d'exprimer explicitement leurs moments en fonction des paramètres $\lambda, \mu, \alpha$ et $\nu$. Sous la condition de stabilité mentionnée dans le corollaire 4.1, le système d'espérances obtenu est inversible, ce qui nous permet de déduire une expression explicite des paramètres en termes des moments du processus $(E_t, I_t, N_t)$ (4.3.3).

Nous pouvons alors estimer les paramètres par plug-in à partir de tout estimateur des moments.

Comme précédemment 3.3.2, la deuxième phase prend en compte le cadre des données manquantes en considérant la chaîne de Markov $(E_{n\Delta t}, I_{n\Delta t}, O_n)_{n \in \mathbb{N}}$ comme un modèle de Markov multi-chaîne cachée (en anglais : Hidden Multi-Chain Markov Model (HMCMM)). Nous avons réécrit les caractéristiques du modèle en tenant compte notre contexte particulier et adapté l'algorithme de Baum Welsh permettant d'estimer la probabilité de transition de la chaîne de Markov cachée couplée à partir de la suite d'observations (4.3.4). Nous estimons ensuite les moments du processus à l'aide de simulations de Monte Carlo et utilisons les formules analytiques obtenues lors de la première étape pour obtenir des estimations des paramètres du modèle exposé-infecté.

### 4.3.3 Expression analytique des paramètres

Le théorème suivant présente les expressions analytiques des paramètres $(\lambda, \mu, \alpha, \nu)$ en termes des moments du processus $(E_t, I_t, N_t)_{t \ge 0}$.

**Théorème 4.3.** Si $\lambda < \mu$, alors les paramètres $\lambda, \mu, \alpha$ et $\nu$ sont donnés par les formules
suivantes

$$
\begin{align}
\lambda &= \frac{N^* \left( \frac{R^*}{E^* I^*} - 1 \right) (E^* + I^*)}{I^* \left( 1 + \left( \frac{R^*}{E^* I^*} - 1 \right) (E^* + I^*) \right)}, \nonumber \\
\mu &= \frac{N^*}{I^*}, \tag{4.7} \\
\alpha &= \frac{N^*}{E^*}, \nonumber \\
\nu &= I^* (\mu - \lambda), \nonumber
\end{align}
$$

où

$$
N^* = \lim_{t \to +\infty} \frac{\mathbb{E}_{(e_0, i_0, n_0)}[N_t]}{t}, \quad E^* = \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[E_t], \quad I^* = \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[I_t],
$$

et

$$
R^* = \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[E_t I_t],
$$

pour $e_0, i_0, n_0 \in \mathbb{N}$.

Preuve. Soit $(E_t, I_t, N_t)_{t \ge 0}$ un processus markovien à valeurs dans $\tilde{S} = \mathbb{N} \times \mathbb{N} \times \mathbb{N}$, de générateur $\tilde{Q}$ donné par : pour tout $e, i, n \ge 0$ :

$$
\begin{align*}
\tilde{Q}_{(e,i,n),(e+1,i,n)} &= \lambda i + \nu, \\
\tilde{Q}_{(e,i,n),(e-1,i+1,n)} &= \alpha e, \\
\tilde{Q}_{(e,i,n),(e,i-1,n+1)} &= \mu i, \\
\tilde{Q}_{(e,i,n),(e,i,n)} &= -((\lambda + \mu)i + \alpha e + \nu).
\end{align*}
$$

Soit $(e_0, i_0, n_0) \in \tilde{S}$, pour toute fonction $f$ positive ou bornée sur $\tilde{S}$ et pour tout $e, i, n \ge 0$,
Nous avons :

$$
\mathbb{E}_{(e_0, i_0, n_0)}[f(E_t, I_t, N_t)] = \sum_{(e,i,n) \in \tilde{S}} p_{(e_0, i_0, n_0)}(e,i,n)(t)f(e, i, n).
$$

En utilisant l'équation de Kolmogorov progressive, on obtient :

$$
\begin{align*}
\frac{d}{dt} \mathbb{E}_{(e_0, i_0, n_0)}[f(E_t, I_t, N_t)]
&= \sum_{(e,i,n) \in \tilde{S}} p'_{(e_0, i_0, n_0)}(e,i,n)(t) f(e, i, n) \\
&= \sum_{(e,i,n) \in \tilde{S}} \left( \sum_{(e',i',n') \in \tilde{S}} p_{(e_0, i_0, n_0)}(e',i',n') \tilde{Q}_{(e',i',n')}(e,i,n) \right) f(e, i, n) \\
&= \sum_{(e',i',n') \in \tilde{S}} p_{(e_0, i_0, n_0)}(e',i',n') \left( \sum_{(e,i,n) \in \tilde{S}} \tilde{Q}_{(e',i',n')}(e,i,n) f(e, i, n) \right) \\
&= \sum_{(e',i',n') \in \tilde{S}} p_{(e_0, i_0, n_0)}(e',i',n') \left( (\lambda i' + \nu) (f(e' + 1, i', n') - f(e', i', n')) \right. \\
&\qquad \left. + \alpha e' (f(e' - 1, i' + 1, n') - f(e', i', n')) \right. \\
&\qquad \left. + \mu i' (f(e', i' - 1, n' + 1) - f(e', i', n')) \right) \\
&= \mathbb{E}_{(e_0, i_0, n_0)}\left[ (\lambda I_t + \nu) (f(E_t + 1, I_t, N_t)) - f(E_t, I_t, N_t) \right] \\
&\quad + \alpha E_t \left( f(E_t - 1, I_t + 1, N_t) - f(E_t, I_t, N_t) \right) + \mu I_t \left( f(E_t, I_t - 1, N_t + 1) \right. \\
&\qquad \left. - f(E_t, I_t, N_t) \right).
\end{align*}
$$
En considérant des fonctions $f(E_t, I_t, N_t)$, telles que $E_t, I_t, N_t, E_t^2, E_t I_t$ et $I_t^2$, on peut établir le système d'équations différentielles suivant :

$$
\left\{
\begin{aligned}
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[E_t] &= -\alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t] + \lambda \mathbb{E}_{(e_0,i_0,n_0)}[I_t] + \nu, \\
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[I_t] &= \alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t] - \mu \mathbb{E}_{(e_0,i_0,n_0)}[I_t], \\
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[N_t] &= \mu \mathbb{E}_{(e_0,i_0,n_0)}[I_t], \\
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[E_t^2] &= (2\nu + \alpha) \mathbb{E}_{(e_0,i_0,n_0)}[E_t] + \lambda \mathbb{E}_{(e_0,i_0,n_0)}[I_t] - 2\alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t^2] \\
&\quad + 2\lambda \mathbb{E}_{(e_0,i_0,n_0)}[E_t I_t] + \nu, \\
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[E_t I_t] &= -\alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t] + \nu \mathbb{E}_{(e_0,i_0,n_0)}[I_t] + \alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t^2] \\
&\quad - (\mu + \alpha) \mathbb{E}_{(e_0,i_0,n_0)}[E_t I_t] + \lambda \mathbb{E}_{(e_0,i_0,n_0)}[I_t^2], \\
\frac{d}{dt} \mathbb{E}_{(e_0,i_0,n_0)}[I_t^2] &= \alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t] + \mu \mathbb{E}_{(e_0,i_0,n_0)}[I_t] + 2\alpha \mathbb{E}_{(e_0,i_0,n_0)}[E_t I_t] - 2\mu \mathbb{E}_{(e_0,i_0,n_0)}[I_t^2].
\end{aligned}
\right.
$$

En résolvant le système (4.8), nous pouvons obtenir les expressions des espérances suivantes :

$$
\mathbb{E}_{(e_0,i_0,n_0)}[E_t] = c_1 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_2 e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\mu\nu}{\alpha(\mu-\lambda)},
$$

$$
\mathbb{E}_{(e_0,i_0,n_0)}[I_t] = c_3 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_4 e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\nu}{\mu-\lambda},
$$

$$
\mathbb{E}_{(e_0,i_0,n_0)}[N_t] = c_5 + c_6 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_7 e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + \frac{\mu\nu}{\mu-\lambda}t,
$$

$$
\mathbb{E}_{(e_0,i_0,n_0)}[E_t I_t] = c_8 e^{-(\mu+\alpha)t} + c_9 e^{-\frac{\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} + c_{10} e^{-\frac{\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda}}{2}t} \\
c_{11} e^{-(\mu+\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda})t} + c_{12} e^{-(\mu-\alpha+\sqrt{(\mu-\alpha)^2+4\alpha\lambda})t} + \frac{\mu\nu((\mu+\alpha)\nu + \alpha\lambda)}{\alpha(\mu-\lambda)^2(\mu+\alpha)},
$$

où $c_i$ pour $i = 1, \dots, 12$ sont des constantes réelles.

Ainsi, si $\lambda < \mu$, alors

$$
E^* := \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[E_t] = \frac{\mu\nu}{\alpha(\mu - \lambda)},
$$

$$
I^* := \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[I_t] = \frac{\nu}{\mu - \lambda},
$$

$$
N^* := \lim_{t \to +\infty} \frac{\mathbb{E}_{(e_0, i_0, n_0)}[N_t]}{t} = \frac{\mu\nu}{\mu - \lambda}, \qquad (4.8)
$$

$$
R^* := \lim_{t \to +\infty} \mathbb{E}_{(e_0, i_0, n_0)}[E_t I_t] = \frac{\mu\nu((\mu + \alpha)\nu + \alpha\lambda)}{\alpha(\mu - \lambda)^2 (\mu + \alpha)}.
$$

Ces quatre équations (4.8) forment un système inversible en fonction des paramètres $(\lambda, \mu, \alpha, \nu)$, car $\lambda \neq \mu$.
$$
\begin{align*}
\lambda &= \frac{N^* \left( \frac{R^*}{E^* I^*} - 1 \right) (E^* + I^*)}{I^* \left( 1 + \left( \frac{R^*}{E^* I^*} - 1 \right) (E^* + I^*) \right)}, \\
\mu &= \frac{N^*}{I^*}, \\
\alpha &= \frac{N^*}{E^*}, \\
\nu &= I^* (\mu - \lambda).
\end{align*}
$$

Ainsi, nous pouvons déduire des estimateurs (de type moment) des paramètres λ, μ, ν et α.

Le processus $N_t$ est observé mais les valeurs $(E_n, I_n)_{n \in \mathbb{N}} := (E_{n\Delta t}, I_{n\Delta t})_{n \in \mathbb{N}}$ sont ca-
chées dans notre contexte, nous devons donc recourir à une méthode adéquate afin d'esti-
mer les valeurs $E^*, I^*$ et $R^*$. C'est la deuxième étape de notre méthodologie, pour laquelle
nous utiliserons un modèle de Markov caché approprié.

4.3.4 Modèle de Markov Multi-Chaîne Cachée

Le modèle de Markov Multi-Chaîne Cachée (en anglais : Hidden Multi-Chain Markov
Model (HMCMM)) est un outil statistique très puissant pour le traitement de données
séquentielles multivariées. Il est considéré comme une extension du modèle de Markov
caché classique en introduisant des probabilités conditionnelles couplées entre les suites
d'états de plusieurs HMM. Cette extension permet de capturer les dépendances complexes
entre les différentes suites d'états, ce qui en fait un choix très utile pour de nombreuses
applications [73, 74, 75, 76].

Le but de cette section est d'adapter les algorithmes du HMM afin d'estimer les pro-
babilités de transition de la chaîne cachée (E_n, I_n)_{n \in \mathbb{N}} à partir des observations (O_n)_{n \in \mathbb{N}}.
Ensuite, nous utiliserons ces estimations pour obtenir les valeurs E*, I* et R* à l'aide d'une
méthode de Monte Carlo. Dans notre modèle de Markov caché exposé-infecté (Fig.4.6),
chaque chaîne représente l'état d'une population (exposée, infectée ou isolée) et où les
chaînes sont couplées entre elles. À chaque instant n, l'état caché des infectés dépend de
l'état de sa propre chaîne et de la chaîne de couplage des exposés à l'instant n - 1. De
même, l'état caché des exposés dépend de l'état de sa propre chaîne et de la chaîne de
couplage des infectés à l'instant n - 1. Ce couplage entre les chaînes reflète l'interaction
entre les populations exposée et infectée dans la propagation de la maladie. Cependant,
la valeur observée O_n à chaque instant dépend non seulement de l'état actuel I_n, mais
aussi de l'état précédent I_{n-1}, ce qui représente un schéma d'émission non standard.
<figure><img src="image_38.png" /><figcaption>FIGURE 4.6 – Modèle de Markov multi chaîne cachée (exposé-infecté).</figcaption></figure>

Nous pouvons ramener notre HMCMM à un HMM avec un schéma d'observations standard en introduisant la chaîne tridimensionnelle ($Z_n)_{n \in \mathbb{N}^* = (E_{n-1}, I_{n-1}, I_n)_{n \in \mathbb{N}^*}$, qui représente une chaîne de Markov cachée à valeurs dans l'espace d'état $\mathbb{N}^3$. Cette transformation permet de démontrer que la suite $(Z_n, O_n)$ est un modèle de Markov caché standard avec les caractéristiques suivantes.

**Lemme 4.2.** Le processus $(Z_n, O_n)_{n \in \mathbb{N}^*}$ est un modèle de Markov caché dont les caractéristiques sont données par le triple $M = (Q, \psi, \rho)$, où

i) La matrice de probabilité de transition $Q$ de la chaîne cachée $(Z_n)$ est

$$Q_{(e,i,j),(e',i',j')} = \frac{p_{(e',i'),(\cdot,j')}p_{(e,i),(e',i')}}{p_{(e,i),(\cdot,j)}} \delta_{i'=j} \quad (e, i, j), (e', i', j') \in \mathbb{N}^3,$$

où $p_{(e,i),(e',i')} = p_{(e,i),(e',i')}(\Delta t)$ est la matrice de transition de la chaîne $(E_n, I_n)_{n \in \mathbb{N}}$ et $p_{(e,i),(\cdot,j)} = \sum_k p_{(e,i),(k,j)}$.

ii) La probabilité d'émission $\psi$ du processus $O$ sachant le processus $Z$ est

$$\psi_{(e,i,j)}(o) = \mathbb{P}(O_n = o | Z_n = (e, i, j)) = \frac{p_{(e,i),(\cdot,j,o)}}{p_{(e,i),(\cdot,j)}} \quad e, i, j, o \in \mathbb{N},$$

où $p_{(e,i),(\cdot,j,o)} = \sum_k p_{(e,i),(k,j,o)}(\Delta t)$ est caractérisée par l'équation **4.6**.

iii) La loi initiale $\rho$ du processus d'état $Z$

$$\rho_{e,i,j} = \mathbb{P}(Z_1 = (e, i, j)) = p_{(e,i),(\cdot,j)}\pi_{(e,i)} \quad e, i, j \in \mathbb{N},$$

où $\pi$ est la loi initiale de la chaîne $(E_n, I_n)_{n \in \mathbb{N}}$.

*Preuve.* i) La probabilité de transition $Q$ : soient $(e, i, j), (e', i', j') \in \mathbb{N}^3$, nous avons

$$\begin{aligned} Q_{(e,i,j),(e',i',j')} &= \mathbb{P}(Z_{n+1} = (e', i', j') | Z_n = (e, i, j)) \\ &= \frac{\mathbb{P}(E_n = e', I_n = i', I_{n+1} = j', E_{n-1} = e, I_{n-1} = i, I_n = j)}{\mathbb{P}(E_{n-1} = e, I_{n-1} = i, I_n = j)}. \end{aligned}$$
En utilisant la propriété de Markov de la chaîne $(E_n, I_n)$, nous obtenons

$$
\begin{align*}
Q_{(e,i,j),(e',i',j')} &= \frac{\mathbb{P}(I_{n+1} = j' | E_n = e', I_n = i') \mathbb{P}(E_n = e', I_n = i' | E_{n-1} = e, I_{n-1} = i)}{\mathbb{P}(I_n = j | E_{n-1} = e, I_{n-1} = i)} \delta_{i'=j} \\
&= \frac{p(e',i'),(\cdot,j')p(e,i),(e',i')}{p(e,i),(\cdot,j)} \delta_{i'=j}.
\end{align*}
$$

ii) La probabilité d'émission $\psi$ : c'est évident par la formule de la probabilité condition-nelle.

iii) La loi initiale $\rho$ : soit $e, i, j \in \mathbb{N}$

$$
\begin{align*}
\rho_{e,i,j} &= \mathbb{P}(Z_1 = (e, i, j)) \\
&= \mathbb{P}(I_1 = j | (E_0, I_0) = (e, i)) \mathbb{P}((E_0, I_0) = (e, i)) \\
&= p_{(e,i),(\cdot,j)} \pi_{(e,i)}.
\end{align*}
$$

□

Dans la suite, nous visons à estimer les paramètres du modèle de Markov caché $(Z_n, O_n)$, à savoir la matrice de transition $p$ de la chaîne $(E_n, I_n)$, en maximisant la fonction de vraisemblance $\mathbb{P}(O_1 = o_1, \dots, O_T = o_T | M)$ à l'aide de l'algorithme de Baum-Welch ou l'algorithme forward-backward. Comme dans le chapitre précédent, nous avons réécrit les fonctions principales de l'algorithme EM, en adaptant les formules récursives à la structure particulière de notre modèle pour obtenir des estimateurs explicites des probabilités de transition $p_{(e,i),(e',i')}$.

Pour implémenter cet algorithme de maximisation de la vraisemblance des observations,
nous devons d'abord définir les probabilités Forward et Backward.

La probabilité Forward peut être calculée par la récursion suivante : pour tout $e, i, j \in \mathbb{N}$ et $2 \le t \le T$

$$
\left\{
\begin{array}{l}
\alpha_{(e,i,j)}(1) = p_{(e,i),(\cdot,j)}\pi_{(e,i)}\psi_{(e,i,j)}(o_1), \\
\alpha_{(e,i,j)}(t) = \psi_{(e,i,j)}(o_t) \sum_{(e',i') \in \mathbb{N}^2} \frac{p_{(e,i),(\cdot,j)}p_{(e',i'),(e,i)}}{p_{(e',i'),(\cdot,i)}}\alpha_{(e',i',i)}(t-1).
\end{array}
\right.
$$

La probabilité Bakward peut être évaluée par la formule récursive : $e, i, j \in \mathbb{N}$ et $1 \le t \le T - 1$

$$
\left\{
\begin{array}{l}
\beta_{(e,i,j)}(T) = 1, \\
\beta_{(e,i,j)}(t) = \sum_{(e',j') \in \mathbb{N}^2} \frac{p_{(e',j),(\cdot,j')}p_{(e,i),(e',j)}}{p_{(e,i),(\cdot,j)}} \psi_{(e',j,j')} (o_{t+1}) \beta_{(e',j,j')}(t+1).
\end{array}
\right.
$$

En utilisant les probabilités Forward et Backward, la vraisemblance peut être calculée en additionnant le produit de ces deux probabilités pour tous les états de la chaîne cachée, comme indiqué par la formule suivante :

$$
\mathbb{P}(O_1 = o_1, \dots, O_T = o_T | M) = \sum_{(e,i,j) \in \mathbb{N}^3} \alpha_{(e,i,j)}(t) \beta_{(e,i,j)}(t).
$$

Le théorème suivant présente le schéma itératif de l'algorithme de Baum-Welch adapté.
**Théorème 4.4.** Étant donné les paramètres $M^n = (Q^n, \psi^n, \rho^n)$, les estimations du maximum de vraisemblance de $M^{n+1} = (Q^{n+1}, \psi^{n+1}, \rho^{n+1})$ étant donnée les observations $(o_1, \dots, o_T)$ sont données, pour $e, i, j, o \in \mathbb{N}$, par

$$
\begin{align*}
Q_{(e,i,j),(e',i',j')}^{n+1} &= \frac{\sum_{t=1}^{T-1} \xi_{(e,i,j),(e',i',j')}^{n}(t)}{\sum_{t=1}^{T} \gamma_{(e,i,j)}^{n}(t)} \delta_{i'=j}, \\
\psi_{(e,i,j)}^{n+1}(o) &= \frac{\sum_{t=1}^{T} \mathbb{1}_{o_t(o)} \gamma_{(e,i,j)}^{n}(t)}{\sum_{t=1}^{T} \gamma_{(e,i,j)}^{n}(t)}, \\
\rho_{e,i,j}^{n+1} &= \gamma_{(e,i,j)}^{n}(1), \\
p_{(e,i),(e',i')}^{n+1} &= \frac{\sum_{j' \in \mathbb{N}} Q_{(e,i,j),(e',j,j')}^{n+1} \rho_{e,i,j}^{n+1}}{\sum_{j \in \mathbb{N}} \rho_{e,i,j}^{n+1}},
\end{align*}
$$

où

$$
\begin{align*}
\xi_{(e,i,j),(e',i',j')}^{n}(t) &= \mathbb{P}(Z_t = (e, i, j), Z_{t+1} = (e', i', j') | O_{1:T} = o_{1:T}, M^n) \\
&= \frac{\alpha_{(e,i,j)}^{n}(t) p_{(e',i'),(\cdot,j')}^{n} p_{(e,i),(e',i')} \psi_{(e',i',j')}^{n}(o_{t+1}) \beta_{(e',i',j')}^{n}(t+1)}{p_{(e,i),(\cdot,j)}^{n} \sum_{(e,i,j) \in \mathbb{N}^3} \alpha_{(e,i,j)}^{n}(t) \beta_{(e,i,j)}^{n}(t)} \delta_{i'=j},
\end{align*}
$$

et

$$
\begin{align*}
\gamma_{(e,i,j)}^{n}(t) &= \mathbb{P}(Z_t = (e, i, j) | O_{1:T} = o_{1:T}, M^n) \\
&= \frac{\alpha_{(e,i,j)}^{n}(t) \beta_{(e,i,j)}^{n}(t)}{\sum_{(e,i,j) \in \mathbb{N}^3} \alpha_{(e,i,j)}^{n}(t) \beta_{(e,i,j)}^{n}(t)},
\end{align*}
$$

avec $\alpha_{(e,i,j)}^n(t), \beta_{(e,i,j)}^n(t)$ sont les probabilités Forward et Backward à l'itération n.

Preuve. Notons par $\mathcal{Z} = \mathbb{N}^3$ l'espace d'état de la chaîne $(Z_n)_{n \in \mathbb{N}^*}$ et $z = (e, i, j) \in \mathcal{Z}$. On cherche à maximiser la vraisemblance d'observations $\mathbb{P}(O_{1:T} = o_{1:T} | M^{n+1})$. En appliquant l'algorithme EM [77] à la maximisation de cette probabilité, on est amené à maximiser la fonction auxiliaire suivante :

$$
\Gamma_{O_{1:T}}(M^n, M^{n+1}) = \sum_{z_{1:T} \in \mathcal{Z}^T} \log \mathbb{P}(O_{1:T} = o_{1:T}, Z_{1:T} = z_{1:T} | M^{n+1}) \mathbb{P}(Z_{1:T} = z_{1:T} | O_{1:T}, M^n).
$$

Sachant que

$$
\mathbb{P}(O_{1:T} = o_{1:T}, Z_{1:T} = z_{1:T} | M^{n+1}) = \rho_{z_1}^{n+1} \left( \prod_{t=1}^{T-1} Q_{z_t, z_{t+1}}^{n+1} \right) \left( \prod_{t=1}^{T} \psi_{z_t}^{n+1}(o_t) \right),
$$

la fonction $\Gamma_{O_{1:T}}$ se ré-écrit

$$
\begin{align*}
\Gamma_{O_{1:T}}(M^n, M^{n+1}) &= \sum_{z_{1:T} \in \mathcal{Z}^T} \log \rho_{z_1}^{n+1} \mathbb{P}(Z_{1:T} = z_{1:T} | O_{1:T}, M^n) \\
&\quad + \sum_{z_{1:T} \in \mathcal{Z}^T} \left( \sum_{t=1}^{T-1} \log Q_{z_t, z_{t+1}}^{n+1} \right) \mathbb{P}(Z_{1:T} = z_{1:T} | O_{1:T}, M^n) \\
&\quad + \sum_{z_{1:T} \in \mathcal{Z}^T} \left( \sum_{t=1}^{T} \log \psi_{z_t}^{n+1}(o_t) \right) \mathbb{P}(Z_{1:T} = z_{1:T} | O_{1:T}, M^n) \\
&= \Gamma_{O_{1:T}}^{\rho}(M^n, M^{n+1}) + \Gamma_{O_{1:T}}^{Q}(M^n, M^{n+1}) + \Gamma_{O_{1:T}}^{\psi}(M^n, M^{n+1}).
\end{align*}
$$
On remarque alors que la fonction $\Gamma_{O_{1:T}}$ peut se décomposer en trois fonctions ayant des paramètres distincts et indépendants. Par conséquent, nous pouvons optimiser analytiquement ces trois fonctions à l'aide de multiplicateurs de Lagrange. Soient $\mathcal{L}(\rho, \theta)$, $\mathcal{L}(Q, \theta)$ et $\mathcal{L}(\psi, \theta)$ les Lagrangiens associés respectivement à $\Gamma_{O_{1:T}}^{\rho}$, $\Gamma_{O_{1:T}}^{Q}$ et $\Gamma_{O_{1:T}}^{\psi}$.

$$ \mathcal{L}(\rho, \theta) = \sum_{(e,i,j)} \log \rho_{(e,i,j)}^{n+1} \mathbb{P}(Z_1 = (e, i, j) | O_{1:T}, M^n) + \theta \left( \sum_{(e',i',j')} \rho_{(e',i',j')}^{n+1} - 1 \right). $$

$$ \mathcal{L}(Q, \theta) = \sum_{(e,i,j)} \sum_{(e',i',j')} \left( \sum_{t=1}^{T-1} \log Q_{(e,i,j),(e',i',j')}^{n+1} \right) \mathbb{P}(Z_t = (e, i, j), Z_{t+1} = (e', i', j') | O_{1:T}, M^n) \\ + \sum_{(e,i,j)} \theta_{(e,i,j)} \left( \sum_{(e',i',j')} Q_{(e,i,j),(e',i',j')}^{n+1} - 1 \right). $$

$$ \mathcal{L}(\psi, \theta) = \sum_{(e,i,j)} \left( \sum_{t=1}^{T} \log \psi_{(e,i,j)}^{n+1}(o) \right) \mathbb{P}(Z_t = (e, i, j) | O_{1:T}, M^n) \mathbb{1}_{o_t(o)} \\ + \sum_{(e,i,j)} \theta_{(e,i,j)} \left( \sum_{o} \psi_{(e,i,j)}^{n+1}(o) - 1 \right). $$

En mettant les dérivées premières des Lagrangiens à zéro, et en procédant de manière similaire à la preuve du théorème 3.3, on peut montrer que $\Gamma_{O_{1:T}}$ est maximisé lorsque

$$ Q_{(e,i,j),(e',i',j')}^{n+1} = \frac{\sum_{t=1}^{T-1} \xi_{(e,i,j),(e',i',j')}^{n}(t)}{\sum_{t=1}^{T} \gamma_{(e,i,j)}^{n}(t)} \delta_{i'=j}, $$

$$ \psi_{(e,i,j)}^{n+1}(o) = \frac{\sum_{t=1}^{T} \mathbb{1}_{o_t(o)} \gamma_{(e,i,j)}^{n}(t)}{\sum_{t=1}^{T} \gamma_{(e,i,j)}^{n}(t)}, $$

$$ \rho_{e,i,j}^{n+1} = \gamma_{(e,i,j)}^{n}(1). $$

D'autre part, en utilisant le lemme 4.2, on peut facilement déduire que

$$ \sum_{j' \in \mathbb{N}} Q_{(e,i,j),(e',j,j')}^{n+1} \rho_{e,i,j}^{n+1} = p_{(e,i),(e',i')}^{n+1} \pi_{(e,i)}^{n+1}. $$

Or

$$ \sum_{j \in \mathbb{N}} \rho_{e,i,j}^{n+1} = \pi_{(e,i)}^{n+1}. $$

Par conséquence :

$$ p_{(e,i),(e',i')}^{n+1} = \frac{\sum_{j' \in \mathbb{N}} Q_{(e,i,j),(e',j,j')}^{n+1} \rho_{e,i,j}^{n+1}}{\sum_{j \in \mathbb{N}} \rho_{e,i,j}^{n+1}}. $$

□

À chaque itération, l'algorithm de Baum-Welch met à jour les paramètres du modèle en utilisant les formules données dans le théorème (4.4) et réévalue la log-vraisemblance des observations, jusqu'à ce qu'une convergence vers les valeurs optimales soit atteinte. L'obtention de la matrice de transition $\hat{p}$ optimale permet de simuler la chaîne de Markov $(E_n, I_n)$ avec une taille suffisamment grande de $n$ pour récupérer les moments $E^*$, $I^*$ et $R^*$. Finalement, le théorème (4.3) peut être utilisé pour estimer les paramètres du modèle exposé-infecté $\lambda, \mu, \alpha$ et $\nu$.
## 4.4 Simulation numérique

Nous allons étudier les performances numériques des estimateurs $\hat{\lambda}^n, \hat{\mu}^n, \hat{\alpha}^n$ et $\hat{\nu}^n$ sur des données synthétiques. Nous commencerons par examiner le cadre où la chaîne bivariée $(E_n, I_n)$ est observé, puis nous nous tournons vers le cadre HMCMM.

### 4.4.1 Observations complètes

Pour estimer les moments $E^*, I^*$ et $R^*$ de la chaîne de Markov observée $(E_n, I_n)$, nous utilisons la méthode de Monte Carlo qui consiste à simuler un grand nombre de trajectoires de la chaîne à partir desquelles nous pouvons estimer les moments empiriques.

Nous simulons plusieurs trajectoires du processus $(E_t, I_t)$ partant d'un état initial $(E_0, i_0) = (0, 0)$ avec des valeurs de paramètres fixées à $\lambda = 0.05, \mu = 0.2, \alpha = 0.1$ et $\nu = 0.015$, sur différents horizons du temps $H$. Pour chaque trajectoire, avec un pas de temps fixé à $\Delta t = 1j$, nous obtenons la chaîne discrète $(E_n, I_n)$ et les observations $O_n$ qui décrivent le nombre de sauts négatifs du processus marginal $I_t$ sur l'intervalle $]n\Delta t, (n + 1)\Delta t]$. Nous utilisons ensuite les moyennes empiriques de la chaîne pour estimer les moments, puis nous obtenons des estimations des paramètres en les injectant dans la formule (4.4). Cette méthode nous permet d'obtenir une distribution des estimateurs (Fig.4.7), à partir de laquelle nous pouvons calculer des intervalles de confiance.

Afin de réduire le temps de calcul des simulations de Monte Carlo, nous utilisons la fonction "dopar" du package "doparallel" [78] pour paralléliser nos calculs. Cette méthode permet de diviser le travail en plusieurs tâches indépendantes qui peuvent être exécutées simultanément, ce qui a permis de réduire considérablement le temps nécessaire pour effectuer les simulations.

Le tableau (4.1) présente les valeurs estimées des paramètres ainsi que leurs intervalles de confiance pour différents horizons $H$ (nombre des observations), pour 10000 trajectoires simulées.

<table><caption>TABLE 4.1 – Estimations des paramètres et leurs intervalles de confiance sur 10000 trajectoires du processus exposé-infecté.</caption><thead><tr><th colspan="3">Nombre des observations n = 1000</th></tr><tr><th>Les paramètres</th><th>Les estimations</th><th>Intervalle de confiance (95%)</th></tr></thead><tbody><tr><td>λ = 0.05</td><td>0.049</td><td>[0.047, 0.052]</td></tr><tr><td>μ = 0.2</td><td>0.1999</td><td>[0.1996, 0.2019]</td></tr><tr><td>α = 0.1</td><td>0.10</td><td>[0.0996, 0.1005]</td></tr><tr><td>ν = 0.015</td><td>0.0147</td><td>[0.0145, 0.0149]</td></tr></tbody></table><table><thead><tr><th colspan="3">Nombre des observations n = 10000</th></tr><tr><th>Les paramètres</th><th>Les estimations</th><th>Intervalle de confiance (95%)</th></tr></thead><tbody><tr><td>λ = 0.05</td><td>0.050</td><td>[0.049, 0.0506]</td></tr><tr><td>μ = 0.2</td><td>0.20</td><td>[0.1998, 0.2001]</td></tr><tr><td>α = 0.1</td><td>0.100</td><td>[0.0998, 0.1002]</td></tr><tr><td>ν = 0.015</td><td>0.015</td><td>[0.0149, 0.01506]</td></tr></tbody></table>
<figure><img src="image_39.png" /><figcaption>FIGURE 4.7 – Distributions des 10000 estimateurs simulés sur une longueur *H* = 10000.</figcaption></figure>

Il convient de noter que la précision des estimations dépend non seulement du nombre
de trajectoires simulées, mais aussi de la longueur de chaque trajectoire. En effet, plus
le nombre de trajectoires simulées est élevé, plus la précision de l'estimation sera élevée
et l'intervalle de confiance sera réduit, ce qui signifie que l'estimation des paramètres
sera plus précise et plus fiable. Cependant, il est important de noter que la longueur de
chaque trajectoire peut également affecter la précision de l'estimation, en particulier si
la longueur est trop courte, ce qui peut entraîner une variance élevée des estimations.
Il est donc important de trouver un compromis entre le nombre de trajectoires et la
longueur de chaque trajectoire pour obtenir des estimations précises des paramètres avec
des intervalles de confiance resserrés.

4.4.2 Observations du cumul des nouveaux cas déclarés

Nous allons maintenant nous concentrer sur le cadre HMCMM, plus adapté à notre
contexte d'étude, où seuls les nouveaux cas déclarés cumulés ($O_n$) sont observés. Pour
cela, nous avons simulé un processus exposé-infecté défini par son générateur infinitésimal
(4.1), à partir du point initial (0,0), avec des paramètres $\lambda = 0.05$, $\mu = 0.2$, $\alpha = 0.1$
et $\nu = 0.015$, sur une période de temps *H* (Fig.4.8). Nous avons ensuite extrait les
observations correspondant aux cumuls des nouveaux sauts -1 du processus du processus
$I_t$ (les cumuls de nouveaux cas déclarés Fig.4.9), afin d'exécuter notre schéma adaptatif
de Baum-Welch et estimer les paramètres recherchés.
<figure><img src="image_40.png" /><figcaption>FIGURE 4.8 – Trajectoire du processus exposé-infecté ($E_t, I_t$) sur un horizon $H = 10000$ jours.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>FIGURE 4.9 – Les cumuls des nouveaux déclarés $O_n$ sur une période de $\Delta t = 1$ jour.</figcaption></figure>
Pour le reste des simulations, Nous supposons que le nombre d'états cachés des processus $E_t$ et $I_t$ est égal à $N + 1$, et nous tronquons les matrices de transition et d'émission respectivement sur les espace d'états finis $\{0, 1, ..., N\} \times \{0, 1, ..., N\}$ et $\{0, 1, ..., N\} \times \{0, 1, ..., N\} \times \{0, 1, ..., M\}$, où M représente l'état maximal observé. Il est important de noter que l'algorithme Baum-Welch est sensible aux choix des paramètres initiaux, car il peut converger vers des optima locaux plutôt que globaux. Si les paramètres initiaux sont trop éloignés des vraies valeurs des paramètres, cela peut entraîner une non-convergence de l'algorithme. Par conséquent, il est important de choisir des paramètres initiaux qui sont proches et dans le même ordre de grandeur que les vraies valeurs des paramètres autant que possible. Cela peut être fait en utilisant des connaissances préalables sur le système étudié ou en effectuant des simulations initiales pour obtenir des valeurs raisonnables pour les paramètres. Dans notre étude, nous avons décidé de sélectionner les paramètres initiaux dans les plages suivantes : $\lambda^{(0)} \in [0.04, 0.07]$, $\mu^{(0)} \in [0.185, 0.25]$, $\alpha^{(0)} \in [0.09, 0.130]$ et $\nu^{(0)} \in [0.013, 0.02]$ avec $length.out = 15$ (en générant 15 valeurs initiales pour chaque paramètre), cela nous permet d'initialiser les caractéristiques de notre modèle HMCMM de la façon suivante : Pour chaque valeur de quadruplet $(\lambda_i^{(0)}, \mu_i^{(0)}, \alpha_i^{(0)}, \nu_i^{(0)})$ pour $i = 1, ..., 15$, nous effectuons des simulations de Monte Carlo parallélisées sur 10000 trajectoires différentes pour estimer les coefficients des matrices $P^{(0)}$ et $\psi^{(0)}$. Ensuite, nous redéfinissons la dernière colonne de chaque matrice en utilisant les formules suivantes :

$$ 
\begin{aligned}
p_{(e,i)(N,N)} &= 1 - \sum_{e'=0}^{N-1} \sum_{i'=0}^{N-1} p_{(e,i)(e',i')}, \\
\psi_{(e,i,j)}(M) &= 1 - \sum_{o=0}^{M-1} \psi_{(e,i,j)}(o),
\end{aligned}
$$

pour $e, i, j = 0, ..., N$.

Pour les autres paramètres du notre schéma itératif, nous fixons $\Delta t = 1$, $maxIterations = 500$ (nombre maximum d'itérations) et $delta = 10^{-9}$ (critère d'arrêt). Nous exécutons l'algorithme de Baum-Welch adapté pour chaque valeur de $(\lambda_i^{(0)}, \mu_i^{(0)}, \alpha_i^{(0)}, \nu_i^{(0)})$ pour $i = 1, ..., 15$, et nous choisissons finalement la matrice de transition qui maximise la vraisemblance des observations.

$$ \hat{p}^m = \max_{\hat{p}^n} \mathbb{P}(O|M^n). $$

Nous utilisons ensuite cette matrice optimale $\hat{p}^m$ pour simuler la chaîne de Markov $(E_n, I_n)$ avec une taille suffisamment grande de $n$ pour récupérer les moments $E^*$, $I^*$ et $R^*$. Enfin, nous avons appliqué le théorème (4.3) pour estimer les paramètres du modèle exposé-infecté $\lambda, \mu, \alpha$ et $\nu$.

L'impact du choix des paramètres de troncature $N$ et de l'horizon temporel $H$ est examiné dans les tableaux 4.2 et 4.3. Ces paramètres sont importants pour la précision des estimations.
<table><caption>TABLE 4.2 – Estimation des paramètres pour différente valeur de N sur une trajectoire d'horizon H = 1000.</caption><thead><tr><th>Nombre des observations</th><th colspan="3">n = 1000</th></tr><tr><th>Troncature de l'espace d'état caché</th><th>N = 3</th><th>N = 4</th><th>N = 5</th></tr></thead><tbody><tr><th>λ<sup>n</sup></th><td>0.059</td><td>0.057</td><td>0.058</td></tr><tr><th>μ<sup>n</sup></th><td>0.245</td><td>0.248</td><td>0.197</td></tr><tr><th>α<sup>n</sup></th><td>0.120</td><td>0.120</td><td>0.092</td></tr><tr><th>ν<sup>n</sup></th><td>0.0182</td><td>0.0185</td><td>0.014</td></tr></tbody></table>

<table><caption>TABLE 4.3 – Estimation des paramètres pour différente valeur de N sur une trajectoire d'horizon H = 10000.</caption><thead><tr><th>Nombre des observations</th><th colspan="3">n = 10000</th></tr><tr><th>Troncature de l'espace d'état caché</th><th>N = 3</th><th>N = 4</th><th>N = 5</th></tr></thead><tbody><tr><th>λ<sup>n</sup></th><td>0.050</td><td>0.049</td><td>0.05</td></tr><tr><th>μ<sup>n</sup></th><td>0.189</td><td>0.196</td><td>0.209</td></tr><tr><th>α<sup>n</sup></th><td>0.098</td><td>0.101</td><td>0.102</td></tr><tr><th>ν<sup>n</sup></th><td>0.0146</td><td>0.015</td><td>0.015</td></tr></tbody></table>

Ces tableaux montrent les résultats d'estimation des paramètres du modèle pour différentes valeurs de N (le nombre d'états cachés) et de H (l'horizon temporel). Les estimations sont basées sur une trajectoire de simulation de longueur H et pour chaque valeur de N, l'algorithme de Baum-Welch est exécuté pour estimer les paramètres du modèle. On observe que les estimations des paramètres ne varient pas de manière significative avec le choix de N, pour H = 10000, les résultats montrent que pour les trois valeurs de N testées, les estimations des paramètres sont assez stables et proches les unes des autres, excepté pour $\mu^n$ qui a une légère variation. Donc, la troncature d'espace d'état n'a pas un impact significatif sur l'estimation, ceci est dû au fait que la majorité des états observés sont des 0 et 1. Cependant, une grande valeur de N augmente également le temps de calcul nécessaire pour obtenir les résultats.

En revanche, le choix du nombre d'observations $n$ (ou de l'horizon temporel $H$) a un impact significatif sur les estimations obtenues. Les résultats de nos simulations montrent qu'un horizon temporel plus long conduit à des estimations plus précises des paramètres. En effet, un horizon temporel plus long fournit une plus grande quantité d'informations qui peuvent être utilisées pour estimer les paramètres. Par conséquent, l'augmentation du nombre d'observations $n$ améliore d'une part la convergence locale de l'algorithme EM vers la vraie valeur de la matrice de transition et d'autre part l'estimation de $N^*$ qui dépend essentiellement par l'horizon $H$ et du nombre d'observations. Toutefois, il est important de noter que l'augmentation du nombre d'observations entraîne également une augmentation significative du temps d'exécution.

Comme précédemment mentionné dans le chapitre 3 (3.4.1), nos expériences numériques ont démontré que l'identifiabilité des moments est garantie avec notre stratégie d'estimation et notre adaptation du modèle. Cela confirme que notre approche permet
d'estimer les paramètres de manière fiable et de retrouver les caractéristiques désirées. Cependant, il convient de noter que nous n'avons pas pu vérifier explicitement si la matrice optimale obtenue par l'algorithme de Baum-Welch adapté est une meilleure approximation de la vraie matrice. Cela est dû au fait que nous n'avons pas pu trouver une expression explicite de la matrice p dans ce contexte spécifique.

Les résultats de notre étude numérique révèlent que la qualité de l'estimation des paramètres du modèle exposé-infecté dépend essentiellement de la quantité de données disponibles. Malheureusement, le jeu de données fourni par l'ARS de Mayotte s'est avéré insuffisant pour obtenir des estimations fiables et précises des paramètres clés du modèle, tels que $λ, μ, α$ et $ν$, et pour saisir pleinement leur influence sur la dynamique de transmission de la fièvre typhoïde à Mayotte.

## 4.5 Conclusion

Dans ce chapitre, nous avons proposé une modélisation de la transmission de la fièvre typhoïde basée sur un modèle à deux compartiments (exposé-infecté), permettant de prendre en compte le facteur de la durée d'incubation de la maladie. Nous avons commencé par décrire notre modèle puis étudié ses propriétés mathématiques telles que son comportement moyen et asymptotique, nous avons également illustré ses propriétés par des simulations numériques permettant de quantifier le comportement qualitatif du modèle. Ensuite, nous avons présenté notre méthode d'estimation qui consiste tout d'abord à exprimer les paramètres du modèle en termes de moments puis nous nous sommes placés dans le cadre de multi-chaîne de Markov cachée afin d'adapter l'algorithme de Baum-Welch à notre schéma d'observations et d'estimer nos paramètres du modèle exposé-infecté. Enfin, des simulations numériques ont été réalisées pour évaluer la performance de notre méthodologie d'estimation.

Cependant, l'application de notre approche au jeu de données de Mayotte fourni par l'ARS n'a pas été possible en raison il y a trop peu de données pour obtenir des résultats fiables. Ces résultats soulignent l'importance d'une surveillance épidémiologique continue et de la collecte de données pour mieux comprendre et surveiller la propagation de la maladie.
# Chapitre 5

## Effet de la pluviométrie sur l'évolution de la fièvre typhoïde à Mayotte

### Sommaire

<table><tr><td>5.1</td><td>Introduction</td><td>109</td></tr><tr><td>5.2</td><td>Description des données de précipitations</td><td>109</td></tr><tr><td>5.3</td><td>Détection de rupture</td><td>110</td></tr><tr><td>5.3.1</td><td>Méthode CUSUM</td><td>111</td></tr><tr><td>5.3.2</td><td>Application aux données de précipitations à Mayotte</td><td>111</td></tr><tr><td>5.4</td><td>Estimation des paramètres d'un modèle de Markov caché à deux régimes</td><td>114</td></tr><tr><td>5.5</td><td>Conclusion</td><td>118</td></tr></table>
**5.1 Introduction**

Dans ce chapitre, nous allons proposer une extension de notre modèle présenté dans le chapitre 3, appliqué au jeu de données de la fièvre typhoïde à Mayotte fourni par ARS. Cette fois, nous allons inclure l'effet de la pluviométrie sur l'évolution de la maladie. En effet, nous avons utilisé des données publiques sur les précipitations journalières à Mayotte [79] et nous allons étudier l'impact de ces conditions météorologiques sur notre modèle d'estimation. Nous avons constaté que la pluviométrie peut avoir un impact important sur la contamination des sources d'eau potable par les eaux de ruissellement, ce qui peut augmenter la transmission de la fièvre typhoïde.

La pluviométrie à Mayotte est relativement importante, avec une moyenne annuelle de précipitations d'environ 2000 mm. La saison des pluies s'étend généralement de novembre à avril, avec des précipitations abondantes et fréquentes. La saison sèche s'étend de mai à octobre, et pendant cette période, les précipitations sont rares. Ces conditions météorologiques peuvent avoir un impact sur la transmission des maladies infectieuses, notamment la fièvre typhoïde. En effet, les précipitations peuvent favoriser la contamination de l'environnement par la bactérie Salmonella Typhi, en particulier dans les zones où les conditions sanitaires sont précaires. C'est pourquoi il est important de considérer cette variable dans notre étude pour identifier son influence sur les paramètres de notre modèle, notamment le taux de contamination par l'environnement $v$.

Dans la section 5.2, nous commençons par décrire les données de précipitations utilisées dans cette étude. Ensuite, nous tentons de modéliser ces données de précipitations à Mayotte (Fig.5.1) en utilisant des modèles classiques de séries chronologiques telles que ARMA (en anglais : AutoRegressive Moving Average) et SARIMA (en anglais : Seasonnal AutoRegressive Integrated Moving Average) ou des techniques de lissage par des splines (voir Annexe E). Malheureusement, ces modèles n'ont pas réussi à capturer la variabilité des données. Nous avons donc décidé d'utiliser ces données uniquement pour détecter les deux régimes saisonniers de Mayotte, à savoir la saison sèche et la saison humide, à l'aide d'une méthode de détection de rupture classique présentée dans la section 5.3. Ensuite, dans la section 5.4, nous décrivons notre problème de Markov caché sous les deux régimes afin d'estimer les paramètres du modèle à partir du même jeu de données présenté dans 3.4.2.

**5.2 Description des données de précipitations**

Les données de précipitations utilisées dans cette étude proviennent de la source historique-meteo.net [79] et elles sont basées sur les enregistrements quotidiens de pluies régions de Mayotte, notamment Mamoudzou, Dapani, Chiconi, Mtsamboro et Bouéni. Ces données fournissent des informations sur les quantités de pluie enregistrées chaque jour dans ces différentes régions.

Afin d'obtenir une mesure agrégée représentative de l'ensemble de l'île, les données de précipitations de chaque région ont été combinées en calculant la moyenne des précipitations quotidiennes. Cette approche permet de déterminer une estimation globale des précipitations sur l'ensemble de Mayotte pour chaque jour.

Une première observation de la Fig.5.1 révèle une variabilité importante dans la quan-
tité de précipitations hebdomadaires à Mayotte, avec des fluctuations régulières présentant des pics et des creux. Cette observation suggère l'existence de deux saisons de précipitations distinctes. De plus, nous remarquons que les épidémies de fièvre typhoïde à Mayotte ont tendance à se produire préférentiellement pendant la saison humide de chaque année, indiquant que la pluviométrie peut jouer un rôle dans la propagation de la maladie.

<figure><img src="image_42.png" /><figcaption>FIGURE 5.1 - RR : La hauteur de précipitation en mm.</figcaption></figure>

- (A) : Moyenne des précipitations hebdomadaires à Mayotte entre 2018 et 2022.

- (B) : Nombre hebdomadaire des nouveaux cas déclarés de la fièvre typhoïde à Mayotte entre 2018 et 2022.

## 5.3 Détection de rupture

L'objectif de cette section est d'identifier les régimes de deux saisons climatiques à Mayotte (sèche et humide), plus précisément de détecter les instants où la série de précipitations présente des changements brusques sur une période donnée. Pour ce faire, nous avons utilisé la statistique CUSUM (en anglais : Cumulative Sum) qui représente une méthode non paramétrique (la distribution de la série est supposée inconnue et n'est pas nécessairement normale) de détection des ruptures.
### 5.3.1 Méthode CUSUM

La méthode CUSUM est une méthode statistique de détection de rupture utilisée pour détecter des changements dans une série chronologique. L'approche consiste à testier de manière séquentielle la présence de rupture en utilisant l'écriture récursive de la statistique de détection de Page [80]. Elle consiste à calculer la somme cumulative des différences entre les moyennes avant et après chaque instant dans la série chronologique. Plus précisément, soit $X_t$ la valeur de la série chronologique à l'instant $t$, et soit $\bar{X}_1$ et $\bar{X}_2$ les moyennes des valeurs de la série avant et après un instant de rupture supposé. La méthode CUSUM calcule la somme cumulative $S_t$ comme suit :

$$S_t = \max(0, S_{t-1} + (X_t - \bar{X}_1) - (X_t - \bar{X}_2)), \quad (5.1)$$

où $S_0 = 0$ est la condition initiale. Si la série ne présente pas de changement, la somme cumulative restera proche de zéro. En revanche, si un changement survient, la somme cumulative augmentera et dépassera un seuil $h$ prédéfini, ce qui indiquera la détection d'une rupture dans la série chronologique. La méthode conventionnelle pour choisir un seuil de détection constant $h$ est basée sur l'inégalité de Wald [81]. Effectivement, selon Egea-Roca et al [82], le seuil de détection $h$ doit être choisi de manière à respecter l'inégalité : $h \le -\log(\alpha)$, où $\alpha$ est le niveau de risque souhaité pour la détection de rupture. En fixant un $\alpha$ donné (par exemple 0.05), on peut donc calculer la valeur maximale possible de $h$. Si le seuil choisi est supérieur à cette valeur maximale, alors le risque de fausse détection augmente. En pratique, il est souvent recommandé de choisir un seuil de détection légèrement inférieur à cette valeur maximale pour réduire ce risque.

Il existe plusieurs variantes de la méthode CUSUM en fonction des statistiques utilisées pour détecter les ruptures. Les deux versions les plus couramment utilisées sont la méthode CUSUM par moyenne et la méthode CUSUM par variance.

La méthode CUSUM par moyenne est basée sur la différence entre la moyenne actuelle et la moyenne précédente de la série de données. Elle est donnée par la formule 5.1.

La méthode CUSUM par variance, quant à elle, est basée sur la différence entre la variance actuelle et la variance précédente de la série de données. Elle est donnée par la formule suivante :

$$S_t = \max(0, S_{t-1} + (X_t - \sigma_1^2) - (X_t - \sigma_2^2)),$$

où $\sigma_1^2$ et $\sigma_2^2$ sont les variances de la série avant et après l'instant $t$. De même, si la statistique CUSUM dépasse le seuil $h$ fixé, cela signifie qu'une rupture a été détectée.

Il existe également une méthode CUSUM combinée qui utilise à la fois la différence entre la moyenne et la différence entre la variance de la série de données pour détecter les ruptures. Cette méthode peut être plus robuste dans certaines situations où les ruptures peuvent affecter à la fois la moyenne et la variance de la série de données.

### 5.3.2 Application aux données de précipitations à Mayotte

Maintenant que nous avons présenté la méthode de détection CUSUM, nous pouvons l'appliquer pour sélectionner les instants des ruptures et ainsi identifier les régimes saisonniers de Mayotte.

En pratique, on utilise le package 'changepoint' du logiciel R [83], il offre principalement trois fonctions pour la détection de rupture, notamment :

* La fonction cpt.mean pour la détection de rupture dans des séries chronologiques variées basées sur le changement de la moyenne.
• La fonction **cpt.var** pour la détection de rupture dans des séries chronologiques univariées basées sur le changement de la variance.

• La fonction **cpt.meanvar** pour la détection de rupture dans des séries chronologiques univariées basées sur le changement de la moyenne et de la variance simultanément.

Pour chaque fonction, Les changements sont déterminés à l'aide de la méthode fournie, qui peut être un point de changement unique (AMOC) ou des points de changement multiples utilisant des méthodes exactes (PELT ou SegNeigh) ou approximatives (Bin-Seg) :

• La méthode **AMOC** (en anglais : Analysis of MOment Change) est une technique utilisée pour localiser un unique point de rupture dans une série d'observations. Elle est aussi connue sous le nom de CUSUM de base. Le point de rupture est alors déterminé comme le point où la somme cumulée atteint son maximum.

• **PELT** (en anglais : Pruned Exact Linear Time) : il s'agit d'une méthode exacte qui cherche à minimiser une fonction de coût basée sur la somme des carrés des écarts entre les données et les prévisions à chaque segment.

• **BinSeg** (en anglais : Binary Segmentation) : il s'agit d'une méthode qui divise la série en sous-segments de taille égale jusqu'à atteindre la partition optimale qui minimise une fonction de coût basée sur les carrés des écarts entre les données et les prévisions à chaque segment.

Le choix de la méthode dépend donc des caractéristiques de la série à analyser et de l'objectif de l'analyse. Il est souvent recommandé de tester plusieurs méthodes (voir Annexe F) pour comparer les résultats et choisir celle qui convient le mieux.

Dans notre étude sur la série de précipitations à Mayotte $R := (R_t)_{0 \le t \le T}$ **5.1**, la méthode la plus performante pour détecter les changements de régime est la fonction **cpt.meanvar** couplée avec la méthode de **PELT**. Les résultats obtenus sont présentés dans la Table **5.1**, qui donne les dates des ruptures $T_i$ (en semaines) sur $J = [0, T]$, ainsi que les moyennes de précipitations associées (**Fig.5.2**).
<figure><img src="image_43.png" /><figcaption>FIGURE 5.2 – Détection des ruptures de données pluviométriques à Mayotte entre 2018 et 2022.</figcaption></figure>

<table><caption>TABLE 5.1 – Les régimes saisonniers à Mayotte et leurs RR moyen associé.</caption><thead><tr><th>Période par semaines J<sub>i</sub> =]T<sub>i</sub>, T<sub>i+1</sub>]</th><th>RR moyen</th><th>Saison</th></tr></thead><tbody><tr><td>[0, 15]</td><td>12.85</td><td>Humide</td></tr><tr><td>[15, 44]</td><td>0.82</td><td>Sèche</td></tr><tr><td>[44, 70]</td><td>8.85</td><td>Humide</td></tr><tr><td>[70, 96]</td><td>0.99</td><td>Sèche</td></tr><tr><td>[96, 117]</td><td>9.80</td><td>Humide</td></tr><tr><td>[117, 147]</td><td>1.11</td><td>Sèche</td></tr><tr><td>[147, 173]</td><td>7.34</td><td>Humide</td></tr><tr><td>[173, 208]</td><td>1.41</td><td>Sèche</td></tr><tr><td>[208, 226]</td><td>15.32</td><td>Humide</td></tr><tr><td>[226, 252]</td><td>0.30</td><td>Sèche</td></tr><tr><td>[252, 260]</td><td>8.19</td><td>Humide</td></tr></tbody></table>

Pour évaluer la significativité des différences entre les moyennes des taux de RR obtenus, nous avons recouru à une analyse de variance (ANOVA) à un facteur. Les résultats de cette ANOVA ont révélé une très faible valeur de p-value ($p < 2.10^{-16}$), indiquant que les moyennes sont statistiquement distinctes. Par la suite, nous avons procédé à des tests post-hoc, afin de cerner précisément quelles paires de groupes présentent des différences significatives entre elles. Parmi les tests post-hoc couramment utilisés, nous avons
employé le test de Tukey. Les résultats de ce test confirment l'existence de différences significatives entre les moyennes des groupes "Humide" et "Sec", comme illustré dans la Figure 5.4.

<figure><img src="image_44.png" /><figcaption>FIGURE 5.3 – Résultats du test Tukey.</figcaption></figure>

## 5.4 Estimation des paramètres d'un modèle de Markov caché à deux régimes

Après avoir identifié les instants des changements dans les données de pluviométrie à Mayotte et défini les différentes périodes correspondant à chaque régime, nous procédons maintenant à l'extension de notre modèle de Markov caché, tel que présenté dans la section 3.3.2. Cette fois-ci, nous incorporons les données de pluviométrie représentant les deux régimes saisonniers, dans le but d'estimer les trois paramètres du modèle spécifiques à chaque régime (λ, μ et ν). Notre objectif est de quantifier l'impact de cette variable météorologique sur la propagation de la fièvre typhoïde à Mayotte.

Nous distinguons deux régimes saisonniers, le régime humide et le régime sec, que
nous notons respectivement $J^H$ et $J^S = J \setminus J^H$. Nous utilisons la fonction $m(t)$ pour
représenter le régime en vigueur à l'instant $t \in J$ :

$$
m(t) = \begin{cases} 1 & \text{si } t \in J^H, \\ 2 & \text{si } t \in J^S. \end{cases}
$$

Afin de réutiliser la même méthodologie d'estimation présentée en détail dans la section **3.3**, nous réécrivons les caractéristiques de notre modèle de Markov caché sous les deux régimes (**Fig.5.4**). Dans ce cas, on note par

$$
p_{i,j}^{m(t)} = \mathbb{P}(I_{t+1}^{m(t+1)} = j | I_t^{m(t)} = i), \quad i,j \in \mathbb{N}, t \in J,
$$

la probabilité de transition de la chaîne $I_t$. Plus précisément,

$$
p_{i,j}^{m(t)} = \begin{cases} p_{i,j}^{(1)}(\lambda_1, \mu_1, \nu_1) & \text{si } t \in J^H, \\ p_{i,j}^{(2)}(\lambda_2, \mu_2, \nu_2) & \text{si } t \in J^S, \end{cases}
$$

où $p_{i,j}^{(1)}(\lambda_1, \mu_1, \nu_1)$ et $p_{i,j}^{(2)}(\lambda_2, \mu_2, \nu_2)$ sont les probabilités de transition pour les régimes humide et sec, respectivement, et sont données par le théorème **3.1**.

<figure><img src="image_45.png" /><figcaption>FIGURE 5.4 – Modèle de Markov caché à deux régimes.</figcaption></figure>

Les paramètres $M = (Q^{(1)}, \psi^{(1)}, \rho^{(1)}, Q^{(2)}, \psi^{(2)}, \rho^{(2)})$ du HMM $(Z_t, O_t, R_t)$ sont données comme suit :

i) La matrice de transition de la chaîne $(Z_t^{m(t)})_{t \ge 1} = (I_{t-1}^{m(t-1)}, I_t^{m(t)})_{t \ge 1}$ :

$$
Q_{(i,j),(i',j')}^{m(t)} = p_{i',j'}^{m(t)} \delta_{i'=j},
$$

pour $(i, j)$, $(i', j') \in \mathbb{N}^2$.

ii) La probabilité d'émission :

$$
\psi_{(i,j)}^{m(t)}(o) = \mathbb{P}(O_t = o | Z_t^{m(t)} = (i,j)) = \frac{p_{i,(j,o)}^{m(t)}}{p_{i,j}^{m(t-1)}},
$$

pour $i, j, o \in \mathbb{N}$, où $p_{i,(j,o)} = p_{i,(j,o)}(\Delta t)$ est caractérisée par l'équation **3.2**.

iii) La distribution initiale :

$$
\rho_{i,j}^{m(1)} = p_{i,j}^{m(0)} \pi_i^{m(0)},
$$
pour $i, j \in \mathbb{N}$ où $\pi_i^{m(0)}$ est la loi initiale de la chaîne $I_t$ sous le régime $m(0)$.

En utilisant ces paramètres, nous pouvons exprimer la loi jointe de $(Z, O)$ par la formule suivante, pour tout $z_{1:T} \in (\mathbb{N}^2)^T$ et $O_{1:T} \in \mathbb{N}^T$ :

$$ 
\begin{aligned}
\mathbb{P}_M(Z_{1:T}^{m(1):m(T)} = z_{1:T}, O_{1:T} = o_{1:T}) &= \rho^{m(1)}(z_1) \prod_{t=1}^{T-1} Q^{m(t)}(z_t, z_{t+1}) \prod_{t=1}^{T} \psi_{z_t}^{m(t)}(o_t) \\
&= \rho^{m(1)}(z_1) \prod_{t \in J^H}^{T-1} Q^{(1)}(z_t, z_{t+1}) \prod_{t \in J^H}^{T} \psi_{z_t}^{(1)}(o_t) \\
& \quad \prod_{t \in J^S}^{T-1} Q^{(2)}(z_t, z_{t+1}) \prod_{t \in J^S}^{T} \psi_{z_t}^{(2)}(o_t).
\end{aligned}
$$

Ainsi, la fonction auxiliaire de l'algorithmé EM s'écrit comme suit :

$$ \Gamma_{O,R}(M, M') = \Gamma_{O,R}(M^{(1)}, M'^{(1)}) + \Gamma_{O,R}(M^{(2)}, M'^{(2)}), $$

où $M^{(i)}$ pour $i = 1, 2$ représente les paramètres du HMM sous le régime humide et sec, respectivement.

Par conséquent, en utilisant la même technique d'optimisation décrite dans la preuve du théorème 3.3, nous pouvons déduire que les ré-estimateurs de Baum-Welch adapté, étant donné le modèle $M = (M^{(1)}, M^{(2)})$, les observations $O_{1:T}$ et les régimes $J^H$ et $J^S$, sont données par :

$$ Q_{(i,j),(i',j')}^{(1)} = \frac{\sum_{t \in J^H} \xi_{(i,j),(i',j')}^{m(t)}(t)}{\sum_{t \in J^H} \gamma_{(i,j)}^{m(t)}(t)} \delta_{i'=j}, $$

$$ Q_{(i,j),(i',j')}^{(2)} = \frac{\sum_{t \in J^S} \xi_{(i,j),(i',j')}^{m(t)}(t)}{\sum_{t \in J^S} \gamma_{(i,j)}^{m(t)}(t)} \delta_{i'=j}, $$

$$ \psi_{(i,j)}^{(1)}(y) = \frac{\sum_{t \in J^H} \mathbb{1}_{o_t = o} \gamma_{(i,j)}^{m(t)}(t)}{\sum_{t \in J^H} \gamma_{(i,j)}^{m(t)}(t)}, $$

$$ \psi_{(i,j)}^{(2)}(y) = \frac{\sum_{t \in J^S} \mathbb{1}_{o_t = o} \gamma_{(i,j)}^{m(t)}(t)}{\sum_{t \in J^S} \gamma_{(i,j)}^{m(t)}(t)}, $$

$$ \rho_{i,j}^{m(1)} = \gamma_{(i,j)}^{m(1)}(1), $$

$$ p_{i,j}^{(1)} = \frac{\sum_{i' \in \mathbb{N}} Q_{(i',i),(i,j)}^{(1)}}{\sum_{i' \in \mathbb{N}} \mathbb{1}_{Q_{(i',i),(i,j)}^{(1)} \neq 0}}, $$

$$ p_{i,j}^{(2)} = \frac{\sum_{i' \in \mathbb{N}} Q_{(i',i),(i,j)}^{(2)}}{\sum_{i' \in \mathbb{N}} \mathbb{1}_{Q_{(i',i),(i,j)}^{(2)} \neq 0}}, $$

où

$$ \gamma_{(i,j)}^{m(t)} = \frac{\alpha_{(i,j)}^{m(t)}(t) \beta_{(i,j)}^{m(t)}(t)}{\sum_{i,j \in \mathbb{N}} \alpha_{(i,j)}^{m(t)}(t) \beta_{(i,j)}^{m(t)}(t)}, $$

$$ \xi_{(i,j),(i',j')}^{m(t)}(t) = \frac{\alpha_{(i,j)}^{m(t)}(t) p_{(i',j')}^{m(t)}(t) \psi_{(i',j')}^{m(t+1)}(o_{t+1}) \beta_{(i',j')}^{m(t+1)}(t+1)}{\sum_{i,j \in \mathbb{N}} \alpha_{(i,j)}^{m(t)}(t) \beta_{(i,j)}^{m(t)}(t)} \delta_{i'=j}, $$
et les probabilités forward et backward $\alpha$ et $\beta$ sont définies récursivement pour $i, j$ dans $\mathbb{N}$ et $2 \le t \le T$ par

$$
\left\{
\begin{array}{l}
\alpha_{(i,j)}^{m(1)}(1) = p_{i,j}^{m(0)} \pi_i^{m(0)} \psi_{(i,j)}^{m(1)}(o_1), \\
\alpha_{(i,j)}^{m(t)}(t) = \psi_{(i,j)}^{m(t)}(o_t) p_{i,j}^{m(t)} \sum_{i' \in \mathbb{N}} \alpha_{(i',i)}^{m(t-1)}(t-1), \\
\beta_{(i,j)}^{m(T)}(T) = 1, \\
\beta_{(i,j)}^{m(t-1)}(t-1) = \sum_{j' \in \mathbb{N}} p_{j,j'}^{m(t)} \psi_{(j,j')}(o_t) \beta_{(j,j')}^{m(t)}(t).
\end{array}
\right.
$$

Enfin, les estimations des paramètres dans chaque régime sont obtenus à partir du théorème **3.2** : dans le régime humide $J^H$ :

$$ (\hat{\lambda}_1, \hat{\mu}_1, \hat{\nu}_1) = g(p_{0,0}^{(1)}, p_{0,1}^{(1)}, p_{1,0}^{(1)}), $$

et dans le régime sec $J^S$ :

$$ (\hat{\lambda}_2, \hat{\mu}_2, \hat{\nu}_2) = g(p_{0,0}^{(2)}, p_{0,1}^{(2)}, p_{1,0}^{(2)}). $$

Pour la partie numérique, nous avons choisi les mêmes intervalles initiaux pour $\lambda_1^{(0)}$ et $\lambda_2^{(0)}$ dans l'intervalle $[0.05, 0.08]$, et $\mu_1^{(0)}$ et $\mu_2^{(0)}$ dans l'intervalle $[0.11, 0.25]$. Cependant, nous avons spécifié des intervalles différents pour $\nu_1^{(0)}$ dans $[0.015, 0.03]$ et $\nu_2^{(0)}$ dans $[0.0015, 0.003]$, car le paramètre $\nu$ (le taux de contamination par l'eau) est directement lié à la pluviométrie et est considéré comme étant le plus représentatif pour quantifier l'impact de cette variable météorologique sur la transmission de la fièvre typhoïde à Mayotte. Les autres paramètres de l'algorithme sont détaillés dans la section **3.4.2**.

Enfin, les estimations des paramètres qui maximisent la vraisemblance des observations pour ce modèle HMM à deux régimes sont les suivantes :

Dans le régime humide :

$$ (\hat{\lambda}_1, \hat{\mu}_1, \hat{\nu}_1) = (0.066, 0.21, 0.05), $$

Dans le régime sec :

$$ (\hat{\lambda}_2, \hat{\mu}_2, \hat{\nu}_2) = (0.061, 0.1, 0.02). $$

Les résultats numériques montrent que dans le régime humide, le taux de contamination de personne à personne est plus élevé ($\hat{\lambda}_1 = 0.066$) que dans le régime sec ($\hat{\lambda}_2 = 0.061$), ce qui peut être expliqué par le fait que l'humidité favorise la survie des bactéries et donc augmente le risque de contamination. De plus, le taux d'isolement ou de déclaration des cas de typhoïde est plus élevé dans le régime humide ($\hat{\mu}_1 = 0.21$) que dans le régime sec ($\hat{\mu}_2 = 0.1$), ce qui peut être dû à une sensibilisation accrue de la population et/ou une surveillance renforcée pendant la saison des pluies. Enfin, le taux de contamination par l'eau est également plus élevé dans le régime humide ($\hat{\nu}_1 = 0.05$) que dans le régime sec ($\hat{\nu}_2 = 0.02$), ce qui suggère que la pluviométrie a un impact significatif sur la transmission de la typhoïde à Mayotte. Cela est cohérent avec la littérature scientifique [84, 85, 86] et les remarques des médecins de la fièvre typhoïde à Mayotte. En effet, lorsque qu'il pleut, l'eau de pluie peut s'écouler dans les rivières, les lacs, les puits ou les réservoirs d'eau. Si ces sources d'eau ne sont pas correctement traitées ou protégées, elles peuvent être contaminées par des bactéries et autres agents pathogènes provenant des eaux usées, des animaux ou des sols.

De plus, la pluie peut également entraîner des inondations, ce qui peut contaminer les
sources d'eau potable par les eaux usées ou les matières fécales des animaux. Dans les zones où l'assainissement et le traitement de l'eau sont insuffisants, les épisodes de pluie peuvent donc augmenter la contamination de l'eau par les bactéries responsables de la typhoïde, ce qui peut expliquer en partie l'impact de la pluie sur la transmission de cette maladie.

Cependant, il est important de noter que ces résultats ne sont basés que sur un échantillon limité de données et peuvent être améliorés avec des données plus nombreuses.

En conclusion, ces résultats montrent l'importance de prendre en compte les variables météorologiques dans l'étude de la transmission des maladies infectieuses hydriques et soulignent la nécessité de renforcer la surveillance et la prévention de la typhoïde pendant la saison des pluies.

## 5.5 Conclusion

Au cours de ce chapitre, nous avons étudié l'impact de la pluviométrie et de la saisonnalité sur la transmission de la fièvre typhoïde à Mayotte en se basant sur une méthode de détection de rupture, à savoir CUSUM et une extension de notre HMM 3.3.2.

Nous avons commencé par une analyse exploratoire des données, en traçant les graphiques de la pluviométrie et des cas de typhoïde au fil du temps. Ensuite, nous avons utilisé la méthode CUSUM pour détecter les changements ou bien les ruptures dans les données de précipitations à Mayotte et identifier les différentes périodes de deux régimes saisonnières, humide et sec.

Nous avons ensuite décrit notre modèle HMM à deux régimes afin d'estimer les trois paramètres $\lambda, \mu$ et $\nu$ dans chaque régime et caractériser et comparer la transmission de la maladie par rapport au changement de saison. Les estimations des paramètres pour ces deux régimes ont montré que le taux de contamination par l'eau était plus élevé dans le régime humide que dans le régime sec, suggérant ainsi que la pluviométrie peut augmenter la contamination de l'eau par les bactéries. Cette observation explique également les pics épidémiques enregistrés pendant la période pluvieuse chaque année.

Enfin, nous avons discuté des résultats obtenus et de leur pertinence dans le contexte de la lutte contre la typhoïde. Nous avons souligné l'importance de prendre en compte la pluviométrie et la saisonnalité dans la conception des stratégies de lutte contre la typhoïde, et avons conclu que les deux méthodes de détection de rupture étudiées peuvent être des outils utiles dans ce domaine.
Chapitre 6

Perspectives

Dans nos perspectives de recherche, nous souhaitons approfondir notre compréhension de l'effet de la variation des niveaux de précipitations sur l'évolution de la fièvre typhoïde. Nous nous pencherons sur l'identification d'un lien potentiel entre les précipitations et l'apparition des pics épidémiques de la maladie. Nos observations (voir Fig.6.1) ont mis en évidence que le niveau moyen des précipitations à Mayotte est variable, avec des années présentant des pics épidémiques plus importants, comme observé en 2022 avec des niveaux de précipitations élevés par rapport aux années précédentes.

Nous avons également remarqué qu'après chaque épisode de fortes pluies, il existe un décalage temporel suivi d'une augmentation des cas de fièvre typhoïde. Ce décalage pourrait s'expliquer par plusieurs facteurs potentiels, tels que la période d'incubation et de contagiOSITÉ de la maladie, ainsi que le processus de développement et de multiplication des bactéries nécessitant un certain laps de temps. Il est possible que les conditions environnementales créées lors des périodes de fortes pluies favorisent la survie et la propagation de la bactérie Salmonella typhi. Ainsi, le délai nécessaire pour que ces conditions environnementales favorables conduisent à une augmentation significative des cas pourrait expliquer le décalage observé. Ces constatations suggèrent donc que d'autres facteurs doivent être pris en compte dans notre analyse statistique.

<figure><img src="image_46.png" /><figcaption>FIGURE 6.1 – Présentation de la moyenne des précipitations et des cumuls de nouveaux cas déclarés à Mayotte de 2018 à 2022.</figcaption></figure>
En outre, il serait intéressant, en fonction des données disponibles dans le futur, d'ex-
plorer l'influence d'autres variables liées à l'accès à l'eau sur la concentration des bactéries
responsables de la fièvre typhoïde, ainsi que sur le processus de transmission par l'eau
contaminée. Parmi ces variables, nous pourrions nous concentrer sur la température, le
pH et la conductivité de l'eau, qui sont des indicateurs importants de sa qualité (Fig.6.2).
La température et le pH de l'eau peuvent jouer un rôle significatif dans la survie et la mul-
tiplication des bactéries. Des études antérieures [87, 88] ont montré que des températures
plus élevées peuvent favoriser la croissance de Salmonella typhi. D'autres études [89, 90]
suggèrent que des valeurs de pH spécifiques peuvent favoriser la survie de Salmonella
typhi.

<figure><img src="image_47.png" /><figcaption>FIGURE 6.2 – Différents facteurs associés à la qualité de l'eau.</figcaption></figure>

Une autre perspective intéressante consisterait à examiner les facteurs socio-économiques et comportementaux qui pourraient influencer la transmission de la fièvre typhoïde. Par exemple, nous pourrions étudier l'impact de la densité de population, des conditions d'assainissement, de l'accès à l'eau potable et des pratiques d'hygiène sur la prévalence de la maladie. Comprendre comment ces facteurs interagissent avec les variables météorologiques pourrait permettre de développer des stratégies de prévention et de contrôle plus efficaces.

Une autre piste de recherche prometteuse consisterait à réaliser une étude théorique
sur la stabilité du modèle PDMP présenté dans la section 2.3. Cette étude permettrait
de déterminer le comportement stationnaire du modèle, d'identifier les différents régimes
possibles et d'étudier les caractéristiques de la population à l'équilibre. Une analyse appro-
fondie de la stabilité nous permettrait de mieux comprendre les dynamiques de la fièvre
typhoïde et d'apporter des informations essentielles pour la prévention et le contrôle de
la maladie.

En outre, nous pourrions envisager la mise en œuvre de procédures de contrôle sto-
chastique optimal dans différents scénarios épidémiques. Ces approches de contrôle per-
mettraient d'identifier les stratégies les plus efficaces pour réduire la propagation de la
maladie. En prenant en compte les paramètres spécifiques du modèle PDMP et en utilisant des techniques d'optimisation adaptées, nous pourrions concevoir des interventions ciblées visant à minimiser l'impact de la fièvre typhoïde dans la population étudiée.

Enfin, il serait intéressant d'étendre notre étude aux autres maladies infectieuses présentes à Mayotte, en incluant d'autres pathologies présentes sur le territoire étudié, telles que les maladies de la peau, afin d'avoir une vision plus complète de la situation sanitaire. L'intégration de paramètres liés à l'exploitation du Système national des données de santé (SNDS) permettrait une généralisation du modèle à d'autres territoires que Mayotte. Cela permettrait d'élargir l'applicabilité de notre recherche et d'adapter les résultats à des contextes différents, contribuant ainsi à une meilleure compréhension et gestion globale des maladies infectieuses.
## B Matrice Jacobienne $Dg$

Nous calculons la matrice jacobienne de $Dg$. Tout d'abord, les dérivées partielles des fonctions $q$, $r$ et $u$ présentées dans la preuve du théorème 3.2 sont données par :

$$ \frac{\partial q}{\partial p_{0,0}} = - \frac{p_{0,1} (\ln(p_{0,0}) + 1)}{(p_{0,0} \ln(p_{0,0}))^2} W \left( \frac{\binom{p_{0,0}}{p_{0,1}} \ln(p_{0,0})}{p_{0,1}} \right) \\ + \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})} (\ln(p_{0,0}) + 1) \left( \frac{1}{p_{0,1}} + \frac{1}{p_{0,0} \ln(p_{0,0})} \right) \frac{W \left( \frac{\binom{p_{0,0}}{p_{0,1}}}{p_{0,1}} \ln(p_{0,0}) \right)}{W \left( \frac{\binom{p_{0,0}}{p_{0,1}}}{p_{0,1}} \ln(p_{0,0}) \right) + 1} $$

Puisque

$$ W \left( \frac{\binom{p_{0,0}}{p_{0,1}} \ln(p_{0,0})}{p_{0,1}} \right) = q \frac{p_{0,0} \ln(p_{0,0})}{p_{0,1}}, \qquad (1) $$

on obtient

$$ \begin{aligned} \frac{\partial q}{\partial p_{0,0}} &= \frac{q (\ln(p_{0,0}) + 1)}{qp_{0,0} \ln(p_{0,0}) + p_{0,1}} \left( 1 + \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})} - \frac{qp_{0,0} \ln(p_{0,0}) + p_{0,1}}{p_{0,0} \ln(p_{0,0})} \right) \\ &= \frac{q (\ln(p_{0,0}) + 1)}{qp_{0,0} \ln(p_{0,0}) + p_{0,1}} (1 - q). \end{aligned} $$

Ensuite, on a

$$ \frac{\partial q}{\partial p_{0,1}} = \frac{1}{p_{0,0} \ln(p_{0,0})} W \left( \frac{\binom{p_{0,0}}{p_{0,1}} \ln(p_{0,0})}{p_{0,1}} \right) \\ - \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})} \left( \frac{p_{0,0} \ln(p_{0,0})}{(p_{0,1})^2} + \frac{1}{p_{0,1}} \right) \frac{W \left( \frac{\binom{p_{0,0}}{p_{0,1}}}{p_{0,1}} \ln(p_{0,0}) \right)}{W \left( \frac{\binom{p_{0,0}}{p_{0,1}}}{p_{0,1}} \ln(p_{0,0}) \right) + 1}, $$

et en utilisant l'équation (1) nous obtenons

$$ \begin{aligned} \frac{\partial q}{\partial p_{0,1}} &= \frac{q}{p_{0,1}} - \frac{q}{qp_{0,0} \ln(p_{0,0}) + p_{0,1}} \left( \frac{p_{0,0} \ln(p_{0,0})}{p_{0,1}} + 1 \right) \\ &= \frac{q}{p_{0,1} \left( q + \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})} \right)} (q - 1), \end{aligned} $$
et enfin on a

$$ \frac{\partial q}{\partial p_{1,0}} = 0. $$

Les expressions des dérivées partielles de $r$ sont données comme suit :

$$ \begin{aligned} \frac{\partial r}{\partial p_{0,0}} &= \frac{1}{p_{0,0} \ln(q)} - \frac{\ln(p_{0,0})}{(\ln(q))^2 q} \frac{\partial q}{\partial p_{0,0}} \\ &= \frac{r}{p_{0,0} \ln(p_{0,0})} - r \frac{(\ln(p_{0,0}) + 1)(1-q)}{\ln(q) (qp_{0,0} \ln(p_{0,0}) + p_{0,1})}, \\ \frac{\partial r}{\partial p_{0,1}} &= \frac{\ln(p_{0,0})}{(\ln(q))^2 q} \frac{\partial q}{\partial p_{0,1}} \\ &= r \frac{(q-1)}{p_{0,1} \ln(q) \left(q + \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})}\right)}, \\ \frac{\partial r}{\partial p_{1,0}} &= 0. \end{aligned} $$

Les dérivées partielles de $u$ s'expriment respectivement comme suit

$$ \begin{aligned} \frac{\partial u}{\partial p_{0,0}} &= \frac{p_{1,0}}{p_{0,0}} = \frac{1-u}{p_{0,0}}, \\ \frac{\partial u}{\partial p_{0,1}} &= 0, \\ \frac{\partial u}{\partial p_{1,0}} &= -\frac{1}{p_{0,0}}. \end{aligned} $$

En définissant les fonctions suivantes

$$ \begin{aligned} \alpha &= \frac{q (\ln(p_{0,0}) + 1) (q - 1)}{\Delta t (q p_{0,0} \ln(p_{0,0}) + p_{0,1}) (q - u)}, \\ \beta &= \frac{1 - q}{\Delta t p_{0,1} \left(q + \frac{p_{0,1}}{p_{0,0} \ln(p_{0,0})}\right) (q - u)}, \\ \eta &= \frac{1 - u}{\Delta t p_{0,0} u (q - u)}, \end{aligned} $$
on peut réécrire de façon concise les dérivées partielles de $\lambda$, $\mu$ et $\nu$ (de la fonction $g$) comme suit :

$$ 
\begin{aligned}
\frac{\partial \lambda}{\partial p_{0,0}} &= \eta(q-1) - \frac{\alpha}{q}(1-q) + \lambda \frac{\alpha \Delta t (q-u)}{1-q} + \lambda \alpha \Delta t + \lambda \frac{1-u}{p_{0,0}(q-u)} \\
&= \lambda(1-u) \left( \frac{\alpha \Delta t}{1-q} + \frac{1}{p_{0,0}(q-u)} \right) - (1-q) \left( \frac{\alpha}{q} + \eta \right),
\end{aligned}
$$

$$ \frac{\partial \lambda}{\partial p_{0,1}} = \lambda \frac{\beta \Delta t q (1-u)}{1-q} - \beta (1-q), $$

$$ \frac{\partial \lambda}{\partial p_{1,0}} = - \frac{\lambda}{p_{0,0}(q-u)} + \eta \frac{1-q}{1-u}, $$

$$ 
\begin{aligned}
\frac{\partial \mu}{\partial p_{0,0}} &= -\eta(1-u) - \frac{\alpha}{q}(1-u) - \frac{\mu}{p_{0,0}} + \mu\alpha\Delta t + \frac{\mu(1-u)}{p_{0,0}(q-u)} \\
&= \mu \left( \alpha\Delta t + \frac{1-q}{p_{0,0}(q-u)} \right) - (1-u) \left( \frac{\alpha}{q} + \eta \right),
\end{aligned}
$$

$$ \frac{\partial \mu}{\partial p_{1,0}} = \mu \beta \Delta t q - \beta (1-u), $$

$$ \frac{\partial \mu}{\partial p_{0,1}} = \mu \frac{1-q}{p_{0,0}(u-1)(q-u)} + \eta, $$

$$ 
\begin{aligned}
\frac{\partial \nu}{\partial p_{0,0}} &= \frac{\nu}{p_{0,0} \ln(p_{0,0})} + \nu \frac{\alpha \Delta t (q-u)}{q \ln(q)} + \nu (1-u) \left( \frac{\alpha \Delta t}{1-q} + \frac{1}{p_{0,0}(q-u)} \right) \\
&\quad - r (1-q) \left( \frac{\alpha}{q} + \eta \right) \\
&= \nu \left( \frac{1}{p_{0,0} \ln(p_{0,0})} + \frac{\alpha \Delta t (q-u)}{q \ln(q)} + (1-u) \left( \frac{\alpha \Delta t}{1-q} + \frac{1}{p_{0,0}(q-u)} \right) \right) \\
&\quad - r (1-q) \left( \frac{\alpha}{q} + \eta \right),
\end{aligned}
$$

$$ \frac{\partial \nu}{\partial p_{0,1}} = \nu \beta \Delta t \left( \frac{q(1-u)}{1-q} - \frac{q-u}{\ln(q)} \right) - r \beta (1-q), $$

$$ \frac{\partial \nu}{\partial p_{1,0}} = - \frac{\nu}{p_{0,0}(q-u)} + r \eta \frac{1-q}{1-u}. $$
C Log-vraisemblance $l(O|M^n)$

<table><caption>TABLE 1 – La log-vraisemblance $l(O|M^n)$ pour différentes valeurs initiales des paramètres.</caption><thead><tr><th>(λ<sup>(0)</sup>, μ<sup>(0)</sup>, ν<sup>(0)</sup>)</th><th>l(Y|M<sup>n</sup>)</th></tr></thead><tbody><tr><td>(0.0500, 0.1100, 0.0150)</td><td>-394.1107</td></tr><tr><td>(0.0516, 0.1174, 0.0157)</td><td>-394.4322</td></tr><tr><td>(0.0532, 0.1247, 0.0166)</td><td>-394.0771</td></tr><tr><td>(0.0546, 0.1321, 0.0174)</td><td>-393.7732</td></tr><tr><td>(0.0563, 0.1395, 0.0182)</td><td>-394.1762</td></tr><tr><td>(0.0579, 0.1468, 0.0189)</td><td>-394.5896</td></tr><tr><td>(0.0595, 0.1542, 0.0197)</td><td>-394.9808</td></tr><tr><td>(0.0611, 0.1616, 0.0205)</td><td>-395.3653</td></tr><tr><td>(0.0626, 0.1689, 0.0213)</td><td>-395.7492</td></tr><tr><td>(0.0642, 0.1763, 0.0221)</td><td>-396.1328</td></tr><tr><td>(0.0658, 0.1837, 0.0229)</td><td>-396.5154</td></tr><tr><td>(0.0674, 0.1911, 0.0237)</td><td>-397.0610</td></tr><tr><td>(0.0689, 0.1984, 0.0245)</td><td>-397.4057</td></tr><tr><td>(0.0705, 0.2058, 0.0253)</td><td>-396.1538</td></tr><tr><td>(0.0721, 0.2132, 0.0261)</td><td>-395.9935</td></tr><tr><td>(0.0737, 0.2205, 0.0268)</td><td>-396.1987</td></tr><tr><td>(0.0753, 0.2279, 0.0276)</td><td>-396.4099</td></tr><tr><td>(0.0768, 0.2353, 0.0284)</td><td>-396.6318</td></tr><tr><td>(0.0784, 0.2426, 0.0292)</td><td>-396.8629</td></tr><tr><td>(0.0800, 0.2500, 0.0300)</td><td>-397.1020</td></tr></tbody></table>
# D Fonction génératrice des moments du processus $(E_t, I_t)_{t \ge 0}$

La fonction génératrice des moments du processus $(E_t, I_t)_{t \ge 0}$ est donnée par : pour tout $(r, s, t) \in \mathbb{R}^2 \times [0, +\infty[$

$$G_{(e,i)}(r, s, t) = \sum_{k,l=0}^{+\infty} r^k s^l p_{(e,i),(k,l)}(t).$$

En utilisant l'équation de Kolmogorov progressive (4.2), on trouve :

$$
\begin{align*}
\frac{\partial G_{(e,i)}(r, s, t)}{\partial t} &= \sum_{k,l=0}^{\infty} r^k s^l \sum_{(e',i')} p_{(e,i),(e',i')} Q_{(e',i'),(k,l)} \\
&= \sum_{k,l=0} r^{k+1} s^l (\lambda l + \nu) p_{(e,i),(k,l)} + \sum_{k=1,l=0} r^{k-1} s^{l+1} \alpha k p_{(e,i),(k,l)} \\
&\quad + \sum_{k=0,l=1} r^k s^{l-1} \mu l p_{(e,i),(k,l)} - \sum_{k,l=0} r^k s^l ((\lambda + \mu)l + \alpha k + \nu) p_{(e,i),(k,l)} \\
&= \lambda s r \sum_{k=0,l=1} r^k l s^{l-1} p_{(e,i),(k,l)} + \nu r \sum_{k,l=0} r^k s^l p_{(e,i),(k,l)} + \alpha s \sum_{k=1,l=0} k r^{k-1} s^l p_{(e,i),(k,l)} \\
&\quad + \mu \sum_{k=0,l=1} r^k l s^{l-1} p_{(e,i),(k,l)} - (\lambda + \mu)s \sum_{k=0,l=1} r^k l s^{l-1} p_{(e,i),(k,l)} \\
&\quad - \alpha r \sum_{k=1,l=0} k r^{k-1} s^l p_{(e,i),(k,l)} - \nu \sum_{k,l=0} r^k s^l p_{(e,i),(k,l)} \\
&= \alpha (s-r) \sum_{k=1,l=0} k r^{k-1} s^l p_{(e,i),(k,l)} + (\lambda s(r-1) - \mu(s-1)) \sum_{k=0,l=1} r^k l s^{l-1} p_{(e,i),(k,l)} \\
&\quad + \nu(r-1) \sum_{k,l=0} r^k s^l p_{(e,i),(k,l)}.
\end{align*}
$$

Ainsi, la fonction génératrice des moments satisfait l'équation différentielle partielle suivante :

$$\frac{\partial G_{(e,i)}(r, s, t)}{\partial t} = \alpha (s - r) \frac{\partial G}{\partial r} + (\lambda s(r - 1) - \mu(s - 1)) \frac{\partial G}{\partial s} \quad (2) \\ + \nu(r - 1) G_{(e,i)}(r, s, t),$$

avec la condition initiale $G_{(e,i)}(r, s, 0) = r^e s^i$.

Nous utilisons ensuite la méthode des caractéristiques pour résoudre cette équation.

Étape 1 : Les caractéristiques de l'équation des dérivées partielles 2 sont données par :

$$\begin{align*}
\frac{dr}{dt} &= -(\lambda s(r-1) - \mu(s-1)), \\
\frac{ds}{dt} &= -\alpha (s-r).
\end{align*}$$

En intégrant, on obtient,

$$\begin{align*}
k &= (s(\lambda + \mu - \lambda r) - \mu) e^{-(\lambda + \mu - \lambda r)t}, \\
k' &= \alpha (r - s) e^{-\alpha t}.
\end{align*}$$

Étape 2 : Tout au long des caractéristiques, $G_{(e,i)}$ satisfait

$$\frac{dG_{(e,i)}(r, s, t)}{dt} = \nu(1 - r)G_{(e,i)}.$$
L'intégration par rapport à s donne

$$
\begin{align}
G_{(e,i)}(r, s, t) &= C(k, k')e^{(1-r)\nu t} \nonumber \\
&= C((s(\lambda + \mu - \lambda r) - \mu)e^{-(\lambda+\mu-\lambda r)t}, \alpha(r-s)e^{-\alpha t})e^{(1-r)\nu t}. \tag{3}
\end{align}
$$

Nous devons maintenant utiliser la condition initiale $G_{(e,i)}(r, s, 0) = r^e s^i$ pour déterminer
la fonction C.

$$
r^e s^i = C((s(\lambda + \mu - \lambda r) - \mu), \alpha(r - s)).
$$

Posons

$$
\begin{cases}
w = s(\lambda + \mu - \lambda r) - \mu, \\
v = \alpha(r - s).
\end{cases}
$$

Ce système est identifiable si et seulement si

$$
r + s \neq \frac{\lambda + \mu}{\lambda} \quad \text{et} \quad s < \frac{\mu^2}{4\lambda(\lambda + \mu)}.
$$

Ainsi, pour $b \in \mathbb{R}$, tel que $0 \le r \le 1$, $0 \le s \le b < \frac{\mu^2}{4\lambda(\lambda + \mu)}$ et $r + s \ne \frac{\lambda + \mu}{\lambda}$, on a

$$
\lambda s^2 - s(\lambda + \mu - \frac{\lambda v}{\alpha}) + (\mu + w) = 0,
$$

ce qui implique,

$$
\begin{align*}
s &= \frac{\left(\lambda + \mu - \frac{\lambda v}{\alpha}\right) - \sqrt{\Delta}}{2\lambda}, \\
r &= \frac{\left(\lambda + \mu + \frac{\lambda v}{\alpha}\right) - \sqrt{\Delta}}{2\lambda},
\end{align*}
$$

avec

$$
\Delta = \left(\lambda + \mu - \frac{\lambda v}{\alpha}\right)^2 - 4\lambda(\mu + w).
$$

Par conséquence,

$$
C(w, v) = \left( \frac{(\lambda + \mu + \frac{\lambda v}{\alpha}) - \sqrt{\Delta}}{2\lambda} \right)^e \left( \frac{(\lambda + \mu - \frac{\lambda v}{\alpha}) - \sqrt{\Delta}}{2\lambda} \right)^i . \quad (4)
$$

En combinant les deux équations 3 et 4, on trouve l'expression finale de $G_{(e,i)}(r, s, t)$,

$$
G_{(e,i)}(r, s, t) = \left( \frac{(\lambda + \mu + \lambda(r-s)e^{-\alpha t}) - \sqrt{\Delta(r,s,t)}}{2\lambda} \right)^e \left( \frac{(\lambda + \mu - \lambda(r-s)e^{-\alpha t}) - \sqrt{\Delta(r,s,t)}}{2\lambda} \right)^i \times e^{(1-r)\nu t},
$$

où

$$
\Delta(r, s, t) = (\lambda + \mu - \lambda(r - s)e^{-\alpha t})^2 - 4\lambda(\mu + (s(\lambda + \mu - \lambda r) - \mu)e^{-(\lambda+\mu-\lambda r)t}).
$$
# E Séries chronologiques

Nous présentons ici différents modèles classiques de séries chronologiques utilisés pour
modéliser de manière paramétrique les précipitations à Mayotte.

## • Modèle SARIMA :

Le modèle SARIMA (*Seasonal Autoregressive Integrated Moving Average*) est une méthode statistique utilisée pour analyser et prévoir les séries chronologiques présentant des tendances saisonnières. Il combine les concepts de l'autorégression (AR), de la moyenne mobile (MA) et de l'intégration (I) pour capturer à la fois les variations saisonnières et les relations entre les valeurs observées successives. Le modèle SARIMA prend en compte les effets saisonniers, les tendances à court terme et les variations aléatoires pour fournir des prévisions précises sur les données temporelles.

Pour l'identification du modèle, nous utilisons l'approche de Box-Jenkins en suivant les étapes suivantes :

1. *Identification de la saisonnalité* : L'étape initiale consiste à examiner les données et à identifier la période de saisonnalité. Cela peut être effectué en analysant les graphiques de la série temporelle et en recherchant des motifs et des régularités saisonnières.

2. *Stationnarité* : Les modèles SARIMA nécessitent que la série temporelle soit stationnaire, c'est-à-dire que ses propriétés statistiques ne varient pas avec le temps. Si la série n'est pas stationnaire, des transformations telles que la différenciation peuvent être appliquées pour obtenir une série stationnaire.

3. *Identification des ordres AR, MA et saisonnier* : À cette étape, les ordres autoregressifs (AR) et de moyenne mobile (MA) sont identifiés en utilisant l'autocorrélation (ACF) et l'autocorrélation partielle (PACF). Ces graphiques aident à déterminer les valeurs des retards significatifs qui contribuent aux modèles AR et MA.

4. *Estimation du modèle* : Une fois que les ordres AR, MA et saisonnier ont été identifiés, le modèle SARIMA est estimé à l'aide de méthodes telles que la méthode des moindres carrés ou la méthode de maximisation de la vraisemblance. L'estimation fournit les coefficients du modèle qui captureront les caractéristiques de la série temporelle.

5. *Diagnostic du modèle* : Après l'estimation, il est important de diagnostiquer le modèle pour vérifier s'il satisfait aux hypothèses de la stationnarité résiduelle et de l'absence d'autocorrélation. Les graphiques des résidus, l'autocorrélation des résidus et les tests statistiques aident à évaluer la qualité du modèle.

6. *Réévaluation et itérations* : Si le modèle ne satisfait pas aux critères de diagnostic, il est nécessaire de réévaluer les ordres et d'ajuster le modèle. Ce processus d'itération est effectué jusqu'à ce qu'un modèle satisfaisant soit obtenu.

En suivant ces étapes, l'approche de Box-Jenkins permet d'identifier et d'estimer un modèle SARIMA approprié pour modéliser et prévoir les séries chronologiques avec des composantes saisonnières.

En pratique, nous utilisons le package fpp3, développé en tandem avec le manuel de Hyndman et Athanasopoulos [91]. Ce package permet l'installation et le chargement automatique d'autres packages essentiels pour la visualisation et l'analyse des séries temporelles.
Pour visualiser simultanément l'autocorrélation (ACF), l'autocorrélation partielle (PACF) et la série chronologique originale, nous utilisons la fonction "tsdisplay".

<figure><img src="image_48.png" /><figcaption>FIGURE 3 – Présentation de la série hebdomadaire des précipitations et de ses graphiques ACF et PACF.</figcaption></figure>

Les lignes bleues tracées représentent les limites de l'intervalle de confiance à 95%. Si une barre de corrélation dépasse ces limites, cela signifie que la corrélation est statistiquement significative, c'est-à-dire qu'elle est significativement différente de zéro.

En observant la présence d'une tendance et d'une composante saisonnière dans la série, on peut conclure qu'elle est non stationnaire. Pour éliminer la tendance, une différenciation d'ordre 1 (d = 1) est nécessaire, et pour supprimer la saisonnalité, une différenciation saisonnière d'ordre 1 (D = 1) est recommandée. Ainsi, il est approprié de transformer la série initiale en effectuant une double différenciation afin d'obtenir une série stationnaire. En examinant l'autocorrélation (ACF) et l'autocorrélation partielle (PACF), nous pouvons identifier les ordres des termes non saisonniers (p, q) et des termes saisonniers (P, Q) du modèle.

Modèle : SARIMA (1,1,1)(5,1,0)[3]
L'utilisation de la fonction "Arima" du package "forecast" permet d'estimer les coefficients du modèle SARIMA ainsi que la variance résiduelle, notée $\sigma_{\epsilon}$, en utilisant la méthode du maximum de vraisemblance.

Pour évaluer la significativité des coefficients du modèle, nous vérifions les p-values des tests de Student.

<figure><img src="image_49.png" /><figcaption>FIGURE 4 – Estimation des coefficients du modèle et leurs intervalles de confiance.</figcaption></figure>

On remarque que les coefficients estimés sont significatifs au seuil 5% ($p$-value < 0.05).

Enfin, nous validons la normalité des résidus en utilisant le test de Blancheur des résidus (Test de Portmanteau).

Le graphique des p-valeurs du test de Ljung et Box Fig.5 (graphique inférieur) montre que l'hypothèse de blancheur des résidus est acceptée, car les p-values associées aux différents décalages ne sont pas statistiquement significatives ($p$-value > 0.05). Cela confirme que les résidus ne présentent pas de corrélation sérielle significative.

L'autocorrélogramme indique que les résidus présentent une forte corrélation pour un décalage de $h = 0$, mais une faible corrélation avec les décalages ultérieurs ($h > 0$). Cela suggère qu'il ne reste pas de structure de corrélation significative dans les résidus après avoir tenu compte du modèle.

Le graphique 6 montre clairement que l'ajustement du modèle proposé est insuffisant pour capturer la structure temporelle des données de manière satisfaisante.
<figure><img src="image_50.png" /><figcaption>FIGURE 5 – Analyse des résidus standardisés.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>FIGURE 6 – série brute vs série ajustée.</figcaption></figure>
## Modèle de Spline cubique de lissage :

C'est une méthode utilisée pour ajuster une courbe ou une fonction à un ensemble de données. Il s'agit d'une approche non paramétrique qui permet de modéliser des relations complexes entre les variables sans imposer de structure spécifique.

Le modèle de spline cubique de lissage divise la courbe en segments plus petits et ajuste un polynôme cubique à chaque segment. Les polynômes cubiques sont choisis de manière à ce que la courbe résultante soit lisse et régulière. Le modèle utilise des noeuds, qui sont des points où les polynômes cubiques se rejoignent, pour contrôler la flexibilité de l'ajustement. Plus il y a de noeuds, plus le modèle est flexible et capable de capturer des variations complexes dans les données. En revanche, trop de noeuds peuvent entraîner un surajustement.

L'estimation des paramètres dans le modèle de spline cubique de lissage est réalisée en minimisant une fonction d'erreur pénalisée, généralement la somme des carrés des résidus. Cette approche permet de trouver les valeurs optimales des coefficients des polynômes cubiques pour obtenir la meilleure correspondance entre le modèle et les données.

$$S(f) = \sum_{i=1}^{n} (y_i - f(x_i))^2 + \lambda \int f''(x^2) dx$$

Le choix de paramètre λ est effectué par la méthode de Validation croisée généralisée.

<figure><img src="image_52.png" /><figcaption>FIGURE 7 – Ajustement du modèle de Spline cubique sur les données de précipitations.</figcaption></figure>

Nous pouvons clairement constaté que la courbe lissée ne parvient pas à capturer les variations brusques présentes dans les données, cela indique que le modèle ne correspond pas adéquatement aux données.
# F Détection de rupture

Dans cette section de l'annexe, nous présentons les différentes méthodes et fonctions utilisées pour détecter les ruptures dans les données de précipitations à Mayotte.

## Détection de ruptures en moyenne :

Nous utilisons la fonction **cpt.mean** du package **changepoint** avec les méthodes suivantes :
1. **PELT**: C'est une méthode de détection de rupture qui offre une approche exacte pour identifier les points de rupture dans une série de données. Elle est robuste aux valeurs aberrantes et peut s'adapter à différentes structures de données.

<figure><img src="image_53.png" /><figcaption>FIGURE 8 – Détection de ruptures en moyenne par la méthode PELT.</figcaption></figure>

2. *BinSeg* : Elle utilise une approche itérative pour trouver la segmentation optimale qui maximise la similarité des valeurs à l'intérieur des segments et minimise la différence entre les segments adjacents. BinSeg est connue pour sa simplicité et sa rapidité de calcul.
<figure><img src="image_54.png" /><figcaption>FIGURE 9 – Détection de ruptures en moyenne par la méthode BinSeg.</figcaption></figure>

Détection de ruptures en moyenne-variance :

Dans ce cas, nous utilisons la fonction cpt.meanvar du package changepoint.

<figure><img src="image_55.png" /><figcaption>FIGURE 10 – Détection de ruptures en moyenne-variance par la méthode PELT.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>FIGURE 11 – Détection de ruptures en moyenne-variance par la méthode BinSeg.</figcaption></figure>

Comme prévu, les résultats obtenus à partir des fonctions mentionnées ci-dessus varient considérablement en fonction des paramètres choisis.

Pour sélectionner l’algorithme de segmentation et les paramètres associés, nous nous baserons sur deux critères :

1. Critères théoriques : Nous évaluerons quel modèle est le plus approprié pour nos données en prenant en compte des critères tels que la constance de la variance, la normalité des résidus et la capacité à identifier des segments de longueurs variables, entre autres.

2. Critères empiriques : Nous choisirons la méthode qui produit les résultats les plus exploitables sur des données réelles.
Bibliographie

[1] D. Bernoulli, "Essai d'une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l'inoculation pour la prévenir," *Histoire de l'Acad., Roy. Sci. (Paris) avec Mem*, pp. 1-45, 1760.

[2] W. O. Kermack and A. G. McKendrick, "A contribution to the mathematical theory of epidemics," *Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character*, vol. 115, no. 772, pp. 700-721, 1927.

[3] D. G. Kendall, "Deterministic and stochastic epidemics in closed populations," in *Proceedings of the third Berkeley symposium on mathematical statistics and probability*, vol. 4, pp. 149-165, University of California Press Berkeley, 1956.

[4] R. Anderson, "May rm infectious diseases of humans," *Dynamics and Control*, vol. 281992, 1991.

[5] N. T. Bailey *et al.*, *The mathematical theory of infectious diseases and its applications*. Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE., 1975.

[6] F. Brauer, C. Castillo-Chavez, and C. Castillo-Chavez, *Mathematical models in population biology and epidemiology*, vol. 2. Springer, 2012.

[7] https://fr.wikipedia.org/wiki/Mayotte.

[8] https://fr.wikipedia.org/wiki/D%C3%A9mographie_de_Mayotte#cite note-diploweb-7.

[9] https://www.insee.fr/fr/statistiques/3713016.

[10] https://www.mayotte.ars.sante.fr/les-bornes-fontaines-monetiques.

[11] "Fièvre typhoïde à mayotte. point épidémiologique," *Point de situation au 18/08/2022*.

[12] R. M. Anderson and R. M. May, *Infectious diseases of humans : dynamics and control*. Oxford university press, 1991.

[13] M. H. Davis, "Piecewise-deterministic markov processes : A general class of non-diffusion stochastic models," *Journal of the Royal Statistical Society : Series B (Methodological)*, vol. 46, no. 3, pp. 353-376, 1984.

[14] K. Pakdaman, M. Thieullen, and G. Wainrib, "Fluid limit theorems for stochastic hybrid systems with application to neuron models," *Advances in Applied Probability*, vol. 42, no. 3, pp. 761-794, 2010.

[15] A. Crudu, A. Debussche, A. Muller, and O. Radulescu, "Convergence of stochastic gene networks to hybrid piecewise deterministic processes," 2012.
[16] M. Costa, "A piecewise deterministic model for a prey-predator community," 2016.

[17] F. Dufour, B. de Saporta, and H. Zhang, *Numerical methods for simulation and optimization of piecewise deterministic Markov processes: application to reliability*. John Wiley & Sons, 2015.

[18] S. Lessard, *Processus stochastiques*. ellipses, 2014.

[19] B. D. Saporta, *Modélisation en écologie, poly de cours, Université de Montpellier,*. 2019-2020.

[20] L. R. Rabiner, "A tutorial on hidden markov models and selected applications in speech recognition," *Proceedings of the IEEE*, vol. 77, no. 2, pp. 257-286, 1989.

[21] K. Hallouli, L. Likforman-Sulem, and M. Sigelle, "Reconnaissance de caractères manuscrits par réseaux bayésiens dynamiques," in *Conférence Internationale Francophone sur l'Écrit et le Document (CIFED 04)*, 2004.

[22] R. Durbin, S. R. Eddy, A. Krogh, and G. Mitchison, *Biological sequence analysis: probabilistic models of proteins and nucleic acids*. Cambridge university press, 1998.

[23] R. Bhar and S. Hamori, *Hidden Markov models: Applications to financial economics*, vol. 40. Springer Science & Business Media, 2004.

[24] P. A. Devijver, "Baum's forward-backward algorithm revisited," *Pattern Recognition Letters*, vol. 3, no. 6, pp. 369-373, 1985.

[25] L. E. Baum and T. Petrie, "Statistical inference for probabilistic functions of finite state markov chains," *The annals of mathematical statistics*, vol. 37, no. 6, pp. 1554-1563, 1966.

[26] D. T. Lauria, B. Maskery, C. Poulos, and D. Whittington, "An optimization model for reducing typhoid cases in developing countries without increasing public spending," *Vaccine*, vol. 27, no. 10, pp. 1609-1621, 2009.

[27] D. Kalajdzievska and M. Y. Li, "Modeling the effects of carriers on transmission dynamics of infectious diseases," *Math. Biosci. Eng.*, vol. 8, no. 3, pp. 711-722, 2011.

[28] S. Mushayabasa, "Impact of vaccines of controlling typhoid fever in kassana-nankana district of upper east region of ghana : insight from a mathematical model," 2011.

[29] H. Abboubakar and R. Racke, "Mathematical modelling and optimal control of typhoid fever," 2019.

[30] V. E. Pitzer, C. C. Bowles, S. Baker, G. Kang, V. Balaji, J. J. Farrar, and B. T. Grenfell, "Predicting the impact of vaccination on the transmission dynamics of typhoid in south asia : a mathematical modeling study," *PLoS neglected tropical diseases*, vol. 8, no. 1, p. e2642, 2014.

[31] https://www.biutils.ch/protocoles/13-analyse-de-l?eau.

[32] D. Sur, L. von Seidlein, B. Manna, S. Dutta, A. K. Deb, B. L. Sarkar, S. Kanungo, J. L. Deen, M. Ali, D. R. Kim, et al., "The malaria and typhoid fever burden in the slums of kolkata, india : data from a prospective community-based study," *Transactions of the Royal Society of Tropical Medicine and Hygiene*, vol. 100, no. 8, pp. 725-733, 2006.

[33] M. I. Khan, A. Pach 3rd, G. M. Khan, D. Bajracharya, S. Sahastrabuddhe, W. Bhutta, R. Tahir, S. Soofi, C. B. Thapa, N. Joshi, et al., "Typhoid vaccine introduction : An evidence-based pilot implementation project in nepal and pakistan," *Vaccine*, vol. 33, pp. C62-C67, 2015.
[34] J. A. Crump, F. G. Youssef, S. P. Luby, M. O. Wasfy, J. M. Rangel, M. Taalat, S. A. Oun, and F. J. Mahoney, "Estimating the incidence of typhoid fever and other febrile illnesses in developing countries," *Emerging infectious diseases*, vol. 9, no. 5, pp. 539-544, 2003.

[35] K. P. Neil, S. V. Sodha, L. Lukwago, S. O-tipo, M. Mikoleit, S. D. Simington, P. Mukobi, S. Balinandi, S. Majalija, J. Ayers, et al., "A large outbreak of typhoid fever associated with a high rate of intestinal perforation in kasese district, uganda, 2008-2009," *Clinical infectious diseases*, vol. 54, no. 8, pp. 1091-1099, 2012.

[36] A. Akullian, E. Ng'eno, A. I. Matheson, L. Cosmas, D. Macharia, B. Fields, G. Bi-gogo, M. Mugoh, G. John-Stewart, J. L. Walson, et al., "Environmental transmission of typhoid fever in an urban slum," *PLoS neglected tropical diseases*, vol. 9, no. 12, p. e0004212, 2015.

[37] A. Gebreyesus, K. Adane, L. Negash, T. Asmelash, S. Belay, M. Alemu, and M. Saravanan, "Prevalence of salmonella typhi and intestinal parasites among food handlers in mekelle university student cafeteria, mekelle, ethiopia," *Food Control*, vol. 44, pp. 45-48, 2014.

[38] J. John, C. J. Van Aart, and N. C. Grassly, "The burden of typhoid and paratyphoid in india : systematic review and meta-analysis," *PLoS neglected tropical diseases*, vol. 10, no. 4, p. e0004616, 2016.

[39] P. Roumagnac, F.-X. Weill, C. Dolecek, S. Baker, S. Brisse, N. T. Chinh, T. A. H. Le, C. J. Acosta, J. Farrar, G. Dougan, et al., "Evolutionary history of salmonella typhi," *Science*, vol. 314, no. 5803, pp. 1301-1304, 2006.

[40] V. E. Pitzer, N. A. Feasey, C. Msefula, J. Mallewa, N. Kennedy, Q. Dube, B. Denis, M. A. Gordon, and R. S. Heyderman, "Mathematical modeling to assess the drivers of the recent emergence of typhoid fever in blantyre, malawi," *Clinical Infectious Diseases*, vol. 61, no. suppl_4, pp. S251-S258, 2015.

[41] O. Diekmann, J. A. P. Heesterbeek, and J. A. Metz, "On the definition and the computation of the basic reproduction ratio r 0 in models for infectious diseases in heterogeneous populations," *Journal of mathematical biology*, vol. 28, pp. 365-382, 1990.

[42] P. Van den Driessche and J. Watmough, "Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission," *Mathematical biosciences*, vol. 180, no. 1-2, pp. 29-48, 2002.

[43] C. T. Codeço, S. Lele, M. Pascual, M. Bouma, and A. I. Ko, "A stochastic model for ecological systems with strong nonlinear response to environmental drivers : application to two water-borne diseases," *Journal of The Royal Society Interface*, vol. 5, no. 19, pp. 247-252, 2008.

[44] A. S. Novozhilov, G. P. Karev, and E. V. Koonin, "Biological applications of the theory of birth-and-death processes," *Briefings in bioinformatics*, vol. 7, no. 1, pp. 70-85, 2006.

[45] S. Nee, "Birth-death models in macroevolution," *Annu. Rev. Ecol. Evol. Syst.*, vol. 37, pp. 1-17, 2006.

[46] J. L. Thorne, H. Kishino, and J. Felsenstein, "An evolutionary model for maximum likelihood alignment of dna sequences," *Journal of Molecular Evolution*, vol. 33, no. 2, pp. 114-124, 1991.
[47] F. W. Crawford, V. N. Minin, and M. A. Suchard, "Estimation for general birth-death processes," *Journal of the American Statistical Association*, vol. 109, no. 506, pp. 730-747, 2014.

[48] P. Moran, "The estimation of the parameters of a birth and death process," *Journal of the Royal Statistical Society : Series B (Methodological)*, vol. 15, no. 2, pp. 241-245, 1953.

[49] J. F. Reynolds, "On estimating the parameters of a birth-death process," *Australian Journal of Statistics*, vol. 15, no. 1, pp. 35-43, 1973.

[50] N. Keiding, "Maximum likelihood estimation in the birth-and-death process," *The Annals of Statistics*, vol. 3, no. 2, pp. 363-372, 1975.

[51] T. Perkins, "Maximum likelihood trajectories for continuous-time markov chains," in *Advances in Neural Information Processing Systems* (Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, eds.), vol. 22, Curran Associates, Inc., 2009.

[52] E. R. Immel, *Problems of estimation and hypothesis testing in connection with birth-and-death stochastic processes*. PhD thesis, Doctoral Thesis, University of California, Los Angeles. Résumé : Ann. Math ?, 1951.

[53] F. W. Crawford and M. A. Suchard, "Transition probabilities for general birth-death processes with applications in ecology, genetics, and evolution," *Journal of mathematical biology*, vol. 65, no. 3, pp. 553-580, 2012.

[54] J. Xu, P. Guttorp, M. Kato-Maeda, and V. N. Minin, "Likelihood-based inference for discretely observed birth-death-shift processes, with applications to evolution of mobile genetic elements," *Biometrics*, vol. 71, no. 4, pp. 1009-1021, 2015.

[55] A. C. Davison, S. Hautphenne, and A. Kraus, "Parameter estimation for discretely observed linear birth-and-death processes," *Biometrics*, vol. 77, no. 1, pp. 186-196, 2021.

[56] I. Bouzalmat, B. de Saporta, and S. M. Manou-Abi, "Parameter estimation for a hidden linear birth and death process with immigration," *arXiv preprint arXiv:2303.00531*, 2023.

[57] J. Fritz, J. LaSalle, and L. Sirovich, *Partial differential equations*. Applied Mathematical Sciences, 1977.

[58] S. Tavaré, "The linear birth-death process : an inferential retrospective," *Advances in Applied Probability*, vol. 50, no. A, pp. 253-269, 2018.

[59] N. T. Bailey, *The elements of stochastic processes with applications to the natural sciences*, vol. 25. John Wiley & Sons, 1991.

[60] S. M. Ross, *Introduction to probability models*. Academic press, 2014.

[61] R. M. Corless, G. H. Gonnet, and D. Hare, "On the lambert w function," *Advances in Computational Mathematics*, vol. 5, pp. 329-359, 1996.

[62] R. M. Corless, G. H. Gonnet, D. E. Hare, D. J. Jeffrey, and D. Knuth, "Lambert's w function in maple," *Maple Technical Newsletter*, vol. 9, no. 1, pp. 12-22, 1993.

[63] I. Teodorescu, "Maximum likelihood estimation for markov chains." *arXiv preprint arXiv:0905.4131*, 2009.

[64] G. M. Goerg, *LambertW : Probabilistic Models to Analyze and Gaussianize Heavy-Tailed, Skewed Data*, 2020. R package version 0.6.6.

[65] M. L. Himmelmann, "Package ?hmm ?," 2022.
[66] K. Athreya and C. Fuh, "Bootstrapping markov chains : countable case," *Journal of Statistical Planning and Inference*, vol. 33, no. 3, pp. 311-331, 1992.

[67] S. B. Y. M. T. D. Y. I. C. V. J. T. G. R. J. B. D. E. S. M. Giorgio Alfredo Spedicato, Tae Seung Kang, "Markovchain : Easy handling discrete time markov chains. 2022," 2022.

[68] C. M. Bishop and N. M. Nasrabadi, *Pattern recognition and machine learning*, vol. 4. Springer, 2006.

[69] K. P. Murphy, *Machine learning : a probabilistic perspective*. MIT press, 2012.

[70] H. Abboubakar and R. Racke, "Mathematical modeling, forecasting, and optimal control of typhoid fever transmission dynamics," *Chaos, Solitons & Fractals*, vol. 149, p. 111074, 2021.

[71] https://www.inrs.fr/publications/bdd/eficatt/fiche.html?refINRS=EFICATT_Fièvretyphoïde.

[72] S. P. Meyn and R. L. Tweedie, "Stability of markovian processes iii : Foster-lyapunov criteria for continuous-time processes," *Advances in Applied Probability*, vol. 25, no. 3, pp. 518-548, 1993.

[73] P. Kumar, H. Gauba, P. P. Roy, and D. P. Dogra, "Coupled hmm-based multi-sensor data fusion for sign language recognition," *Pattern Recognition Letters*, vol. 86, pp. 1-8, 2017.

[74] S. Ghosh, J. Li, L. Cao, and K. Ramamohanarao, "Septic shock prediction for icu patients via coupled hmm walking on sequential contrast patterns," *Journal of bio-medical informatics*, vol. 66, pp. 19-31, 2017.

[75] R. Sarno and K. R. Sungkono, "Coupled hidden markov model for process discovery of non-free choice and invisible prime tasks," *Procedia Computer Science*, vol. 124, pp. 134-141, 2017.

[76] M. Brand, "Coupled hidden markov models for modeling interacting processes," tech. rep., Citeseer, 1997.

[77] J. A. Bilmes et al., "A gentle tutorial of the em algorithm and its application to parameter estimation for gaussian mixture and hidden markov models," *International computer science institute*, vol. 4, no. 510, p. 126, 1998.

[78] R. Calaway, R. Analytics, S. Weston, D. Tenenbaum, and M. R. Calaway, "Package ?doparallel?," 2015.

[79] https://www.historique-meteo.net/afrique/mayotte/.

[80] E. S. Page, "Continuous inspection schemes," *Biometrika*, vol. 41, no. 1/2, pp. 100-115, 1954.

[81] A. Wald, "Statistical decision functions which minimize the maximum risk," *Annals of Mathematics*, pp. 265-280, 1945.

[82] D. Egea-Roca, G. Seco-Granados, and J. A. López-Salcedo, "Comprehensive overview of quickest detection theory and its application to gnss threat detection," *Gyroscopy and Navigation*, vol. 8, pp. 1-14, 2017.

[83] R. Killick and I. Eckley, "changepoint : An r package for changepoint analysis," *Journal of statistical software*, vol. 58, no. 3, pp. 1-19, 2014.
[84] M. D. Sobsey, C. E. Stauber, L. M. Casanova, J. M. Brown, and M. A. Elliott, "Point of use household drinking water filtration : a practical, effective solution for providing sustained access to safe drinking water in the developing world," *Environmental science & technology*, vol. 42, no. 12, pp. 4261-4267, 2008.

[85] K. Cann, D. R. Thomas, R. Salmon, A. Wyn-Jones, and D. Kay, "Extreme water-related weather events and waterborne disease," *Epidemiology & Infection*, vol. 141, no. 4, pp. 671-686, 2013.

[86] O. Oguntoke, O. Aboderin, and A. Bankole, "Association of water-borne diseases morbidity pattern and water quality in parts of ibadan city, nigeria," *Tanzania Journal of health research*, vol. 11, no. 4, 2009.

[87] M. Hatta and H. L. Smits, "Detection of salmonella typhi by nested polymerase chain reaction in blood, urine, and stool samples," *The American journal of tropical medicine and hygiene*, vol. 76, no. 1, pp. 139-143, 2007.

[88] G. C. Buckle, C. L. F. Walker, and R. E. Black, "Typhoid fever and paratyphoid fever : Systematic review to estimate global morbidity and mortality for 2010.," *Journal of global health*, vol. 2, no. 1, 2012.

[89] I. Gantois, R. Ducatelle, F. Pasmans, F. Haesebrouck, R. Gast, T. J. Humphrey, and F. Van Immerseel, "Mechanisms of egg contamination by salmonella enteritidis," *FEMS microbiology reviews*, vol. 33, no. 4, pp. 718-738, 2009.

[90] T. Akritidou, S. Akkermans, S. Gaspari, N. D. Azraini, C. Smet, T. Van de Wiele, and J. F. Van Impe, "Effect of gastric ph and bile acids on the survival of listeria monocytogenes and salmonella typhimurium during simulated gastrointestinal digestion," *Innovative Food Science & Emerging Technologies*, vol. 82, p. 103161, 2022.

[91] R. J. Hyndman and G. Athanasopoulos, *Forecasting : principles and practice*. OTexts, 2018.