基于Dirichlet 模型的國內(nèi)COVID-19 疫情數(shù)據(jù)的Bayesian 估計

楊云源， 楊新平， 張 潔

（1.楚雄師范學院地理科學與旅游管理學院，云南楚雄675000； 2.楚雄師范學院數(shù)學與統(tǒng)計學院，云南楚雄675000；3.楚雄醫(yī)藥高等專科學校，云南楚雄675000）

截至2020 -04 -10，全國 COVID -19 確診83 369例，海外確診1 537 039 例.疫情在美國、意大利、西班牙、德國等地急劇爆發(fā).通過合理的數(shù)學模型，分別計算我國COVID -19 4 個概率指標（治愈概率、重癥概率、病死概率、輕癥概率）和4 個總量指標（累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)）的點估計及區(qū)間估計分別進行計算，以期為國內(nèi)和海外COVID -19 疫情監(jiān)測、指標估計和疫情阻斷提供參考.

與治愈率、重癥率、病死率、輕癥率相關的研究，主要集中在疫情數(shù)據(jù)描述統(tǒng)計和疫情相關指標預測（估算）上.

1）疫情數(shù)據(jù)描述統(tǒng)計.如中華預防醫(yī)學會COVID-19 防控專家組［1］基于截至 2020 -02 -10 的42 780 例確診病例數(shù)據(jù)統(tǒng)計，指出COVID -19 病死率為2.38%，合并基礎疾病的老年男性病死率較高.中國疾病預防控制中心COVID -19 應急響應機制流行病學組［2］對截至2020 -02 -11 的72 314例4 類病例的年齡、性別、旅居史等進行回顧性分析，揭示疫情自首次報告病例日后30 d 蔓延至31個?。▍^(qū)/市），粗病死率為2.3%.

2）疫情相關指標預測（估算），主要圍繞再生數(shù)、感染規(guī)模、倍增時間等估算展開.如Wu 等［3］使用馬爾科夫-蒙特卡洛法（MCMC）估算武漢COVID-19 再生數(shù)為2.68，流行倍增時間為6.4 d.陳端兵等［4］使用蒙特卡洛法（MC）依據(jù)少量患者的癥狀出現(xiàn)時間對有效再生數(shù)進行估計，表明我國在防控措施出臺兩周內(nèi)，COVID-19 傳播得到了有效遏制.其他方法如：使用 SEIR 或修正 SEIR 模型［5-7］、自回歸移動平均模型（ARIMA）［8］、隨機傳輸模 型［9］、最 大 似 然 法［10］、自 助 抽 樣 法 （Bootstrap）［11］進行 COVID-19 疫情估算.以上文獻中，文獻［1 -2］基于病例數(shù)據(jù)的統(tǒng)計分析較全面，描述統(tǒng)計結(jié)果具有一定的可信度；模型使用中，主要針對再生數(shù)、感染規(guī)模等進行估算；參數(shù)確定直接關乎模型精度，如文獻［5 -6］對相同研究區(qū)（湖北省、武漢市）使用相同方法，不同參數(shù)得到不同估算結(jié)果；文獻［7］的疫情估算使用數(shù)據(jù)、方法較全面，首次使用了AI技術，預測國內(nèi)疫情于4 月底結(jié)束.對于COVID-19 的治愈率、重癥率、病死率、輕癥率等重要指標，主要通過描述統(tǒng)計來量化，通過模型估計的較少.

Bayesian方法在醫(yī)學運用方面，高文龍等［12］、張繼巍等［13］分別簡要介紹了 GeoBugs、OpenBugs在醫(yī)學中的應用方法.貝葉斯方法可以量化給定數(shù)據(jù)條件下的不確定性參數(shù)值，借助先驗信息和數(shù)據(jù)信息可提高數(shù)據(jù)分布參數(shù)的量化精度［14］.如 Yarnell等［15］將貝葉斯分析用于危重病醫(yī)學中.Thall［16］將Bayesian方法用于癌癥術后發(fā)病率估算與藥物劑量確定.Li等［17］使用Bayesian方法來確定不同非藥物療法對失眠癥的療效和安全性的差異.

對疫情指標采用概率意義下的可信區(qū)間估計不僅能得到誤差限而且精度高，有一定的實際意義.本研究根據(jù)國家衛(wèi)建委2020 -01 -20 至03 -28 公報數(shù)據(jù)，建立 Bayesian 框架下的 Dirichlet 模型，通過MCMC 方法進行求解，得到逐日4 個概率指標（治愈概率、重癥概率、病死概率、輕癥概率）的點估計與區(qū)間估計，進一步計算逐日4 個總量指標（累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)）的點估計與區(qū)間估計，并進行誤差分析，給出了估計誤差限.以期研究方法與結(jié)果能為國內(nèi)和海外COVID -19 疫情監(jiān)測、指標估計和疫情阻斷提供參考.

1 數(shù)據(jù)和方法

收集國家衛(wèi)建委2020 -01 -20 至03 -28 公布的累計確診病例數(shù)n、累計治愈病例數(shù)c、現(xiàn)有重癥病例數(shù)e、累計病死病例數(shù)d數(shù)據(jù).累計輕癥病例為 m，則

將 n、c、e、d、m 視為總體，第 t天的觀測值分別記nt、ct、et、dt、mt.2020 -01 -20 至 01 -22 共 3 d 公報數(shù)據(jù)中的e1＝0、d1＝0，并且 3 d 的新增治愈數(shù)為0.所以本研究使用的樣本數(shù)據(jù)為2020 -01 -23 至03 -28 共66 d的日報數(shù)據(jù).

COVID-19 確診病例逐日病情診斷、治療康復出院前核算檢測都需要一定的時間，加之前期檢測試劑、醫(yī)療資源短缺，上報全國的疫情的各項日報指標值工作量較大.數(shù)據(jù)的誤差可能是來自隨機因素，也可能是系統(tǒng)誤差，該誤差不能消除，只能設法減?。?8］，采用可信區(qū)間對COVID -19 疫情指標進行估計就變得更有實際意義.用帶有一定誤差的數(shù)據(jù)對治愈率、重癥率、病死率、輕癥率等相對指標進行計算，誤差會累積，而Bayesian框架下的Dirichlet分布是解決此類問題的有力的工具.將國家衛(wèi)健委公布的累記確診病例數(shù)視為當日的獨立試驗次數(shù)（當日每個確診病例病情變化視為一次實驗過程），結(jié)合試驗的結(jié)果（日報值），以多項分布為基礎，建立Bayesian框架下的Dirichlet模型，借助數(shù)據(jù)信息及先驗信息，通過Gibbs抽樣方法得到總體參數(shù)的后驗樣本，采用MCMC方法得到治愈概率pt1、重癥概率pt2、病死概率pt3及輕癥概率pt4為參數(shù)的Bayesian 估計，進一步計算總量指標（治愈病例數(shù)、重癥病例數(shù)、病死病例數(shù)、輕癥病例數(shù)）的Bayesian點估計和區(qū)間估計.使用軟件為Eviews、R語言和OpenBUGS.

2 Bayesian框架下的Dirichlet模型的建立

2.1 模型的建立將第t 天累計確診病例數(shù)nt視為nt次獨立試驗，當日每個確診病例的病情變化視為一次實驗，每次試驗有4 個互不相容的結(jié)果：At1＝“治愈”，At2＝ “重癥”，At3＝ “病死”及At4＝“輕癥”，且只有一種結(jié)果發(fā)生（無癥狀感染者日報起始于 2020 -04 -01），P（Ati）＝ pti，i ＝ 1，2，3，4.記待估參數(shù)為

日報數(shù)據(jù)記為

由多項分布的定義有：

datat～ Multinomial（θt，nt），相應的似然函數(shù)為模型（1）

模型（1）中，

因為參數(shù)θt沒有任何先驗信息，按照 Nguyen的理論［19］，選擇 Dirichlet 分布為其先驗分布，第 t天先驗分布模型為模型（2）

記 α ＝ （α1，α2，α3，α4）T，從而得到θt的后驗精確分布

θt～ Dirichlet（α1+ct，α2+et，α3+dt，α4+mt），

相應的后驗密度函數(shù)為模型（3）

模型（3）中，參數(shù)θt沒有任何先驗信息，根據(jù)Lunn等的理論［20］，對于先驗分布 p（θt），選擇 α ＝（1，1，1，1）T，該分布等價于四維的均勻分布作為θt先驗分布.分別求出各個參數(shù)的滿條件分布p（θti｜θt（-i），datat，α），其中，θt（-i）表示θt中 4 個分量去掉第i個分量后剩下的3 個分量構成的向量.比如

參數(shù)θt1相應的滿條件分布為模型（4）

模型（4）是一個標準分布，其他3 個滿條件分布與之類似，可以采用Gibbs 抽樣得到全部參數(shù)（計 66 ×4 ＝264 個參數(shù)）的后驗樣本，形成 264 條MC鏈，完成Bayesian推斷.

2.2 模型的求解在OpenBUGS 環(huán)境下進行模型求解.由于參數(shù)個數(shù)較多，設定隨機數(shù)種子后，由軟件自動生成參數(shù)的迭代初始值.為降低后驗樣本間的自相關性、同時又確保馬氏鏈（MC 鏈）收斂，取間隔抽樣步長為100，燃燒期為5 000，進行105次迭代Gibbs抽樣，得到264 個待估參數(shù)的后驗樣本（264 條MC鏈）.扔掉燃燒期的前4 999 個樣本后，用剩余的95 001 個樣本通過MC 方法完成對參數(shù)θt的Bayesian推斷，相應的參數(shù)估計值（表1）包括：均值估計標準差（σ）、MC 誤差（εMC）、中位值估計（Median）及95% CI（可信區(qū)間）.同時將逐日 4 個參數(shù)（pt1、pt2、pt3及pt4各 66 個 Bayesian 估計）列表顯示，表中“…”表示中間的行省略.

表1 參數(shù)θt的Bayesian估計及部分統(tǒng)計量匯總表Tab. 1 Summary of Bayesian estimation and partial statistics of θt

2.3 模型診斷對于模型（3）中的264 個參數(shù)，分別繪制History-strace-plot圖、核密度估計曲線圖和 Autocorrelation 圖（以p1，1和p66，4為例，見圖 1，其余的262 個圖略，下同）. 264 個參數(shù)的 History -strace-plot 顯示，前4 999 個燃燒值樣本舍棄后，264 條MC鏈均收斂，且極限分布都是各自的后驗分布，適宜用相應的MC 鏈作為后驗樣本進行統(tǒng)計推斷.通過后驗樣本，分別求出參數(shù)pt，1至pt，4，t ＝ 1，2，…，66 的后驗核密度估計（圖 1 C -D，以p1，1和p66，4為例），結(jié)果顯示參數(shù)的分布呈對稱性，以正態(tài)分布為近似分布，后驗均值與中位值幾乎相等.每一條 MC 鏈的自相關圖（圖 1 E -F，以p1，1和p66，4為例）顯示，在滯后期≥2 以后，自相關系數(shù)接近于0，所以每一條鏈均可視為獨立同分布的MC 鏈，樣本均值成為一致、有效、無偏的估計. 264 個參數(shù)的εMC均小于5.37 × 10-5. pt1、pt2、pt3及pt4最大εMC對應的相應數(shù)據(jù)信息見表 2.表 2 中（max｛εMC｝／σ）×100%最大值僅為0.774 55%，說明模型（3）具有較高的精度［19］.最大的εMC發(fā)生的時點均是2020 -01 -23，最大εMC對應4 個概率全部集中在2020 -01 -23；而且在2020-02-13前，樣本計算得到的指標數(shù)據(jù)波動較大，后期逐漸光滑.說明2020 -02 -13前的日報數(shù)據(jù)誤差明顯大于后面的日報數(shù)據(jù)誤差，再次說明區(qū)間估計的重要性.由各自的可信區(qū)間計算得到pt1、pt2、pt3及pt4點估計的最大誤差限分別是0.013 535、0.027 7、0.011 695 及 0.030 7，預測誤差最大發(fā)生的時點也是2020 -01 -23.

圖 1 參數(shù)p1，1和p66，4的 History-strace-plot圖、核密度估計曲線和 Autocorrelation 圖Fig. 1 History-strace-plot，Kernel density estimation curve and Autocorrelation of parameters p1，1 and p66，4

表2 模型（3）的各條MC鏈的最大MC誤差及相應數(shù)據(jù)信息匯總表Tab. 2 Summary of maximum MC error and corresponding data information of each MC chain of Model（3）

采用概率意義下的可信區(qū)間對疫情指標進行估計，精準性更高.它不僅給出每個日報指標的下限，而且給出上限及誤差信息，同時，可信區(qū)間能覆蓋指標的直實數(shù)據(jù)值.為此，利用模型求解過程中得到了264 條MC鏈（每條鏈都收斂于各自的后驗分布），采用MC 方法得到相應參數(shù)的估計及統(tǒng)計量.然后提取逐日4 個概率的2.5%的后驗分位數(shù)和97.5%的分位數(shù)繪制各自95%的可信區(qū)間，從而得到逐日4 個概率的95%的可信區(qū)間帶.

3 結(jié)果與分析

3.1 治愈概率 、重癥概率、病死概率及輕癥概率估計本文2.3 模型診斷結(jié)果，充分說明模型（3）對于疫情的相對指標（4 個概率）估計的合理性與精確性，為了進一步將模型（3）的估計結(jié)果（概率）和樣本計算的結(jié)果（頻率）作對比分析.記fti（i ＝1，2，3，4）分別表示第 t 天樣本計算的治愈頻率、重癥頻率、病死頻率和輕癥頻率.第 t 天相應概率Dirichlet模型（3）計算得到的中位值估計（ptif），數(shù)值結(jié)果見表 1，分別用qtil和qtiu（i ＝1，2，3，4）分別表示概率ptif的95%的可信區(qū)間的下限和上限，并用計算結(jié)果繪制相應的頻率，概率曲線及95%可信區(qū)間帶（圖2）.

圖2 治愈頻率、重癥頻率、病死頻率、輕癥頻率及對應的后驗概率估計Fig. 2 Cure frequency，severe frequency，fatal frequency，mild frequency and the corresponding posterior probability estimation

結(jié)果顯示：1）同一天的4 個概率和為1；2）除病死概率的可信區(qū)間帶略寬，預測誤差略大外，其他三者的可信區(qū)間帶較窄，幾乎與預測值曲線重疊，精度較高；3）4 個概率曲線在02 -14 后均呈現(xiàn)出光滑變化特性；4）治愈概率、重癥概率、病死概率及輕癥概率基本上呈現(xiàn)3 個階段的變化特征.第一階段（2020 -01 -23 至02 -13）治愈概率經(jīng)歷短期下降后緩慢上升，在2020 -02 -11 呈現(xiàn)小幅振蕩；重癥概率在幅度為2.8% ～24.2%的范圍內(nèi)振蕩；病死概率在幅度為1.84% ～4.39%的范圍內(nèi)劇烈振蕩；輕癥概率在幅度為68% ～82%的范圍內(nèi)振蕩.第二階段（2020 -02 -14 至03 -26）治愈概率呈S型曲線逐漸上升，重癥概率在2020 -02 -16經(jīng)過小幅振蕩后逐步下降，前期（02 -14 至02 -21）下降速度慢度，后期（02 -21 以后）下降速度稍快；病死概率前期（02 -14 至02 -23）上升速度比后期（02 - 23 以后）快，變化范圍為 2. 18% ～4.18%，可信帶稍寬；輕癥概率呈反S 型曲線逐漸下降.第三階段（03 -26后），治愈概率上升到92%以上；因為疫情得到控制，大量治愈病人出院，醫(yī)院滯留病人快速減少，重癥概率降到1.2%以下，且呈減少趨勢；病死概率穩(wěn)定在4%左右；輕癥概率下降到3%以下，并趨于穩(wěn)定.

3.2 累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)估計為簡化公式，記

由Dirichlet分布的性質(zhì)可知

由K-Pearson替換原理，將各自的總體期望用后驗均值估計替代，可得第t 天總體總量指標的Bayesian估計，分別記為累計治愈病例數(shù)估計ctf，現(xiàn)有重癥病例數(shù)估計etf，累計病死病例數(shù)估計dtf，累記輕癥病例數(shù)估計（記為mtf）.并借助表1 中概率的可信區(qū)間上下限，進一步求各個總體總量指標的95%的可信區(qū)間（下限、上限分別在相應符號后加l、u表示.如：累記治愈數(shù)95%的可信區(qū)間記為95%CI：（ctl，ctu），下同）.將 4 個總量指標點估計及區(qū)間估計繪制曲線圖（見圖3）.

圖3 4 個總量指標估計及對應的95%可信區(qū)間Fig. 3 Four total quantity indexes and corresponding 95% confidence interval

結(jié)果顯示：

1）除累記病死病例的可信區(qū)間帶略寬，預測誤差略大外，其他三者的可信區(qū)間帶較窄，幾乎與估計值曲線重疊，精度較高，而且上述四條可信區(qū)間帶全部覆蓋了樣本值.

2）累計治愈人數(shù)和累計病死人數(shù)呈S 型曲線變化.表明醫(yī)院治愈病例快速增加，后期進入對重癥集中救治攻堅階段.兩個指標趨于穩(wěn)定，略呈上升趨勢.重癥病例數(shù)以2020 -02 -18 為分界，前期為增長階段，含 3 個振蕩時點（02 -07、02 -12、02 -16），間隔期分別為5 d.后期為下降階段，下降速度慢于前期的增長速度.輕癥病例數(shù)以2020 -02 -16 為分界，前期為增長階段，含一個跳躍時點（02 -12）.該值是第四版診療方案向第五版診療方案過渡，臨床病例并入確診病例產(chǎn)生.后期為下降階段，下降速度慢于前期的增長速度.

3）將圖3 中的4 組數(shù)據(jù)分別提取出來進行對比分析，發(fā)現(xiàn)絕大多數(shù)時點的Bayesian 點估計和樣本觀測值幾乎一樣，誤差也較小.

為了反映總量指標的預測精度，將累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)和累記輕癥病例數(shù)預測誤差限最大的4 個時點及對應的值提取出來匯總成表（表3），其他時點的預測誤差限均比該時點的誤差限小.

結(jié)果顯示：這4 個指標的點估計的最大預測誤差限為277.141 例，相應的指標是2020 -02 -25的輕癥病例數(shù)（樣本值36 852例）.

表3 4 個總量指標預測誤差限最大的時點及對應值Tab. 3 Time and corresponding values of the maximum prediction error limit of four total quantity indexes

66 d逐日4 個指標，共計264 個點估計值，528個區(qū)間估計值（下限，上限），無法在一個表中呈現(xiàn).以樣本截止日2020 -03 -28 的計算結(jié)果為例（表4）進行說明.結(jié)果顯示：點估計和樣本值幾乎相等，95%CI的上下限之差，點估計誤差限都很小，充分說明用可信區(qū)間進行疫情的估計更為合理.

表4 2020 -03 -28 的4 個總量指標Bayesian估計及誤差Tab. 4 Bayesian estimation and error of four aggregate indicators on March 2020 -03 -28

4 結(jié)論與展望

建立Bayesian框架下的Dirichlet 模型，借助樣本和先驗信息，通過Gibbs 抽樣方法得到總體參數(shù)的后驗樣本，采用MCMC方法計算得到264 個模型參數(shù)的Bayesian 估計，即逐日4 個概率指標（治愈概率、重癥概率、病死概率及輕癥概率）的點估計及區(qū)間估計；使用Dirichlet 分布的性質(zhì)進一步計算4個總量指標（累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死死病例數(shù)、現(xiàn)有輕癥病例數(shù)）的點估計及區(qū)間估計.在Dirichlet模型求解過程中，264 個參數(shù)對應264 條 MC 鏈均收斂，最大εMC僅為標準差的0.775%，且極限分布都是各自的后驗分布，適宜用相應的MC 鏈作為后驗樣本進行統(tǒng)計推斷.計算了2020 -01 -23 至03 -28 逐日4 個概率（治愈概率、病死概率、重癥概率、輕癥概率）的點估計（共264個值）和95%區(qū)間估計（共528 個值）；治愈概率、病死概率總體呈現(xiàn)上升后趨于平穩(wěn)趨勢，重癥概率、輕癥概率總體呈現(xiàn)逐漸下降趨勢；2020 -03 -26 后4 個概率指標都趨于平穩(wěn).計算了2020 -01 -23 至03 -28 逐日4 個總量指標（累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)）的點估計與95%區(qū)間估計.4 個總量指標點估計的最大預測誤差限為277.141 例，對應的指標是2020 -02 -25 的現(xiàn)有輕癥病例數(shù)（樣本值368 52例）.截至2020 -03 -28，累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)點估計誤差限分別為 146. 595、53. 100 8、109.948、85.515 例，對應的95%區(qū)間估計分別為：（75 301.20，75 594.39）、（689.881 0，796.082 6）、（3 189.937，3 409.833）、（1 864.045，2 035.075）.

用Bayesian框架下的Dirichlet模型計算逐日4個概率（治愈概率、病死概率、重癥概率、輕癥概率）的點估計與區(qū)間估計；計算逐日4 個總量指標（累計治愈病例數(shù)、現(xiàn)有重癥病例數(shù)、累計病死病例數(shù)、現(xiàn)有輕癥病例數(shù)）的點估計與區(qū)間估計，各自的計算精度均很高.Bayesian 方法還可廣泛運用于更廣泛的領域.