[PYTHON] [Suite] La nouvelle Corona est-elle vraiment une menace? Létalité calculée avec Stan

Le post précédent "La nouvelle Corona est-elle vraiment une menace? Je l'ai vérifié auprès de Stan (était)" était indigeste sur les deux points suivants:

  1. Après tout, je ne sais pas quel est le taux de létalité du COVID-19.
  2. Enfin, Stan ne fonctionne pas bien

Donc, cette fois, pour me venger, je vais essayer de trouver le "taux de létalité véritable" en utilisant Stan.

Je n'ai pas expliqué à plusieurs reprises ce que j'ai expliqué dans Article précédent, veuillez donc vous y référer également: bow:

Avertissement

Je ne suis pas un expert en modélisation bayésienne ou en Stan, il peut donc y avoir des erreurs et de meilleures techniques dans le contenu de cet article.

De plus, si vous trouvez quelque chose qui ne va pas, je vous serais reconnaissant de bien vouloir me le faire savoir. : arc:

supposition

Dans mon dernier article, j'ai soutenu que l'enquête sur l'infection [^ 1] sur le problème de Diamond Princess était basée sur l'hypothèse qu'elle était la plus précise au monde. Nous suivrons cette politique également cette fois.

Cependant, le nombre de cas est faible avec le Diamond Princess seul, et il n'est pas possible de connaître le taux de létalité par âge. Par conséquent, nous utiliserons les données [^ 2] centrées sur la lutte contre la maladie en Chine.

Comme les populations de ces deux données sont complètement différentes, je pose des hypothèses et les intègre d'une manière ou d'une autre. Ici, nous supposons que le nombre d'infections dans les données chinoises a été obtenu auprès de toutes les personnes réellement infectées avec une certaine probabilité de $ c $ quel que soit le groupe d'âge ou la mortalité.

Cette hypothèse est clairement inexacte, car les groupes d'âge varient dans les tendances symptomatiques et, par conséquent, la proportion d'individus infectés capturés peut varier. Par conséquent, les conclusions de cet article sont inexactes en raison de cette hypothèse. Cependant, selon les données de l'Institut national des maladies infectieuses [^ 1], la proportion de ceux qui présentent des symptômes (cas confirmés symptomatiques) par rapport à ceux qui n'en ont pas (cas confirmés asymptomatiques) est plutôt plus faible chez la jeune génération. Des rapports [^ 3] indiquent également que le taux de gravité des jeunes est étonnamment élevé. De là, je pense que c'est une hypothèse plus absurde que l'intuition.

La létalité dans le Diamond Princess combine les deux données en supposant qu'elle suit la létalité du groupe d'âge calculée en utilisant le taux de capture de $ c $ supposé ici.

En revanche, on suppose que la capture des morts se fait sans exception dans les statistiques chinoises.

Ensuite, le nombre de morts sur le Diamond Princess était de six. Le septième décès a été confirmé le 7 mars, mais il est possible qu'il n'ait pas été infecté à ce moment-là, car plus de deux semaines se sont écoulées depuis les données de l'Institut national des maladies infectieuses [^ 1]. Il est basé sur le fait qu'il est élevé, qu'il s'agisse d'un décès lié au COVID-19 ou non, n'a pas été annoncé et que le nombre de personnes infectées a augmenté depuis le moment des données au 7 mars.

Enfin, il existe des différences selon la race, des différences dans la répartition des personnes infectées (en fait, la répartition des personnes infectées en Chine et la répartition des personnes infectées dans le monde devraient être différentes) et des différences dans l'environnement médical (en particulier lorsque le nombre de patients augmente rapidement, un effondrement médical se produit. J'ai décidé de ne pas considérer (ce qui devrait augmenter le taux de létalité). Cependant, la différence de nombre de personnes infectées entre la Chine et le Diamond Princess est exceptionnellement prise en compte.

Les hypothèses peuvent être résumées comme suit.

modèle

Ci-dessous, $ N $ est le nombre de personnes infectées, $ D $ est le nombre de décès et $ p $ est le taux de létalité. Sauf indication contraire, l'exposant indique le groupe d'âge et l'indice indique l'emplacement ($ {\ rm dp} $ = Diamond Princess, $ {\ rm ch} $ = Chine). L'absence d'exposants représente tous les groupes d'âge.

Par exemple, $ N_ {{\ rm dp}} ^ {60-69} $ indique "le nombre de personnes infectées âgées de 60 à 69 ans identifiées sur le Diamond Princess". (En fait, je ne mentionnerai plus un groupe d'âge spécifique à partir de maintenant, mais j'écrirai $ N_ {{\ rm dp}} ^ {a} $ et j'écrirai "Infection du groupe d'âge $ a $ confirmée sur le Diamond Princess" Il est utilisé pour représenter le nombre de personnes.)

Supposons que la létalité estimée est de $ p $ si les hypothèses sont correctes. De plus, cette létalité estimée $ p $ est divisée par groupe d'âge et exprimée en $ p ^ {60-69} $ (sans indice).

Le nombre de décès parmi les personnes infectées de la tranche d'âge $ a $ en Chine suit une distribution binaire.

D_{{\rm ch}}^a \sim {\rm binomial}(N_{{\rm ch}}^a, p_{{\rm ch}}^a)

Peut être modélisé comme.

De même, le nombre de décès parmi les personnes infectées dans le Diamond Princess est également le même.

D_{{\rm dp}} \sim {\rm binomial}(N_{{\rm dp}}, p_{{\rm dp}})

Ce sera. Veuillez noter que nous ne modélisons pas par groupe d'âge car nous n'avons pas la répartition par âge des décès de Diamond Princess.

Reliez les deux avec le taux de capture $ c $ mentionné dans le chapitre Hypothèses (#Assumptions).

Puisqu'il n'y a aucune information sur le taux de capture, nous supposons une distribution uniforme.

c \sim {\rm uniform}(0, 1)

Nous avons supposé que $ c $ est constant indépendamment du groupe d'âge et de la létalité, donc il est indépendant de $ p_ {{\ rm ch}} ^ a $. Donc, la létalité estimée du groupe d'âge $ a $ est

p^a = c \times p_{{\rm ch}}^a

À tous âges

p = c \times {1 \over N_{{\rm ch}}} \sum_a{p_{{\rm ch}}^a \times N_{{\rm ch}}^a} = {1 \over N_{{\rm ch}}} \sum_a{p^a \times N_{{\rm ch}}^a}

Ce sera.

Nous avons également supposé que la létalité du Diamond Princess suit la létalité estimée ($ p_ {{\ rm dp}} ^ a = p ^ a $), donc $ p_ {{\ rm dp}} $ est

p_{{\rm dp}} = {1 \over N_{{\rm dp}}} \sum_a{p^a \times N_{{\rm dp}}^a}

Ce sera.

la mise en oeuvre

Maintenant que nous avons tous les acteurs, nous allons l'implémenter dans Stan.

Des données d'entrée

Les noms des variables sont «6x» (c'est-à-dire $ N {{\ rm ch}} ^ {60-69} $ = «N_ch_6x») pour les groupes d'âge 60-69 et «80_up» (c'est-à-dire $ N) pour les 80 ans et plus. C'est presque le même que le symbole [Model Chapter](# model) sauf qu'il est exprimé par {{\ rm ch}} ^ {80-} $ = N_ch_80_up).

Veuillez noter que les données Diamond Princess [^ 1] contiennent des données pour les années 90, mais les données chinoises [^ 2] concernent les personnes de plus de 80 ans.

data {
    //Nombre de morts sur le Diamond Princess
    int D_dp;

    //Nombre de personnes infectées par groupe d'âge sur le Diamond Princess
    int N_dp_0x;
    int N_dp_1x;
    int N_dp_2x;
    int N_dp_3x;
    int N_dp_4x;
    int N_dp_5x;
    int N_dp_6x;
    int N_dp_7x;
    int N_dp_8x;
    int N_dp_9x;

    //Nombre de décès par groupe d'âge en Chine
    int D_ch_0x;
    int D_ch_1x;
    int D_ch_2x;
    int D_ch_3x;
    int D_ch_4x;
    int D_ch_5x;
    int D_ch_6x;
    int D_ch_7x;
    int D_ch_80_up;

    //Nombre de personnes infectées par groupe d'âge en Chine
    int N_ch_0x;
    int N_ch_1x;
    int N_ch_2x;
    int N_ch_3x;
    int N_ch_4x;
    int N_ch_5x;
    int N_ch_6x;
    int N_ch_7x;
    int N_ch_80_up;
}

Étant donné que Stan peut définir les données traitées pour plus de commodité, nous calculerons les valeurs pour tous les groupes d'âge.

transformed data {
    //Nombre de personnes infectées sur le Diamond Princess
    int N_dp;

    //Nombre de personnes infectées en Chine
    int N_ch;

    N_dp = N_dp_0x + N_dp_1x + N_dp_2x + N_dp_3x + N_dp_4x + N_dp_5x + N_dp_6x + N_dp_7x + N_dp_8x + N_dp_9x;
    N_ch = N_ch_0x + N_ch_1x + N_ch_2x + N_ch_3x + N_ch_4x + N_ch_5x + N_ch_6x + N_ch_7x + N_ch_80_up;
}

Paramètres cibles estimés

Spécifiez les paramètres à estimer.

parameters {
    //Taux de capture c
    real<lower=0, upper=1> c;

    //Létalité en Chine
    real<lower=0, upper=1> p_ch_0x;
    real<lower=0, upper=1> p_ch_1x;
    real<lower=0, upper=1> p_ch_2x;
    real<lower=0, upper=1> p_ch_3x;
    real<lower=0, upper=1> p_ch_4x;
    real<lower=0, upper=1> p_ch_5x;
    real<lower=0, upper=1> p_ch_6x;
    real<lower=0, upper=1> p_ch_7x;
    real<lower=0, upper=1> p_ch_80_up;
}

Si vous l'écrivez à nouveau, certaines personnes pensent que le taux de létalité en Chine peut être calculé par $ D_ {{\ rm ch}} \ sur N_ {{\ rm ch}} $, donc cela peut être une cible d'estimation. Je pense.

En fait, la létalité en Chine n'a pas encore été déterminée avec ce modèle. Par exemple, si vous lancez un dé deux fois et obtenez un 6 les deux fois, vous ne pouvez pas dire que le dé est un dé avec un 100% 6, donc même si vous recevez le nombre de personnes infectées et le nombre de décès, le taux de létalité est immédiatement déterminé. Vous ne pouvez pas le faire.

L'avantage de la modélisation bayésienne est que vous pouvez facilement gérer des choses aussi incertaines sans être incertain, ce qui est le bon point de Stan.

Maintenant, comme les données, les parameters peuvent convertir des valeurs, alors faisons-le.

transformed parameters {
    //Létalité estimée par âge
    real<lower=0, upper=1> p_0x;
    real<lower=0, upper=1> p_1x;
    real<lower=0, upper=1> p_2x;
    real<lower=0, upper=1> p_3x;
    real<lower=0, upper=1> p_4x;
    real<lower=0, upper=1> p_5x;
    real<lower=0, upper=1> p_6x;
    real<lower=0, upper=1> p_7x;
    real<lower=0, upper=1> p_80_up;

    //Létalité estimée
    real<lower=0, upper=1> p;

    //Létalité chez Diamond Princess
    real<lower=0, upper=1> p_dp;

    p_0x    = c * p_ch_0x;
    p_1x    = c * p_ch_1x;
    p_2x    = c * p_ch_2x;
    p_3x    = c * p_ch_3x;
    p_4x    = c * p_ch_4x;
    p_5x    = c * p_ch_5x;
    p_6x    = c * p_ch_6x;
    p_7x    = c * p_ch_7x;
    p_80_up = c * p_ch_80_up;

    p = (p_0x    * N_ch_0x +
         p_1x    * N_ch_1x +
         p_2x    * N_ch_2x +
         p_3x    * N_ch_3x +
         p_4x    * N_ch_4x +
         p_5x    * N_ch_5x +
         p_6x    * N_ch_6x +
         p_7x    * N_ch_7x +
         p_80_up * N_ch_80_up
         ) / N_ch;

    p_dp = (p_0x    * N_dp_0x +
            p_1x    * N_dp_1x +
            p_2x    * N_dp_2x +
            p_3x    * N_dp_3x +
            p_4x    * N_dp_4x +
            p_5x    * N_dp_5x +
            p_6x    * N_dp_6x +
            p_7x    * N_dp_7x +
            p_80_up * (N_dp_8x + N_dp_9x)
            ) / N_dp;
}

modèle

Il n'y a rien de spécial à dire car c'est comme écrit dans Model Chapter. Cette expressivité est ce qui rend Stan si bon.

model {
    c ~ uniform(0, 1);

    D_ch_0x    ~ binomial(N_ch_0x,    p_ch_0x);
    D_ch_1x    ~ binomial(N_ch_1x,    p_ch_1x);
    D_ch_2x    ~ binomial(N_ch_2x,    p_ch_2x);
    D_ch_3x    ~ binomial(N_ch_3x,    p_ch_3x);
    D_ch_4x    ~ binomial(N_ch_4x,    p_ch_4x);
    D_ch_5x    ~ binomial(N_ch_5x,    p_ch_5x);
    D_ch_6x    ~ binomial(N_ch_6x,    p_ch_6x);
    D_ch_7x    ~ binomial(N_ch_7x,    p_ch_7x);
    D_ch_80_up ~ binomial(N_ch_80_up, p_ch_80_up);

    D_dp ~ binomial(N_dp, p_dp);
}

Courir

Les données

J'ai entré les données suivantes. Celles-ci sont basées sur les données décrites dans le chapitre Hypothèses (#Assumptions).

data = {
    #Nombre de morts sur le Diamond Princess
    'D_dp': 6,

    #Nombre de personnes infectées par groupe d'âge sur le Diamond Princess
    'N_dp_0x':   1,
    'N_dp_1x':   5,
    'N_dp_2x':  28,
    'N_dp_3x':  34,
    'N_dp_4x':  27,
    'N_dp_5x':  59,
    'N_dp_6x': 177,
    'N_dp_7x': 234,
    'N_dp_8x':  52,
    'N_dp_9x':   2,

    #Nombre de décès en Chine
    'D_ch_0x':      0,
    'D_ch_1x':      1,
    'D_ch_2x':      7,
    'D_ch_3x':     18,
    'D_ch_4x':     38,
    'D_ch_5x':    130,
    'D_ch_6x':    309,
    'D_ch_7x':    312,
    'D_ch_80_up': 208,

    #Nombre de personnes infectées par groupe d'âge en Chine
    'N_ch_0x':     416,
    'N_ch_1x':     549,
    'N_ch_2x':    3619,
    'N_ch_3x':    7600,
    'N_ch_4x':    8571,
    'N_ch_5x':   10008,
    'N_ch_6x':    8583,
    'N_ch_7x':    3918,
    'N_ch_80_up': 1408,
}

Environnement d'exécution

Comme pour Last time, en utilisant PyStan Je l'ai couru.

Cliquez ici pour le code source complet

résultat

J'ai défini beaucoup de paramètres, donc quand je l'exécute, j'obtiens d'énormes résultats.

             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
c            0.21  8.7e-4   0.08   0.08   0.15    0.2   0.25   0.38   7792    1.0
p_ch_0x    2.4e-3  3.0e-5 2.3e-3 7.1e-5 7.1e-4 1.7e-3 3.2e-3 8.3e-3   5850    1.0
p_ch_1x    3.6e-3  3.3e-5 2.6e-3 4.3e-4 1.7e-3 3.0e-3 4.9e-3   0.01   6364    1.0
p_ch_2x    2.2e-3  9.7e-6 8.0e-4 9.5e-4 1.6e-3 2.1e-3 2.7e-3 4.0e-3   6894    1.0
p_ch_3x    2.5e-3  6.3e-6 5.7e-4 1.5e-3 2.1e-3 2.5e-3 2.9e-3 3.7e-3   8234    1.0
p_ch_4x    4.6e-3  8.7e-6 7.4e-4 3.2e-3 4.1e-3 4.5e-3 5.0e-3 6.1e-3   7259    1.0
p_ch_5x      0.01  1.3e-5 1.1e-3   0.01   0.01   0.01   0.01   0.02   7562    1.0
p_ch_6x      0.04  2.4e-5 2.0e-3   0.03   0.03   0.04   0.04   0.04   7096    1.0
p_ch_7x      0.08  5.2e-5 4.3e-3   0.07   0.08   0.08   0.08   0.09   6650    1.0
p_ch_80_up   0.15  1.1e-4 9.2e-3   0.13   0.14   0.15   0.15   0.17   7552    1.0
p_0x       4.8e-4  7.1e-6 5.3e-4 1.2e-5 1.3e-4 3.1e-4 6.5e-4 1.9e-3   5530    1.0
p_1x       7.4e-4  8.5e-6 6.3e-4 7.0e-5 3.0e-4 5.8e-4 9.8e-4 2.4e-3   5462    1.0
p_2x       4.5e-4  3.2e-6 2.4e-4 1.3e-4 2.8e-4 4.0e-4 5.7e-4 1.1e-3   5637    1.0
p_3x       5.2e-4  2.7e-6 2.3e-4 1.8e-4 3.4e-4 4.8e-4 6.4e-4 1.1e-3   7336    1.0
p_4x       9.4e-4  4.6e-6 3.9e-4 3.6e-4 6.5e-4 8.8e-4 1.2e-3 1.9e-3   7248    1.0
p_5x       2.7e-3  1.2e-5 1.0e-3 1.1e-3 1.9e-3 2.5e-3 3.3e-3 5.0e-3   7484    1.0
p_6x       7.4e-3  3.1e-5 2.8e-3 3.0e-3 5.4e-3 7.0e-3 9.0e-3   0.01   7796    1.0
p_7x         0.02  6.9e-5 6.1e-3 6.6e-3   0.01   0.02   0.02   0.03   7770    1.0
p_80_up      0.03  1.3e-4   0.01   0.01   0.02   0.03   0.04   0.06   7937    1.0
p          4.7e-3  2.0e-5 1.8e-3 1.9e-3 3.5e-3 4.5e-3 5.8e-3 8.6e-3   7853    1.0
p_dp         0.01  4.7e-5 4.2e-3 4.5e-3 8.3e-3   0.01   0.01   0.02   7865    1.0
lp__        -4215    0.06   2.35  -4220  -4216  -4214  -4213  -4211   1641    1.0

Si vous extrayez la partie importante

             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
c            0.21  8.7e-4   0.08   0.08   0.15    0.2   0.25   0.38   7792    1.0
p          4.7e-3  2.0e-5 1.8e-3 1.9e-3 3.5e-3 4.5e-3 5.8e-3 8.6e-3   7853    1.0

Ce sera.

Autrement dit, ** si l'hypothèse est correcte **, le taux de capture des personnes infectées en Chine est d'environ 20%, et le taux de létalité en considérant cela est d'environ 0,5% et de 0,19% à 0,86% dans l'intervalle de confiance de 95%. A été estimé.

Comme mentionné précédemment, en réalité ** cette hypothèse n'est pas exacte **, donc aucun intervalle de confiance à 95% ne peut être fiable, mais le taux de létalité sur le Diamond Princess, où la majorité des passagers avaient plus de 60 ans. Cependant, comme la mesure réelle est d'environ 1%, je m'attends à ce que cette estimation ne soit pas largement manquée.

Résumé

Cette fois, contrairement à précédent, nous avons pu estimer la létalité du nouveau virus corona en utilisant Stan.

Je pense que j'ai pu exprimer la bonté de Stan plus qu'avant. (Je ne pourrais pas trop le faire)

En conséquence, nous avons obtenu une létalité estimée de 0,19% à 0,86%, ce qui correspond à la situation infectieuse en Chine, bien que nous ayons dû cracher les sourcils car nous avons fait une hypothèse forte.

De côté

Dans le post précédent, [j'ai fait une estimation en convertissant la distribution binomiale en une distribution bêta](https://qiita.com/akeyhero/items/894dd3b5c206325191ce#%E3%83%99%E3%83%BC% E3% 82% BF% E5% 88% 86% E5% B8% 83% E3% 81% A7% E7% 84% A1% E7% 90% 86% E3% 81% 8F% E3% 82% 8A% E8% A7% A3% E3% 81% 84% E3% 81% A6% E3% 81% BF% E3% 82% 8B), mais cette fois encore, le nombre de décès parmi les personnes infectées de la tranche d'âge $ a $ en Chine La modélisation

D_{{\rm ch}}^a \sim {\rm binomial}(N_{{\rm ch}}^a, p_{{\rm ch}}^a)

ne pas,

p_{{\rm ch}}^a \sim {\rm beta}(D_{{\rm ch}}^a + 1, N_{{\rm ch}}^a - D_{{\rm ch}}^a + 1)

Vous pouvez en déduire de la même manière en le remplaçant par. Si vous êtes intéressé, essayez-le.

Recommended Posts

[Suite] La nouvelle Corona est-elle vraiment une menace? Létalité calculée avec Stan
La nouvelle Corona est-elle vraiment une menace? Validé avec Stan (était)
[New Corona] Le prochain pic est-il en décembre? J'ai essayé l'analyse des tendances avec Python!