[PYTHON] [Fortsetzung] Ist die neue Corona wirklich eine Bedrohung? Berechnete Letalität mit Stan

Der vorherige Beitrag "Ist die neue Corona wirklich eine Bedrohung? Ich habe sie mit Stan überprüft (war)" war in den folgenden zwei Punkten unverdaulich:

  1. Schließlich weiß ich nicht, wie hoch die Letalitätsrate von COVID-19 ist.
  2. Schließlich arbeitet Stan nicht gut

Aus Rache werde ich diesmal versuchen, mit Stan die "wahre Letalitätsrate" zu ermitteln.

Ich habe nicht wiederholt erklärt, was ich in Vorheriger Artikel erklärt habe. Bitte beachten Sie auch dies: bow:

Haftungsausschluss

Ich bin kein Experte für Bayes'sche Modellierung oder Stan, daher kann der Inhalt dieses Artikels Fehler und bessere Techniken enthalten.

Auch wenn Sie etwas falsch finden, würde ich es begrüßen, wenn Sie mich wissen lassen könnten. : Bogen:

Annahme

In meinem letzten Beitrag habe ich argumentiert, dass die Infektionsuntersuchung [^ 1] zum Thema Diamond Princess auf der Annahme beruhte, dass sie die genaueste der Welt sei. Wir werden diese Politik auch diesmal befolgen.

Die Anzahl der Fälle ist jedoch bei der Diamond Princess allein gering, und es ist nicht möglich, die Letalitätsrate nach Alter zu ermitteln. Daher werden wir die Daten [^ 2] verwenden, die sich auf Chinas Kontrolle der Krankheitskontrolle konzentrieren.

Da die Populationen dieser beiden Daten völlig unterschiedlich sind, setze ich Annahmen und integriere sie irgendwie. Hier nehmen wir an, dass die Anzahl der Infektionen in Chinas Daten von allen wirklich infizierten Personen mit einer bestimmten Wahrscheinlichkeit von $ c $ unabhängig von Altersgruppe oder Mortalität erhalten wurde.

Diese Annahme ist eindeutig ungenau, da die Altersgruppen unterschiedliche symptomatische Tendenzen aufweisen und daher der Anteil der gefangenen infizierten Personen variieren kann. Daher sind die Schlussfolgerungen in diesem Artikel aufgrund dieser Annahme ungenau. Nach Angaben des Nationalen Instituts für Infektionskrankheiten [^ 1] ist der Anteil der Patienten mit Symptomen (symptomatisch bestätigte Fälle) an denen ohne Symptome (asymptomatisch bestätigte Fälle) in der jüngeren Generation jedoch eher geringer. Es gibt auch Berichte [^ 3], dass der Schweregrad junger Menschen unerwartet hoch ist. Aus diesem Grund denke ich, dass es eine absurdere Annahme als Intuition ist.

Die Letalität in der Diamond Princess kombiniert die beiden Daten unter der Annahme, dass sie der Letalität der Altersgruppe folgt, die unter Verwendung der hier angenommenen Erfassungsrate von $ c $ berechnet wurde.

Andererseits wird davon ausgegangen, dass die Gefangennahme der Toten in der chinesischen Statistik ausnahmslos erfolgt.

Als nächstes gab es sechs Todesfälle bei der Diamond Princess. Ein siebter Tod wurde am 7. März bestätigt, aber es ist möglich, dass er zu diesem Zeitpunkt nicht infiziert war, da seit den Daten des Nationalen Instituts für Infektionskrankheiten mehr als zwei Wochen vergangen sind [^ 1]. Es basiert auf der Tatsache, dass es hoch ist, es wurde nicht bekannt gegeben, ob es sich um einen COVID-19-bedingten Tod handelt, und dass die Anzahl der infizierten Personen höher war als zum Zeitpunkt der Daten vom 7. März.

Schließlich gibt es Unterschiede nach Rasse, Unterschiede in der Verteilung infizierter Menschen (tatsächlich sollte die Verteilung infizierter Menschen in China und die Verteilung infizierter Menschen in der Welt unterschiedlich sein) und Unterschiede im medizinischen Umfeld (insbesondere wenn die Anzahl der Patienten schnell zunimmt, tritt ein medizinischer Kollaps auf. Ich habe beschlossen, nicht in Betracht zu ziehen (was die Letalitätsrate erhöhen sollte). Der Unterschied in der Anzahl der Infizierten zwischen China und der Diamantprinzessin wird jedoch ausnahmsweise berücksichtigt.

Die Annahmen können wie folgt zusammengefasst werden.

Modell-

Unten ist $ N $ die Anzahl der Infizierten, $ D $ die Anzahl der Todesfälle und $ p $ die Letalitätsrate. Sofern nicht anders angegeben, gibt der hochgestellte Index die Altersgruppe und der tiefgestellte Index den Ort an ($ {\ rm dp} $ = Diamond Princess, $ {\ rm ch} $ = China). Das Fehlen von hochgestellten Zeichen repräsentiert alle Altersgruppen.

Zum Beispiel gibt $ N_ {{\ rm dp}} ^ {60-69} $ "die Anzahl der infizierten Personen im Alter von 60-69 an, wie auf der Diamond Princess angegeben". (Eigentlich werde ich von nun an keine bestimmte Altersgruppe mehr erwähnen, sondern $ N_ {{\ rm dp}} ^ {a} $ schreiben und "Infektion der Altersgruppe $ a $ auf der Diamond Princess bestätigt" schreiben. Es wird verwendet, um die Anzahl der Personen darzustellen.)

Nehmen wir an, dass die geschätzte Letalität $ p $ beträgt, wenn die Annahmen korrekt sind. Außerdem wird diese geschätzte Letalität $ p $ nach Altersgruppen geteilt und als $ p ^ {60-69} $ (ohne Indizes) ausgedrückt.

Die Anzahl der Todesfälle bei infizierten Menschen in der Altersgruppe $ a $ in China folgt einer binären Verteilung.

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

Kann modelliert werden als.

In ähnlicher Weise ist auch die Anzahl der Todesfälle bei infizierten Menschen in der Diamond Princess gleich.

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

Es wird sein. Bitte beachten Sie, dass wir nicht nach Altersgruppen modellieren, da wir nicht die Altersverteilung der Todesfälle bei Diamond Princess haben.

Verbinden Sie die beiden mit der Erfassungsrate $ c $, die im Kapitel Annahmen (#Annahmen) angegeben ist.

Da keine Informationen zur Erfassungsrate vorliegen, gehen wir von einer gleichmäßigen Verteilung aus.

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

Wir haben angenommen, dass $ c $ unabhängig von Altersgruppe und Letalität konstant ist, also unabhängig von $ p_ {{\ rm ch}} ^ a $. Die geschätzte Letalität der Altersgruppe $ a $ ist also

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

In jedem Alter

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}

Es wird sein.

Wir haben auch angenommen, dass die Letalität der Diamantprinzessin der geschätzten Letalität folgt ($ p_ {{\ rm dp}} ^ a = p ^ a $), also ist $ p_ {{\ rm dp}} $

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

Es wird sein.

Implementierung

Jetzt, da wir alle Schauspieler haben, werden wir es in Stan implementieren.

Eingabedaten

Der Variablenname lautet "6x" (dh $ N {{\ rm ch}} ^ {60-69} $ = "N_ch_6x"), wenn die Altersgruppe 60-69 Jahre alt ist, und "80_up" (dh $ N), wenn die Altersgruppe 80 Jahre oder älter ist. Es ist fast dasselbe wie das Symbol [Modellkapitel](# Modell), außer dass es als {{\ rm ch}} ^ {80-} $ = N_ch_80_up) ausgedrückt wird.

Bitte beachten Sie, dass die Diamond Princess-Daten [^ 1] Daten für die 90er Jahre enthalten, die chinesischen Daten [^ 2] jedoch für Personen über 80 Jahre.

data {
    //Anzahl der Todesfälle bei der Diamond Princess
    int D_dp;

    //Anzahl der Infizierten nach Altersgruppen auf der 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;

    //Anzahl der Todesfälle nach Altersgruppen in China
    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;

    //Anzahl der Infizierten nach Altersgruppen in China
    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;
}

Da Stan der Einfachheit halber verarbeitete Daten definieren kann, berechnen wir die Werte für alle Altersgruppen.

transformed data {
    //Anzahl der infizierten Personen auf der Diamond Princess
    int N_dp;

    //Anzahl der Infizierten in China
    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;
}

Geschätzter Zielparameter

Geben Sie die zu schätzenden Parameter an.

parameters {
    //Erfassungsrate c
    real<lower=0, upper=1> c;

    //Tödlichkeit in China
    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;
}

Wenn Sie es noch einmal schreiben, denken einige Leute, dass die Letalitätsrate in China durch $ D_ {{\ rm ch}} \ über N_ {{\ rm ch}} $ berechnet werden kann, so dass es sich möglicherweise um ein Schätzziel handelt. Ich denke.

Tatsächlich wurde die Letalität in China mit diesem Modell noch nicht bestimmt. Wenn Sie beispielsweise zweimal würfeln und beide Male eine 6 erhalten, können Sie nicht sagen, dass es sich bei den Würfeln um Würfel mit einer 100% 6 handelt. Selbst wenn Sie die Anzahl der infizierten Personen und die Anzahl der Todesfälle erhalten, wird die Letalitätsrate sofort bestimmt. Du kannst es nicht tun.

Das Gute an der Bayes'schen Modellierung ist, dass Sie mit solchen unsicheren Dingen leicht umgehen können, ohne unsicher zu sein, was der gute Punkt von Stan ist.

Jetzt können Parameter wie die Daten Werte konvertieren, also lass es uns tun.

transformed parameters {
    //Geschätzte Letalität nach Alter
    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;

    //Geschätzte Letalität
    real<lower=0, upper=1> p;

    //Tödlichkeit in der Diamantprinzessin
    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;
}

Modell-

Es gibt nichts Besonderes zu sagen, da es wie in Modellkapitel geschrieben ist. Diese Ausdruckskraft macht Stan so gut.

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);
}

Lauf

Daten

Ich habe folgende Daten eingegeben. Diese basieren auf den im Kapitel Annahmen (#Annahmen) beschriebenen Daten.

data = {
    #Anzahl der Todesfälle bei der Diamond Princess
    'D_dp': 6,

    #Anzahl der Infizierten nach Altersgruppen auf der 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,

    #Anzahl der Todesfälle in China
    '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,

    #Anzahl der Infizierten nach Altersgruppen in China
    '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,
}

Ausführungsumgebung

Wie Letztes Mal mit PyStan Ich habe es ausgeführt.

Klicken Sie hier für den gesamten Quellcode

Ergebnis

Ich habe viele "Parameter" eingestellt, daher erhalte ich beim Ausführen einige große Ergebnisse.

             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

Wenn Sie den wichtigen Teil extrahieren

             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

Es wird sein.

Das heißt, ** wenn die Annahme richtig ist **, beträgt die Erfassungsrate infizierter Personen in China etwa 20% und die Letalitätsrate unter Berücksichtigung von etwa 0,5%, 0,19% bis 0,86% im 95% -Konfidenzintervall. Wurde geschätzt.

Wie bereits erwähnt, ist diese Annahme in der Realität ** nicht korrekt **, sodass kein 95% -Konfidenzintervall als vertrauenswürdig eingestuft werden kann, sondern die Letalitätsrate der Diamond Princess, bei der die Mehrheit der Passagiere über 60 Jahre alt war. Da die tatsächliche Messung jedoch etwa 1% beträgt, gehe ich davon aus, dass diese Schätzung nicht weitgehend verfehlt wird.

Zusammenfassung

Dieses Mal konnten wir im Gegensatz zu vorherigen die Letalität des neuen Koronavirus mithilfe von Stan abschätzen.

Ich glaube, ich konnte die Güte von Stan mehr als zuvor ausdrücken. (Ich konnte es nicht zu viel tun)

Als Ergebnis erhielten wir eine geschätzte Letalität von 0,19% bis 0,86%, was der Infektionssituation in China entspricht, obwohl wir auf unsere Augenbrauen spucken mussten, weil wir eine starke Annahme machten.

Beiseite

Im vorherigen Beitrag [habe ich eine Schätzung vorgenommen, indem ich die Binomialverteilung in eine Beta-Verteilung konvertiert habe](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), aber auch diesmal die Zahl der Todesfälle bei infizierten Menschen der Altersgruppe $ a $ in China Modellieren

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

nicht,

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

Sie können auf die gleiche Weise schließen, indem Sie es durch ersetzen. Wenn Sie interessiert sind, probieren Sie es bitte aus.

Recommended Posts

[Fortsetzung] Ist die neue Corona wirklich eine Bedrohung? Berechnete Letalität mit Stan
Ist die neue Corona wirklich eine Bedrohung? Validiert mit Stan (war)
[New Corona] Ist der nächste Höhepunkt im Dezember? Ich habe die Trendanalyse mit Python versucht!