実効再生産数Rt計算 エクセルで簡単EpiEstimは何を計算しているか
御無沙汰しています。
すっかり忘れてしまいそうなので以前理解した内容を書いておきます。
エクセルでも簡単に実効再生産数Rtを計算できるEpiRstimについて
EpiEstimの動作は、Step1、Step2からなる
Step1 SIの分布(ws)を計算。
ここで、SI(Sireal Interval)は、一人目の発症からその人がうつした人が
発症するまでの期間。
Step2 Step1の結果を用いてRtを計算する。
過去の発症者の数(I(t-τ), τ=1,2,..s)にwsを畳み込み積分し(sはwsの最大長)、
I(t)をその値で割ったものがRt
とってもシンプル♪
ちなみに、ここまで紹介してきた以下の設定では、Step1は決め打ちにしています。
なので、実際に計算しているのはStep2のみです。
キーワード コロナ 新型コロナ COVID-19
【新型コロナ】実効再生産数Rtをエクセルで計算したら1ヶ月後が見えた?(WHOお墨付きEpiEstimで)(3)
一昨日、昨日に続き、WHOお墨付きEpiEstimのexl版で実効再生産数を計算しました。
実効再生産数Rtは1人目から2人目にうつした倍率である事を利用すると、1ヶ月後が見えてくるかも知れない。
さっそく、エクセルで計算だ!
全国のRtはEpiEstimで0.9だったので、もしその値がずっと変わらないと想定して計算してみると、1ヶ月後のpcrポジティブは激減する。入院治療を要する人、重篤患者数も同様に減少して、1ヶ月後をこのようにざっくり試算できる。
2020/9/27までの直近2週間ほどのデータを用いて、2020/9/24のRtを推定した。
この時のグラフはこのようになる。あくまで、1ヶ月間Rtが全く変わらない想定のグラフです。。ありえないけれども。
試しに、Rtが1以上になる沖縄のデータで同様に計算すると、1ヶ月後に感染者が倍になる事になる。
日本全国都道府県のRtを計算した結果が以下のようになる。
水色がRt<1で縮小傾向、黄色は1以上2.5未満で黄色信号、ピンクは2.5以上でR0より大きい値の県になる。といっても、元々感染者数が少ない県では少し感染者が発生しただけで大きな値になってしまうため、信頼度は疑問。でも、中国地方が目立つね。何か特別なことがあったんでしょうか?
1ヶ月後の入院治療を要する患者数をざっくり推定するとこのようになった。ピンクは病床数を上回る。黄色はかなりの割合を占めなくてはならないことになる(他の病菌お患者さんもいるので現実的には難しいと予想される。)。島根、広島、鹿児島は今の状態で感染者数が増加すると患者を収容できなくなるという事になる。
キーワード PC xls COVID-19
【新型コロナ】実効再生産数Rtをエクセルで計算してみた(WHOお墨付き!?)(2)
昨日の続きです(ただし、使用するデータを2020/02/27-09/26に更新しました)。
WHOお墨付きのEpiEstimのexl版で出力したデータをエクセル上で加工します。
(1) 新しいエクセルシートにEpiEstimのOutput2 R estimatesタブからRのmeanの計算値をコピーし、Time periodの値に気をつけて、対応する年月日を貼り付ける。
(2) こんなグラフが簡単に描ける。
(2) こんなグラフが簡単に描ける。
(3) 左右の0.05quantileと0.95quantaileも使ってあげると、もっとカッコイイ(?)
【参考文献】
WHO(2020) Public health criteria to adjust public health and social measures in the context of COVID-19
Anne Cori(2013) A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics
【新型コロナ】実効再生産数Rtをエクセルで計算してみた(WHOお墨付き!?)
前回の書き込みから御無沙汰しております。
猛反省して新型コロナ(COVID-19)の論文を調べ、紆余曲折してはやウンヶ月。PythonとRのインストールにからかわれておりました。^^;
世界保健機構(WHO)報告(3月12日)から、再生産数の定義は以下
Rt (the effective number of secondary cases per infectious case in a population)
1人の感染者が2人目にうつす数。
SIは、1人目が発症してから、うつした人が発症するまでの期間。
そこにはRプログラムが添付されていて(https://CRAN.R-project.org/package=EpiEstim)
その作者Anne Coriさんの論文A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemicsでは、EXCELもリリースされていました。
http://tools.epidemiology.net/EpiEstim.xls
つまり、WHOお墨付き!?
この論文、2013年です。専門家の方ならとうにご存じですね。orz
と思いつつご紹介
使い方もいたってシンプル。
EpiEstim.xlsをダウンロード
↓
DataタブのIncidenceに読み込みたいデータを貼り付ける
↓
Estimate R!
↓
結果がFigureに表示される
計算内容は、この論文もさることながら本人が監修しつつ他の方が書いている論文が簡潔で分かり良いです。こちらも機会があればご紹介します。
さっそく、日本の新型コロナのRtを計算してみましょう。
(1) http://tools.epidemiology.net/EpiEstim.xlsからダウンロード
(2) コンテンツの有効化
(3) Dataタブを開くと
(4) オリジナルは他の感染症のデータになっているので、COVID19Serial Intervalの値に書き換える。ちょっと古いかも知れないがWHOで5-6日と言っていたので、mean5.5日、SD1くらいにしてみよう。
(5) 日本のデータをIncidenceに貼り付ける。このとき、データ長を変えた場合には、Timeを変更してMax Timeもそれに合わせて変更する。
厚生労働省のデータ(https://www.mhlw.go.jp/stf/covid-19/open-data.html)を使おうとすると、”陽性者数pcr_positive_daily.csv”と”入院治療を要する者の数”cases_total.csvのどちらを使うべきか迷う事になる。Rtは1st case→2nd caseの倍率なので、必要なのはcases_dailyであるはずだ。しかし、ここでcases_total.csvからdailyを作ろうとすると、マイナスの値が続々発生して、、どうも定義が違うらしいorzということで断念して、”陽性者数pcr_positive_daily.csv”を使ってみた。
(6) Estimate R!
(7) Figuresタブが開いて、Incidence、Serial interval distribution、R averaged over time periods (posterior median and 95%CrI)の3つのグラフが表示される。
Rの計算の最初の方は少々おかしいが、まぁご愛敬ということで。
日本の2020/3/1~2020/9/25の結果
論文についているPythonやらRやらのプログラムを動かそうとして、インストールだのバージョン問題だのに振り回されて毎度1ヶ月以上要していた事を考えると、Excelは神
ただし、中はいじるなとのことで、自分の研究の参考にするのはちょっと難しそうかな。
【参考文献】
WHO(2020) Public health criteria to adjust public health and social measures in the context of COVID-19
Anne Cori(2013) A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics
キーワード 実効再生産数
【新型コロナ】何も対策しない場合の死者数がなんと!?
前回の反省から、ちゃんと修正しようと思っていろいろ調べておりましたが
思いの外はまってしまい、このままでは飽きてしまって風化させる恐れがあるので
これまで分かったところだけでも書こうと思います。
「何もしなければ41万人死亡」というのは4月に北大の西浦氏が発表された内容で、ショッキングだったので噂になりましたが、少し前の事なのでもしかするとご本人も何かアップデートされているかも知れません。
私は疫学については専門家ではないので、、とお断りした上で、
WHOの3/4付けのCOVID-19とインフルエンザの比較記事のデータを元に、
基本的な疫学モデルであるSIRモデル(1930年頃に発表されその後も脈々とマイナーチェンジされている)に入れて、エクセルで計算しました。
SIRモデルは、S(感染していない)、I(感染者)、R(回復、もしくは死亡)という群を想定しています。
絵にするとこんな感じです。一度感染して回復した人は、免疫がついて再度感染しないという事がミソになっています。
ここで、WHOの基本再生産数R0は、感染する強さを表し、発症間隔(SI)は、感染を次に引き継ぐ時間を表します。
そして、パラメータを入れてSIRモデルでS、I、Rをプロットするとこんな感じ。
ここで、Rというのは一度感染した人の累積になり、最終的に94%が感染する結果になっています。
一昨日の厚生労働省のデータから、感染者に対する死亡率は約5.7%と計算できる(この値がそこそこ変動する点に注意が必要)ので、
人口1億2千万の94%の5.7%といえば、死者数670万人という事になる、、えっっ何か間違ってない?これ!??
まぁそんなとこで。。^^;
ピークもCOVID-19はインフルエンザの2.5倍なので医療崩壊必至な様子です。
ところで、実際のインフルエンザの感染がこの結果と同じかというとそうではありません。インフルエンザでは、ワクチンが投与されているからです。しかし、ワクチンは、インフルエンザを流行させないように投与されているはずなのにも関わらず毎年1000万人も感染しているっておかしくないですか?
そう、単純なSIRモデルでは、インフルエンザの感染者数は計算できません。ここまでのSIRモデルは平均値で話をしていますが、値にはばらつきがあるので、幅を持っている。だから、なんちゃら分布にフィッティングという話が始まるんでしょう。
、、にしても、1000万人って多くない?
そう、この1000万人が、しかも、毎年推定1000万人とか言っており、生データが拾えなかった、、orz
↑それではまっていたんですね。(TT)
ほとんどの人が保険診療なのだから生データがありそうなものだけれど、保険データは見つからず、都道府県別データは、定点観測とかいう昭和初期から続いているようなデータで、、
インフルエンザのデータは都道府県別だけど、インフルエンザの種類別データは都道府県別ではないとか(そもそも、SIRモデルは単一の病原体でなければ想定がおかしいでしょう)、、国立感染症研究所はデータの提供方法がスマートでありがたかったのですが、、内容が結果的に使えなくて残念だった。それとも私の探し方が悪かったのか。
もしかして、これもあれかしら?
日教組と文部科学省が仲が悪いように、厚生労働省と保険機構も仲が悪いのかしら?
と、余計なことを考える今日この頃。
新型コロナ対策を考えるにあたって、インフルエンザの都道府県別生データがあれば参考になるかと思います。定点観測データは、専門家の方には有用なデータなのかと思いますが、独特な値なので、おそらく海外で通用しないのではないかと思いました。
最後に、謝辞
多くのCOVID-19関連の論文や周辺論文が広く無料でアクセスできる状態で提供されており、本当に、人類一丸となってコロナと戦おうとしているのだと、感動しました。
関係者の皆様方、本当にありがとうございます。
【新型コロナ】レムデシビルの医療崩壊期間の短縮効果をエクセルで計算してみた。
*お断り:申し訳ありません。この計算に用いたデータは古いデータ(4月22日時点)のため、汎用性はありません。
先日認可が下りた新型コロナの薬レムデシビルと、今後1年以内の認可を目指すイベルメクチンについて、事前投与した場合、推定死者数に対する効果を計算した(レムデシブルとイベルメクチンの効果をエクセルで計算してみた(重症化前投与))。
ここでは、より現実的な投与方法を想定して、重症化後の薬の投与を想定した計算をして比較を行う。
<従来法>
・従来データを用いる
・集中施設占有期間は2週間
<レムデシビル条件>
・重症化前投与条件:重症化数を70%に改善
・重症化後投与条件:死者数を70%に改善
・集中施設占有期間を70%に改善
* 重症化数:重症患者数+死者数
<イベルメクチン条件>
・重症化前投与条件:重症化数を70%に改善
・重症化後投与条件:死者数を70%に改善
・集中施設占有期間は2週間据え置き
重症化後投与の結果は以下になる。
結果として、おおまかな推定死者数はほぼ変わらず、従来は41万人、レムデシビル使用時に約29万人、イベルメクチン使用時に約12万人である。
新型コロナ感染者に対する医療崩壊率も、大きな違いは無いが、ここでは、医療崩壊による他の疾患患者への影響が考えられていないためである。
この計算では、インフルエンザの感染者の季節推移を参考にして、ピーク期には1週間に年間感染者数の1/10が発生し、中シーズンには年間感染者数の1/100が発生すると仮定している。医療崩壊の条件は、日医総研のデータより各都道府県の集中治療室の和を重症患者数+死者数の和が上回る場合とした。
この仮定において、レムデシビルの医療崩壊期間短縮効果は以下のようになる。
この重症化後投与(レムデシビル(後))と重症化前投与(レムデシビル(前))の差は、山梨県や岩手県では中シーズンに医療崩壊の瀬戸際に立つため、レムデシビルを陽性患者や重症化しやすい患者に事前投与する事によって、医療崩壊を免れる可能性があるという事を示す。その後のピーク期には医療崩壊が起こるため、医療崩壊期間の短縮効果として現れる。しかしそれは、薬の副作用との兼ね合いになるのだろう。
ピーク期には全都道府県で医療崩壊が起こる一方、中シーズンにおいて医療崩壊する都道府県は以下のように予測される。
*新潟~広島の16県:広島、三重、滋賀、大分、栃木、和歌山、山形、愛知、静岡、島根、岐阜、徳島、青森、岩手、山梨、新潟
*注意:疫学の非専門家による大胆な仮定に基づく計算です。お気づきの点がございましたらご指摘頂けますとありがたいです。
キーワード
COVID-19 コロナウイルス エクセル excel 統計
【新型コロナ】レムデシビルでも死者30万人!エクセルで計算してみた
*お断り:申し訳ありません。この計算に用いたデータは古いデータ(4月22日時点)のため、汎用性はありません。
AFPによれば、レムデシビルの効果は治療期間を31%減らし、致死率も31%減らすとのことだ。この値はアメリカでの治験によるので日本でも同じかどうかは分からないが、仮に同じだとして、従来の死亡数推定数から効果を見積もってみよう。
また、先の計算では集中治療を1週間と見積もっていたが、報道発表では2週間であったのでその修正も行った。結果、レムデシビル使用時は約10日になる。
一方、報道発表によれば、イベルメクチンは致死率を8.5%から1.4%に改善したとの事だ。これも同様に見積もってみよう。
レムデシビルは既に5月7日に認可が下り、アビガンは5月中に認可が下りる予定で、イベルメクチンは1年以内の認可を目指しているそうだ。
早速、エクセルで計算してみよう。
<前提条件の変更>
・集中施設占有期間は2週間 (←先の計算では1週間)
<レムデシビル条件>
・重症化数を70%に改善 *重症化する前の投与を前提とした。
・集中施設占有期間を70%に改善
* 重症化数:重症患者数+死者数
<イベルメクチン条件>
・重症化数を30%に改善 *重症化する前の投与を前提とした。
・集中施設占有期間は2週間据え置き
重症化する前の投与を想定したため、新型コロナ感染者に対する医療崩壊率は、従来、レムデシビル、イベルメクチンでそれぞれ92.5%、90.2%、90%となった。
結果は、以下のようになる。
レムデシビルを用いても推定死者数は約30万人となり、薬が前評判通りに効いたとしても、推定死者数はまだまだ多い印象だ。
キーワード
COVID-19 コロナウイルス エクセル excel 統計