COVID-19潛伏期分布估計的統計學方法比較*

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

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

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

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

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

資料與方法

1.數據收集

本文數據來源于Lauer等[11]對COVID-19潛伏期分布估計的早期研究,該研究納入2020-01-04至2020-02-24中國湖北以外確診的181例COVID-19患者,這些患者的基本信息以及感染暴露和癥狀出現的時間區間均可從網絡新聞或公共衛生報告中獲取。

2.統計學方法

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

圖1 潛伏期觀測數據示意圖

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

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

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

我們關注的是潛伏期Ti的分布F(t),記S(t)=1-F(t)為Ti的生存函數。以下簡述基于區間刪失觀測的潛伏期分布估計方法(表1)。

表1 潛伏期估計方法匯總

(1)單區間刪失方法

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

(2)雙區間刪失方法

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

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

(3)隨機過程方法

圖2 COVID-19潛伏期估計的更新過程方法模型

3. 模型評價

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

4. 統計軟件

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

結 果

1. 基線特征

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

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

2. 潛伏期分布估計

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

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

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

圖4 基于雙區間刪失數據的COVID-19潛伏期分布估計

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

圖5 基于更新過程的COVID-19潛伏期分布估計

3.模型評價

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

表3 不同分布假設及分析方法獲得的COVID-19潛伏期估計結果

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

討 論

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

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

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

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

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