COVID-19潛伏期分布估計(jì)的統(tǒng)計(jì)學(xué)方法比較*

劉 裕 卓冰婷 陳俊宏 杜志成 郝元濤

【提 要】 目的 回顧和評(píng)估COVID-19潛伏期分布估計(jì)的統(tǒng)計(jì)學(xué)方法,為有效、快速、準(zhǔn)確地收集和分析潛伏期數(shù)據(jù)提供參考和借鑒。方法 利用COVID-19疫情早期發(fā)表的數(shù)據(jù),比較分析單區(qū)間刪失、雙區(qū)間刪失和隨機(jī)過(guò)程三類方法不同分布假設(shè)下獲得的COVID-19潛伏期分布最大似然估計(jì)和貝葉斯估計(jì)。結(jié)果 同類方法不同分布假設(shè)間,非參數(shù)方法要比參數(shù)方法擬合效果更好,但非參數(shù)方法存在較多的跳躍點(diǎn),且無(wú)法獲得估計(jì)的95%置信區(qū)間;同類方法相同分布假設(shè)條件下,最大似然估計(jì)與貝葉斯估計(jì)結(jié)果和擬合效果相近;同類方法的對(duì)數(shù)正態(tài)假設(shè)條件下獲得的潛伏期分布的大分位數(shù)(>90%分位數(shù))可能較大地偏離非參數(shù)估計(jì)結(jié)果;從數(shù)據(jù)利用的角度,雙區(qū)間刪失方法對(duì)數(shù)據(jù)的利用率最高;由于數(shù)據(jù)收集和利用的差異,不同方法得到的潛伏期分布估計(jì)可能存在較大差異。結(jié)論 采用雙區(qū)間刪失觀測(cè)的參數(shù)模型獲取傳染病潛伏期分布的最大似然估計(jì),可提高數(shù)據(jù)的收集、利用和分析效率;仔細(xì)比較不同分布假設(shè)下參數(shù)模型和非參數(shù)模型的結(jié)果,并謹(jǐn)慎解釋潛伏期大分位數(shù)的估計(jì)結(jié)果,將有利于作出正確的防控決策。

傳染病的潛伏期是指宿主首次暴露于傳染源到其首次出現(xiàn)疾病相關(guān)臨床表現(xiàn)(體征或癥狀)的時(shí)間間隔[1]。掌握潛伏期分布對(duì)病例的定義、傳染源的追溯、接觸者追蹤隨訪期和隔離期的設(shè)置、入境篩查隔離策略的制定、無(wú)癥狀感染人群醫(yī)學(xué)觀察期的確定等,以至疫情規(guī)模和傳播潛力的測(cè)算,都具有重要意義[2-4]。

然而,潛伏期分布的準(zhǔn)確估計(jì)并非易事,我們以新型冠狀病毒感染(COVID-19)為例加以說(shuō)明。首先,COVID-19的感染暴露時(shí)間無(wú)法直接觀測(cè),往往只能知道感染暴露是在某個(gè)時(shí)間段發(fā)生的,也就是說(shuō),它是一種區(qū)間刪失觀測(cè)(interval censored data)[5]。這也是盡管截至2020年1月22日已報(bào)告422例COVID-19確診患者,但Linton等[6]只納入10例具有明確暴露日期和發(fā)病日期的數(shù)據(jù)預(yù)估COVID-19潛伏期分布的可能原因。其次,感染以后患者出現(xiàn)癥狀的時(shí)間經(jīng)常不能準(zhǔn)確回憶,也就是說(shuō),患者出現(xiàn)癥狀(即發(fā)病)的時(shí)間也可能是區(qū)間刪失觀測(cè)。由于COVID-19疫情發(fā)生在冬季,疫情早期人們對(duì)該病知之甚少,其臨床癥狀與呼吸道感染重疊,COVID-19確診患者對(duì)首次出現(xiàn)新型冠狀病毒(SARS-Cov-2)感染癥狀的回憶往往摸棱兩可。此時(shí),調(diào)查得到的暴露感染和癥狀出現(xiàn)時(shí)間經(jīng)常都是區(qū)間刪失的情況,即我們獲得的是雙區(qū)間刪失觀測(cè)(doubly interval censored data)[7]。再者,對(duì)區(qū)間刪失尤其是雙區(qū)間刪失觀測(cè)數(shù)據(jù)的潛伏期分布估計(jì)遠(yuǎn)比精確觀測(cè)復(fù)雜,結(jié)果穩(wěn)定性也可能更差[8]。自2019年12月發(fā)生COVID-19疫情以來(lái),研究人員采用不同的統(tǒng)計(jì)學(xué)方法分析各自收集的數(shù)據(jù)估計(jì)COVID-19的潛伏期分布,得到其潛伏期中位數(shù)在4.0天到7.8天之間[6,9-13],相差較大。而他們獲得的COVID-19潛伏期大分位數(shù)估計(jì)差異更大,這體現(xiàn)在對(duì)潛伏期超過(guò)14天的患者比例的估計(jì)。例如,Bi等[9]估計(jì)潛伏期超過(guò)14天的COVID-19患者在5%左右,而Qin等[12]的結(jié)果顯示這個(gè)數(shù)值超過(guò)10%。

分布假設(shè)和估計(jì)方法對(duì)潛伏期的分布估計(jì)具有深刻影響[14]。為了加深傳染病潛伏期分布的理解,提升潛伏期分布監(jiān)測(cè)中數(shù)據(jù)收集和利用的效率,我們有必要對(duì)現(xiàn)有的分析模型進(jìn)行評(píng)估。本研究旨在綜述潛伏期分布估計(jì)方法,采用COVID-19疫情早期Lauer等[11]收集的數(shù)據(jù)對(duì)這些方法進(jìn)行比較,以期為有效、快速、準(zhǔn)確地預(yù)估潛伏期分布提供參考和借鑒。

資料與方法

1.數(shù)據(jù)收集

本文數(shù)據(jù)來(lái)源于Lauer等[11]對(duì)COVID-19潛伏期分布估計(jì)的早期研究,該研究納入2020-01-04至2020-02-24中國(guó)湖北以外確診的181例COVID-19患者,這些患者的基本信息以及感染暴露和癥狀出現(xiàn)的時(shí)間區(qū)間均可從網(wǎng)絡(luò)新聞或公共衛(wèi)生報(bào)告中獲取。

2.統(tǒng)計(jì)學(xué)方法

考慮包含n個(gè)獨(dú)立樣本的研究,假設(shè)樣本i(i=1,2,…,n)感染暴露和出現(xiàn)癥狀的時(shí)間分別為Ei和Oi(Oi>Ei),則該樣本的潛伏期為Ti=Oi-Ei。然而,在實(shí)踐中我們往往只知道感染暴露或癥狀出現(xiàn)落在某個(gè)可能的區(qū)間,也就是說(shuō),我們一般獲取如下形式的雙區(qū)間刪失觀測(cè)(圖1):

圖1 潛伏期觀測(cè)數(shù)據(jù)示意圖

Xi={(EiL,EiR],(OiL,OiR]}

其中,Ei∈(EiL,EiR],Oi∈(OiL,OiR],而且,EiL≤EiR,OiL≤OiR;特別地,當(dāng)區(qū)間的左端點(diǎn)與右端點(diǎn)相等時(shí)(EiL=EiR?Ei或OiL=OiR?Oi),表示觀測(cè)到的是確切的感染暴露或癥狀出現(xiàn)時(shí)間。如果能夠獲取Ei或Oi確切的觀測(cè)時(shí)間,則

Ti∈(OiL-Ei,OiR-Ei],Ti∈(Oi-EiR,Oi-EiL]或Ti=Oi-Ei。

我們關(guān)注的是潛伏期Ti的分布F(t),記S(t)=1-F(t)為Ti的生存函數(shù)。以下簡(jiǎn)述基于區(qū)間刪失觀測(cè)的潛伏期分布估計(jì)方法(表1)。

表1 潛伏期估計(jì)方法匯總

(1)單區(qū)間刪失方法

假定所有樣本的癥狀出現(xiàn)時(shí)間都是已知的,即對(duì)任意i,OiL=OiR?Oi,此時(shí),Ti∈(Oi-EiR,Oi-EiL]?(TiL,TiR]。這樣,潛伏期Ti的分布估計(jì)就簡(jiǎn)化為單個(gè)區(qū)間刪失數(shù)據(jù)的分析。令{sj}mj=0為{0,TiL,TiR:i=1,2,…,n}的唯一有序排列;記αij=I(sj∈(TiL,TiR])(I是示性函數(shù)),pj=F(sj)-F(sj-1),則似然函數(shù)可以表示為:

(2)雙區(qū)間刪失方法

類似地,如果假定潛伏期Ti服從某種特定的分布且可以表示成上述線性模型的形式,則我們同樣可以通過(guò)AFT模型來(lái)刻畫潛伏期的分布。令δi=I(EiL

其中,gφ和fθ分別為感染暴露時(shí)間和潛伏期的概率密度函數(shù),φ和θ分別為各自的分布參數(shù)。通常,假定感染暴露時(shí)間在觀測(cè)區(qū)間(EiL,EiR]均勻分布,這樣,我們最大化似然函數(shù)就可以獲得潛伏期分布參數(shù)θ的MLE估計(jì),從而得到潛伏期分布的估計(jì)。與單區(qū)間刪失方法一樣,我們也可以通過(guò)貝葉斯方法獲得潛伏期的分布估計(jì)。

(3)隨機(jī)過(guò)程方法

圖2 COVID-19潛伏期估計(jì)的更新過(guò)程方法模型

3. 模型評(píng)價(jià)

對(duì)于最大似然估計(jì),我們計(jì)算負(fù)對(duì)數(shù)似然函數(shù)值進(jìn)行同類方法內(nèi)的比較;同樣計(jì)算貝葉斯估計(jì)的負(fù)對(duì)數(shù)似然函數(shù)值,并與最大似然估計(jì)進(jìn)行比較。此外,我們對(duì)各種方法的數(shù)據(jù)利用情況及影響傳染病防控政策制定的潛伏期分位數(shù)估計(jì)(2.5%、25%、50%、75%、90%、95%、97.5%和99%分位數(shù))進(jìn)行仔細(xì)比較。

4. 統(tǒng)計(jì)軟件

本研究所有數(shù)據(jù)處理和建模過(guò)程均通過(guò)R軟件實(shí)現(xiàn)。其中,單區(qū)間刪失方法的NPMLE估計(jì)采用survival程序包,而MLE估計(jì)和Bayes估計(jì)采用icenReg程序包;雙區(qū)間刪失方法的NPMLE估計(jì)采用doubcens程序包,而MLE估計(jì)和Bayes估計(jì)采用coarseDataTools程序包;隨機(jī)過(guò)程方法基于Qin等[12]提供的R代碼實(shí)現(xiàn)。

結(jié) 果

1. 基線特征

研究數(shù)據(jù)來(lái)源于2020-01-04至2020-02-24中國(guó)湖北以外確診的COVID-19患者,總共181例。這些患者來(lái)自以亞洲為主的五大洲;年齡跨度較大,從2歲到80歲,平均年齡為46.0(±15.4)歲;108例(61.0%)為男性;159例(90.9%)有武漢旅居史,137例(75.7%)有明確的癥狀出現(xiàn)日期。具體信息見(jiàn)表2。

表2 研究對(duì)象的基本特征[n(%)]

2. 潛伏期分布估計(jì)

將137例具有明確癥狀出現(xiàn)日期的COVID-19患者納入單區(qū)間刪失方法分析。圖3的結(jié)果顯示,Turnbull的NPMLE存在較多的跳躍“階梯”;對(duì)于參數(shù)模型,相同的分布假設(shè)下,MLE估計(jì)與Bayes估計(jì)結(jié)果相近;盡管各模型對(duì)COVID-19的中位潛伏期的估計(jì)接近(5.4～6.0天),但對(duì)于大分位數(shù)(如>95%分位數(shù))的估計(jì)與Turnbull的NPMLE估計(jì)相比差別有變大趨勢(shì),置信區(qū)間變長(zhǎng),尤其是潛伏期的對(duì)數(shù)正態(tài)假設(shè)下,其MLE估計(jì)和Bayes估計(jì)與Turnbull的NPMLE估計(jì)差距最大,且99%分位數(shù)估計(jì)超過(guò)14天。

圖3 COVID-19潛伏期分布的單區(qū)間刪失方法分析

納入所有181例數(shù)據(jù)的雙區(qū)間刪失方法分析結(jié)果見(jiàn)圖4。可見(jiàn),NPMLE估計(jì)存在較多的“跳躍”點(diǎn);相同分布假設(shè)下的參數(shù)模型,其MLE估計(jì)與Bayes估計(jì)接近;相比之下,不同分布假設(shè)的參數(shù)模型估計(jì)的結(jié)果差別要大,估計(jì)的中位潛伏期在5.0～5.5天之間;對(duì)于潛伏期大分位數(shù)(如>95%分位數(shù))的估計(jì)與NPMLE估計(jì)差距變大,95%置信區(qū)間變寬,在對(duì)數(shù)正態(tài)假設(shè)下尤為明顯。這些結(jié)果都與單區(qū)間刪失方法得到的結(jié)果類似。

圖4 基于雙區(qū)間刪失數(shù)據(jù)的COVID-19潛伏期分布估計(jì)

從更新過(guò)程的角度,研究數(shù)據(jù)中包含59例2020-01-19至2020-01-21期間離開(kāi)武漢并在武漢以外確診且獲得確切癥狀出現(xiàn)日期(即前向復(fù)發(fā)時(shí)間明確)的患者,得到COVID-19潛伏期分布的MLE估計(jì)(圖5)。除威布爾分布假設(shè)下潛伏期分布的小于50%分位數(shù)估計(jì)明顯偏離其他兩種分布假設(shè)(對(duì)數(shù)正態(tài)分布和伽馬分布)外,其他各分位數(shù)估計(jì)接近,而且潛伏期中位數(shù)估計(jì)在4.0天左右。

圖5 基于更新過(guò)程的COVID-19潛伏期分布估計(jì)

3.模型評(píng)價(jià)

為了進(jìn)行模型間的比較,我們計(jì)算各模型擬合結(jié)果的負(fù)對(duì)數(shù)似然函數(shù)值。盡管Bayes估計(jì)目標(biāo)函數(shù)的優(yōu)化采用的是后驗(yàn)分布函數(shù),但本研究的結(jié)果顯示,相同分析方法和分布假設(shè)條件下,按潛伏期分布的Bayes估計(jì)計(jì)算得到的負(fù)對(duì)數(shù)似然函數(shù)值,略大于MLE估計(jì)的結(jié)果(表3),數(shù)值非常接近,提示Bayes估計(jì)與MLE估計(jì)吻合度很高。因此,這里僅比較不同模型的MLE估計(jì)。

表3 不同分布假設(shè)及分析方法獲得的COVID-19潛伏期估計(jì)結(jié)果

本研究的結(jié)果顯示,無(wú)論是單區(qū)間刪失方法、雙區(qū)間刪失方法,還是隨機(jī)過(guò)程的角度,各種方法不同分布假設(shè)條件下,其MLE估計(jì)的負(fù)對(duì)數(shù)似然函數(shù)值都非常接近,且都大于非參數(shù)方法。這提示,從擬合優(yōu)度的角度,非參數(shù)方法的結(jié)果優(yōu)于參數(shù)方法。如果我們以非參數(shù)模型結(jié)果為基準(zhǔn),無(wú)論是單區(qū)間刪失方法還是雙區(qū)間刪失方法,對(duì)數(shù)正態(tài)分布假設(shè)條件下的潛伏期大分位數(shù)(≥95%)估計(jì)更傾向于偏離非參數(shù)模型;而隨機(jī)過(guò)程方法在三個(gè)分布假設(shè)條件下的潛伏期大分位數(shù)估計(jì)基本一致。從數(shù)據(jù)利用的角度,由于受諸多假設(shè)條件的限制,隨機(jī)過(guò)程方法能夠利用的樣本數(shù)目(n=59)明顯少于單區(qū)間刪失方法(n=137)和雙區(qū)間刪失方法(n=181)。

討 論

本研究首先回顧了COVID-19潛伏期分布的統(tǒng)計(jì)估計(jì)方法,即單區(qū)間刪失方法,雙區(qū)間刪失方法和隨機(jī)過(guò)程方法,從收集數(shù)據(jù)的結(jié)構(gòu)、數(shù)學(xué)符號(hào)化過(guò)程到模型的構(gòu)建和實(shí)現(xiàn),以及模型的評(píng)價(jià),逐一進(jìn)行了詳細(xì)介紹;其次,利用Lauer等[11]收集的181例確診患者感染暴露和出現(xiàn)癥狀的信息,對(duì)三種方法的MLE估計(jì)和Bayes估計(jì)結(jié)果進(jìn)行了比較。我們的比較結(jié)果顯示,同類方法不同分布假設(shè)間,非參數(shù)方法要比參數(shù)方法擬合效果更好,但非參數(shù)方法存在較多的跳躍點(diǎn),且無(wú)法獲得估計(jì)的95%置信區(qū)間;同類方法相同分布假設(shè)條件下,MLE估計(jì)與Bayes估計(jì)結(jié)果和擬合效果相近;同類方法的對(duì)數(shù)正態(tài)假設(shè)條件下獲得的潛伏期分布的大分位數(shù)(>90%分位數(shù))可能較大地偏離非參數(shù)估計(jì)結(jié)果;從數(shù)據(jù)利用的角度,雙區(qū)間刪失方法對(duì)數(shù)據(jù)的利用率最高;由于數(shù)據(jù)收集和利用的差異,不同方法得到的潛伏期分布估計(jì)可能存在較大差異。

區(qū)間刪失數(shù)據(jù)的NPMLE估計(jì)被認(rèn)為是分析該類數(shù)據(jù)的金標(biāo)準(zhǔn)[7]。但非參數(shù)方法依賴于對(duì)潛伏期可能取值點(diǎn)的“猜測(cè)”,一般只能從樣本數(shù)據(jù)獲得,對(duì)于樣本數(shù)據(jù)以外的取值點(diǎn),在估計(jì)結(jié)果則體現(xiàn)為無(wú)信息的“水平線”或“線性插值”,這就是我們看到NPMLE存在較多“跳躍”點(diǎn)的原因(圖1A和圖3A)。另外,因?yàn)镹PMLE估計(jì)不需要任何的分布假設(shè)條件,從而無(wú)法進(jìn)行統(tǒng)計(jì)推斷,也就沒(méi)法計(jì)算估計(jì)的置信區(qū)間?；诖?研究人員普遍選擇的是潛伏期分布的參數(shù)模型估計(jì)[6,9-13]。然而,由于我們難于像非刪失數(shù)據(jù)估計(jì)方法那樣方便地檢查統(tǒng)計(jì)分布假設(shè)的準(zhǔn)確性(如殘差),我們完全有必要先獲得區(qū)間刪失數(shù)據(jù)的NPMLE,并將參數(shù)模型結(jié)果與之比較,只有在參數(shù)模型并未嚴(yán)重偏離NPMLE結(jié)果情況,才能有理由相信我們的參數(shù)模型結(jié)果的有效性和可靠性[16]。

在我們的研究里,同類方法相同分布假設(shè)條件下的MLE估計(jì)與Bayes估計(jì)結(jié)果和擬合效果相近。但是,一般模型的Bayes估計(jì),通常以MLE估計(jì)為初始估計(jì),采用模擬算法(如MCMC方法)通過(guò)最大化后驗(yàn)函數(shù)獲得。前期關(guān)于COVID-19潛伏期分布估計(jì),Backer等[10]和Linton等[6]利用stan語(yǔ)言[20]實(shí)現(xiàn),而且一般需要額外計(jì)算留一法交叉驗(yàn)證(leave-one-out cross validation,LOO-CV)或泛化信息量準(zhǔn)則(widely applicable information criterion,WAIC)參數(shù)[21]進(jìn)行模型比較,模型的收斂性有時(shí)難以保證。因此,盡管Bayes估計(jì)有其優(yōu)勢(shì)[22],但無(wú)論是從理論還是計(jì)算的復(fù)雜度而言,基于區(qū)間刪失數(shù)據(jù)的潛伏期分布Bayes估計(jì)不如其MLE估計(jì)直接和便捷。

不同分析方法之間,數(shù)據(jù)利用的效率差異較大,結(jié)果的變異也較大。顯然,基于雙區(qū)間刪失方法利用了所有收集的181例數(shù)據(jù),顯示了最高的數(shù)據(jù)利用效率。理論上,雙區(qū)間刪失方法對(duì)單區(qū)間刪失數(shù)據(jù)同樣適用,為此,我們采用雙區(qū)間刪失方法對(duì)137例單區(qū)間刪失數(shù)據(jù)重新進(jìn)行了分析,結(jié)果與單區(qū)間刪失方法完全一致。而Qin等[12]提出的隨機(jī)過(guò)程方法,雖然最終分析計(jì)算過(guò)程比較簡(jiǎn)單,而且在一定程度上可糾正數(shù)據(jù)收集過(guò)程中的回憶偏倚,但其假設(shè)條件較多,導(dǎo)致滿足條件的數(shù)據(jù)較少,從而產(chǎn)生樣本選擇偏倚,使其計(jì)算結(jié)果與區(qū)間刪失方法得到的估計(jì)差別較大。另外,在新發(fā)傳染病流行早期,數(shù)據(jù)采集和分析利用效率直接影響防控決策及其效果。因此,基于雙區(qū)間刪失數(shù)據(jù)分析方法是潛伏期的分布估計(jì)較好的選擇。

綜上所述,采用雙區(qū)間刪失數(shù)據(jù)的最大似然法估計(jì)傳染病潛伏期分布,可以提高數(shù)據(jù)的收集、利用和分析效率,減少樣本的選擇偏倚;潛伏期分布估計(jì)過(guò)程中,除了比較不同分布假設(shè)下的估計(jì)結(jié)果,還要與非參數(shù)模型估計(jì)進(jìn)行比較,并在不同數(shù)據(jù)集之間驗(yàn)證結(jié)果的可靠性;對(duì)潛伏期分布大分位數(shù)的估計(jì)和解釋要謹(jǐn)慎,僅依賴于模型擬合優(yōu)度統(tǒng)計(jì)量獲得的“最佳”估計(jì),有可能高估最長(zhǎng)潛伏期。