[PYTHON] J'ai essayé de prédire le nombre de personnes infectées au niveau national de la nouvelle corona avec un modèle mathématique

Aperçu

J'ai trouvé l'article de The Lancet [^ 1], qui prédit le nombre de personnes infectées avec un modèle mathématique, alors j'ai implémenté un modèle mathématique pour prédire le nombre de personnes infectées au Japon et je l'ai comparé à d'autres maladies infectieuses.

Résumé

Le nombre de personnes infectées a été prédit par le modèle SEIR. ―― Le nombre de reproduction de base de la nouvelle corona est d'environ 2,9, ce qui est un peu plus élevé que celui de la grippe.
Si le contrôle de l'infection n'est pas efficace, le nombre maximal de personnes infectées au pays est d'environ GW. ――Le nombre de décès domestiques l'année prochaine sera supérieur à 10 000 (lorsque le taux de létalité est de 0,01%).
Les restrictions d'immigration pour les étrangers n'ont aucun sens.

table des matières

[Objectif de cette fois](# 1-Objectif de cette fois)
[Modèle mathématique des maladies infectieuses](# 2-Modèle mathématique des maladies infectieuses)
[Estimation des paramètres](# 3-Estimation des paramètres)
[Implémenté et illustré en python](# 4-Implémenté et illustré en python)
[Supplément](# 5-Supplément)

1. Objet de cette période

Prédire le nombre de futures infections domestiques du nouveau virus corona (COVID-19, ci-après dénommé le nouveau corona) à l'aide d'un modèle mathématique.
Comparez les nombres de reproduction de base de la nouvelle couronne avec d'autres maladies infectieuses.
Mise en œuvre de calculs scientifiques avec Scipy.

2. Modèle mathématique des maladies infectieuses

"Formulez le processus de propagation d'une maladie infectieuse, la durée pendant laquelle une personne infectée se développe et devient grave." [^ 2]

Dans le modèle actuariel des maladies infectieuses, la population est divisée selon le stade de l'infection, et l'augmentation ou la diminution de chaque groupe est exprimée par une équation différentielle. Je vais vous expliquer le modèle SEIR simple et la version étendue utilisée cette fois.

2.1. Modèle SEIR

$ S $: Nombre de personnes sans immunité (Susceptible)
$ E $: Nombre de personnes dans la période post-exposition (Exposées)
$ I $: Nombre de personnes pendant la période infectieuse (infectieuse)
$ R $: Nombre de personnes qui se sont rétablies d'une infection et ont acquis une immunité (récupérées)
$ \ beta $: Taux d'infection, $ \ gamma $: Taux de récupération (isolement), $ \ sigma $: Probabilité d'infectivité, $ N $: Population totale

Ensuite, les équations différentielles simultanées suivantes sont obtenues.

\begin{align}
\frac{dS(t)}{dt} &= -\beta \frac{I(t)}{N}S(t) \\

\frac{dE(t)}{dt} &= \beta \frac{I(t)}{N}S(t) - \sigma E(t) \\

\frac{dI(t)}{dt} &= \sigma E(t) - \gamma I(t) \\

\frac{dR(t)}{dt} &= \gamma I(t) \\

N &= S(t) + E(t) + I(t) + R(t)
\end{align}

Cependant, il ignore le taux de natalité et le taux de mortalité, et dit qu'il ne sera plus jamais infecté car il se rétablira et acquiert une immunité après le début.

Si vous donnez une valeur initiale appropriée et que vous la résolvez, vous pouvez prédire le nombre de personnes infectées à l'avenir. (Je mettrai le chiffre [^ 3] pour référence)

Pour une explication et une implémentation détaillées, reportez-vous à l'article de qiita [^ 4].

Ici, $ R_0 = \ beta / \ gamma $ est appelé le nombre de reproduction de base, et "(tous les individus sont initialement sensibles. Cela signifie "le nombre de personnes infectées secondaires produites par personne infectée (en état d'avoir)". [^ 5] ** Le fait est que plus $ R_0 $ est élevé, plus l'infectivité est grande. ** **

Vous pouvez également utiliser $ R_0 $ pour calculer dans quelle mesure la vaccination généralisée peut prévenir une épidémie.

En plus du SEIR, divers modèles peuvent être créés en subdivisant la population entière par lieu de résidence et par âge, et en reflétant les caractéristiques des maladies infectieuses dans des équations différentielles. Par exemple, il y a SIR qui ne considère pas E, SEIS qui ne considère pas l'acquisition d'immunité et MSEIR qui considère l'immunité passive [^ 6].

2.2. Modèle SEIRD étendu

Le SEIR ci-dessus étant trop simple, nous l'étendrons cette fois à un modèle qui considère les entrées et sorties de population et le taux de mortalité avec les pays étrangers en se référant à l'article [^ 1].

$ R_0 : Numéro de reproduction de base ( R_0 = \ beta / \ gamma $)
$ T_E : Période moyenne post-exposition (jours) ( \ gamma ^ {-1} = T_I $)
$ T_I : Période d'infectivité moyenne (jours) ( \ sigma ^ {-1} = T_E $)
$ N_ {J, I} $: nombre de voyageurs du Japon à l'étranger (personnes / jour)
$ N_ {I, J} $: nombre de voyageurs d'outre-mer vers le Japon (personnes / jour)
$ f $: létalité
$ D (t) $: Nombre de décès (personnes)

Ensuite, les équations différentielles simultanées suivantes sont également obtenues.

\begin{align}
\frac{dS(t)}{dt} &= - \frac{R_0}{D_I} \frac{I(t)}{N(t)} S(t) - N_{J,I} + 
 \left(1 - \frac{E(t) + I(t)}{N(t)} \right)N_{I,J} \\

\frac{dE(t)}{dt} &= \frac{R_0}{D_I} \frac{I(t)}{N(t)} S(t) - \frac{E(t)}{T_E} + \frac{E(t) + I(t)}{N(t)} N_{I,J}\\

\frac{dI(t)}{dt} &= \frac{E(t)}{T_E} - \frac{I(t)}{T_I}\\

\frac{dR(t)}{dt} &= (1-f) \frac{I(t)}{T_I}\\

\frac{dD(t)}{dt} &= f \frac{I(t)}{T_I} \\

N(t) &= S(t) + E(t) + I(t) + R(t)
\end{align}

Cependant, en ce qui concerne $ N_ {J, I} $ et $ N_ {I, J} $, on suppose que tous ceux qui partent du Japon à l'étranger appartiennent à $ S $, et ceux qui viennent au Japon de l'étranger appartiennent à $ S, E $. Je vais.

Comme dans l'article [^ 1], les paramètres autres que $ R_0 $ sont remplacés par des valeurs appropriées. $ T_I et T_E $ sont mis à $ T_E = 6,0 et T_I = 2,4 $ en référence à SARS Corona. $ N_ {J, I}, N_ {I, J} $ est défini sur $ N_ {J, I} = 40000, N_ {I, J} = 70000 $ en référence au nombre de voyageurs en février 2019.

3. Estimation des paramètres

Le paramètre inconnu du modèle est $ R_0 $. Il existe diverses méthodes d'estimation telles que la méthode du carré minimum, l'estimation de probabilité maximale, l'estimation MAP, l'estimation bayésienne [^ 7], mais cette fois, la vraisemblance est calculée comme le processus de Poisson non stationnaire [^ 8] et $ R_0 $ est calculé par l'estimation de probabilité maximale. ..

3.1. Estimation la plus probable

Tout le monde aime l'estimation la plus probable. Tout d'abord, exprimez la probabilité (probabilité) des données observées en fonction de $ R_0 $, et trouvez $ R_0 $ qui la maximise. Lorsque le nombre de personnes infectées est modélisé comme un processus de Poisson non stationnaire en référence à l'article [^ 1] et que la probabilité est calculée,

L(R_0) = \prod_{d=D_s}^{D_e} \frac{e^{-\lambda_d}\lambda_d^{x_d}}{x_d!}

Cependant, $ D_s $: premier jour d'observation, $ D_s $: dernier jour d'observation, $ x_d $: nombre de personnes infectées à $ d $

\lambda (t, R_0) = E(t) + I(t) \\
\lambda _d (R_0) = \int_{d-1}^{d} \lambda (u, R_0) du

ça ira. $ \ lambda $ représente le nombre moyen de personnes infectées lorsqu'elles suivent la distribution de Poisson.

En excluant les termes non liés à $ R_0 $ de la vraisemblance logarithmique, $ R_0 $ peut être estimé ci-dessous.

\newcommand{\argmin}{\mathop{\rm arg~min}\limits}

\argmin_{R_0 > 0} - \log L(R_0) = \argmin_{R_0 > 0} \sum_{d=D_s}^{D_e} (\lambda_d(R_0) - x_d \log \lambda_d(R_0))

4. Implémenté et illustré en python

Implémentation de la classe SEIR

import scipy
from scipy.integrate import odeint
from scipy.optimize import minimize
from sklearn.metrics import mean_absolute_error, mean_squared_error


class SIR(object):

    def __init__(self, r_0, t_i, dt, init_S, init_I, init_R):
        self.r_0 = r_0
        self.t_i = t_i
        self.dt = dt
        self.init_state_list = [init_S, init_I, init_R]

    def _get_param(self, params, key, default, t=None):
        """Correspond à des paramètres variant dans le temps
        """
        if isinstance(params, dict):
            if key in list(params.keys()):
                param = params[key]
                if isinstance(param, list):
                    return param[np.clip(int(t / self.dt), 0, len(param)-1)]
                elif isinstance(param, np.ndarray):
                    return param
                else:
                    return default
            else:
                return default
        else:
            return default

    def _ode(self, state_list, t=None, params=None):
        """Définir des équations différentielles simultanées
        """
        r_0 = self._get_param(params, 'r_0', self.r_0, t=t)
        t_i = self._get_param(params, 't_i', self.t_i, t=t)
        S, I, R = state_list
        N = S + I + R
        dstate_dt = list()
        dstate_dt[0] = - (r_0 / t_i) * (I / N) * S
        dstate_dt[1] = (r_0 / t_i) * (I / N) * S - I / t_i
        dstate_dt[2] = I / t_i
        return dstate_dt

    def solve_ode(self, len_days=365, params=None):
        """Résoudre des équations différentielles
        """
        t = np.linspace(0, len_days, int(len_days / self.dt), endpoint=False)
        args = (params,) if params else ()
        return odeint(self._ode, self.init_state_list, t, args=args)


class CustomizedSEIRD(SIR):

    def __init__(self, r_0=None, t_e=6.0, t_i=2.4, n_i_j=70000, n_j_i=40000, f=0.0001, dt=1,
            init_S=126800000, init_E=0, init_I=1, init_R=0, init_D=0):
        self.r_0 = r_0
        self.t_e = t_e
        self.t_i = t_i
        self.n_i_j = n_i_j
        self.n_j_i = n_j_i
        self.f = f
        self.dt = dt
        self.init_state_list = [init_S, init_E, init_I, init_R, init_D]

    def _ode(self, state_list, t=None, params=None):
        """Définir des équations différentielles simultanées
        """
        r_0 = self._get_param(params, 'r_0', self.r_0, t=t)
        t_e = self._get_param(params, 't_e', self.t_e, t=t)
        t_i = self._get_param(params, 't_i', self.t_i, t=t)
        n_i_j = self._get_param(params, 'n_i_j', self.n_i_j, t=t)
        n_j_i = self._get_param(params, 'n_j_i', self.n_j_i, t=t)
        f = self._get_param(params, 'f', self.f)
        S, E, I, R, D = state_list
        N = S + E + I + R
        dstate_dt = list()
        dstate_dt.append(- (r_0 / t_i) * (I / N) * S - n_j_i + (1 - ((E + I) / N)) * n_i_j)
        dstate_dt.append((r_0 / t_i) * (I / N) * S - E / t_e + ((E + I) / N) * n_i_j)
        dstate_dt.append(E / t_e - I / t_i)
        dstate_dt.append((1 - f) * I / t_i)
        dstate_dt.append(f * I / t_i)
        return dstate_dt

    def _calc_neg_log_likelihood_r0(self, r_0, X):
        """Probabilité du journal(R_Seule la partie liée à 0)Calculer
        """
        solution = self.solve_ode(len_days=len(X), params=dict(r_0=r_0))
        lambda_arr = solution[int(1/self.dt)-1::int(1/self.dt), 2] # I
        return - np.sum(- lambda_arr + X * np.log(lambda_arr))
    
    def _calc_error(self, r_0, X, metric=mean_absolute_error):
        """Calculer l'erreur absolue moyenne et l'erreur quadratique moyenne
        """
        solution = self.solve_ode(len_days=len(X), params=dict(r_0=r_0))
        e_arr = solution[int(1/self.dt)-1::int(1/self.dt), 2]  # I
        return metric(X, e_arr)

    def exec_point_estimation(self, init_r_0, X, project='mle'):
        """Paramètres d'estimation ponctuelle
        """
        if project == 'mle':
            result = minimize(self._calc_neg_log_likelihood_r0, init_r_0, args=(X,), method='Nelder-Mead')
        elif project == 'lad':
            result = minimize(self._calc_error, init_r_0, args=(X, mean_absolute_error), method='Nelder-Mead')
        elif project == 'ls':
            result = minimize(self._calc_error, init_r_0, args=(X, mean_squared_error), method='Nelder-Mead')
        else:
            print(f'Invalid project: {project}')
            return None
        if self.r_0 is None:
            self.r_0 = result.x[0]
        return result

    def exec_map(self):
        """Estimation MAP
        """
        pass

    def exec_mcmc(self):
        """Méthode MCMC
        """
        pass

if __name__ == '__main__':
    # 1/16 à 3/Nombre de personnes infectées jusqu'à 8
    X = np.array([
        1, 1, 1, 1, 1, 1, 1, 1, 2, 3,  # 1/25
        4, 4, 7, 8, 14, 17, 20, 20, 20, 23,  # 2/4
        35, 45, 86, 90, 96, 161, 163, 203, 251, 259,  # 2/14
        338, 414, 520, 616, 705, 728, 743, 743, 838, 851, # 2/24
        862, 891, 919, 938, 947, 961, 980, 999, 1035, 1056, # 3/5
        1113, 1157, 1190 # 3/8
    ])
    model = CustomizedSEIRD()
    result = model.exec_point_estimation(init_r_0=1, X=X, project='mle')
    print(result.x[0])

2.8806152343750036

Le «nombre de personnes infectées» est difficile à gérer, et seuls ceux qui se rendent à l'hôpital sont effectivement comptés, et personne ne connaît le vrai $ x_d $ compte tenu des critères de diagnostic et de la spécificité de sensibilité du test. Hmm. Ici, il est estimé en fonction du nombre de personnes infectées [^ 9], y compris les navires de croisière.

Prédiction des futures infections domestiques par le modèle SEIR étendu (I est le nombre d'infections)

C'est un graphique lorsque $ R_0 $ estimé par la ligne bleue ne change pas dans le futur. Si nous continuons comme il est, le nombre de personnes infectées atteindra son maximum juste lorsque la GW sera terminée. C'est un graphique lorsque la ligne rouge et la ligne bleue changent $ R_0 $ après 3/9. Puisque $ R_0 $ est "le nombre de personnes infectées secondaires produites par personne infectée", il sera plus petit en raison des contre-mesures contre les infections et des restrictions de mouvement. Si $ R_0 $ peut être rendu suffisamment petit, il peut être contenu comme une ligne verte.

Le chiffre ci-dessous est un chiffre publié par le ministère de la Santé et du Bien-être social [^ 9], mais la forme du graphique du nombre de personnes infectées correspond exactement.

厚生省コロナ感染者数.png

De plus, par exemple, si vous modifiez $ N_ {i, j} $, vous pouvez calculer si le nombre de personnes infectées changera en raison des restrictions d'immigration des étrangers. En calculant avec le même modèle, même si $ N_ {i, j} = 0 $ après 3/9, si $ R_0 $ est le même, le nombre de personnes infectées ne changera pas. Cela signifie que les restrictions à l'immigration n'ont aucun sens (si ce modèle est correct).

Ensuite, comparons le nombre de reproduction de base et la létalité avec d'autres maladies infectieuses [^ 10]. (Le taux de létalité de la nouvelle couronne est d'environ 0,1% à 5% du taux de létalité de chaque pays)

Il semble que le nombre de reproduction de base et le taux de létalité de la nouvelle couronne ne soient pas significativement plus élevés que ceux des autres maladies infectieuses. (Pour le risque de maladies infectieuses, il est nécessaire de prendre en compte la gravité des symptômes et la présence ou l'absence de vaccins, et je pense que de faibles taux de reproduction de base et de faibles taux de létalité ne signifient pas qu'ils sont sûrs, mais les détails sont Demandez à quelqu'un qui résiste aux maladies infectieuses)

5. Supplément

Ministère de la santé, du travail et de la protection sociale

Société japonaise pour les maladies infectieuses

WHO

CDC

Site sympa où vous pouvez vérifier le nombre de personnes infectées (Johns Hopkins CSSE)

** Je ne suis ni un spécialiste des maladies infectieuses, ni un service de santé publique, ni un statisticien. La maladie infectieuse réelle n'est pas aussi simple que ce modèle, donc cette prédiction est probablement fausse. Nous ne sommes pas responsables des dommages causés par le contenu de cet article. ** **

[^ 2]: Arrêtez les maladies infectieuses avec un modèle mathématique

[^ 4]: Début du modèle mathématique des maladies infectieuses: aperçu du modèle SEIR par Python et introduction à l'estimation des paramètres

[^ 5]: Prédiction des épidémies de maladies infectieuses: problèmes quantitatifs dans les modèles mathématiques de maladies infectieuses

[^ 7]: Estimation la plus probable, estimation MAP, estimation Bayes apprise dans Fujii 4ème Dan

[^ 8]: modèle SIR et processus de Poisson non stationnaire

[^ 9]: À propos de la nouvelle infection à coronavirus