Joinpoint回歸模型及其在傳染病流行趨勢分析中的應用

廣東省疾病預防控制中心 廣東省公共衛(wèi)生研究院(511430)

曾四清

時間序列數(shù)據(jù)趨勢分析的經典方法包括移動平均模型、回歸模型、差分自回歸移動平均模型等,常用的回歸模型包括線性模型、指數(shù)模型、對數(shù)模型等。傳統(tǒng)回歸分析主要反映全局數(shù)據(jù)總體趨勢，可能無法揭示局部數(shù)據(jù)的特定趨勢。因此，分段回歸模型應運而生[1]，但如何分段又成為新的問題。Kim等提出的Joinpoint Regression(JPR)模型提供了解決方法[2]。近年來，JPR模型在癌癥和慢性病流行病學趨勢研究領域得到廣泛應用，主要用于分析發(fā)病率和死亡率時間趨勢變化特征，也用于分析其年齡趨勢、空間趨勢、疾病負擔，建立時空預測模型等[3-7]；也有少數(shù)用于傷害、傳染病等發(fā)病率和死亡率流行病學趨勢變化研究，以及文獻計量學研究等[8-12]。其分析結果可用于探討疾病流行趨勢變化規(guī)律及其影響因素，還可用于疾病預測、疾病負擔和干預效果評價等[3-12]。本文主要介紹JPR模型的基本原理和分析方法、應用條件、分析步驟和實現(xiàn)程序，并以2014年廣東省登革熱大流行疫情為例進行分析，探討該方法用于傳染病大流行疫情時間趨勢變化特征分析及預測方面的適用性，以便為相關專業(yè)人員正確選用該分析模型提供參考。

原理和方法

1．JPR模型的基本原理和分析方法

JPR模型分析的基本思想是通過模型擬合將一個長期趨勢線分成若干有統(tǒng)計學意義的趨勢區(qū)段，各段用連續(xù)的線性進行描述[2]。JPR模型分為線性數(shù)據(jù)模型和對數(shù)線性數(shù)據(jù)模型兩類[13]。

假設有一序列觀察值(x1,y1),…,(xn,yn),其中，x1≤…≤xn,線性數(shù)據(jù)模型表達式為：

E[yi|xi]=β0+β1xi+δ1(xi-τ1)++…+δk(xi-τk)+

(1)

對數(shù)線性數(shù)據(jù)模型表達式為：

E[yi|xi]=eβ0+β1 xi+δ1(xi-τ1)++…+δk( xi-τk)+

(2)

式中，yi(i=1，2，…，n)為因變量，xi(i=1，2，…，n)為自變量；β0表示不變參數(shù)，β1表示斜率參數(shù)(回歸系數(shù))；δk=βn+1,1-βn,1表示分段函數(shù)的回歸系數(shù)；τk為未知轉折點，k為轉折點個數(shù)，當(xi-τk)>0時，(xi-τk)+=(xi-τk)，否則，(xi-τk)+=0。

進行JPR分析時，要解決三個關鍵問題[13-14]，即選擇模型類別，確定轉折點數(shù)量、位置及模型參數(shù)，優(yōu)選模型。

(1)模型類別選擇基本原則和方法

如果系列數(shù)據(jù)樣本量在100以上，經過分布一致性檢驗，因變量服從正態(tài)分布或近似正態(tài)分布，宜選用模型(1)分析；如果因變量服從泊松分布或指數(shù)分布，則選用模型(2)分析[13-14]。如果分布不明確，宜同時采用模型(1)和(2)進行預分析，比較二者的均方差(mean squared errors,MSE)，選擇MSE較小者為優(yōu)選的擬合模型[13]。

(2)轉折點數(shù)量、位置及模型參數(shù)分析方法

采用網格搜索法(grid search method,GSM)，也可采用Hudson’s法，本文只介紹GSM。GSM是將參數(shù)所處的空間劃分成網格，每一個交點對應一個規(guī)劃方案，然后在設定的區(qū)間內以固定步長逐點計算對應方案的性能指標，以確定最優(yōu)參數(shù)[1]。進行JPR模型分析時，GSM創(chuàng)建“轉折點”全部可能位點的網格，計算每種情況下的誤差平方和(the sum of squared errors,SSE)及MSE，選擇MSE最小時的網格位點為最優(yōu)轉折點。JPR模型的回歸系數(shù)β0、β1、δ1、…、δk等參數(shù)估計值均采用GSM計算。

(3)優(yōu)選模型分析方法

采用Monte Carlo置換檢驗(permutation test)、貝葉斯信息準則(Bayesian information criterion,BIC)或修正BIC方法(modified BIC,MBIC)。

①置換檢驗：基本原理是：設定轉折點的最小值為MIN，最大值為MAX；無效假設H0設轉折點個數(shù)為ka，備選假設H1設轉折點個數(shù)為kb，其中ka

以線性數(shù)據(jù)模型為例說明如下[2]：假定某時間序列觀察值(x1,y1),…,(xn,yn)最多有2個轉折點。首先，無效假設 H0為沒有轉折點，即回歸模型E[yi|xi]=β0+β1xi；備選假設H1為有2個轉折點τ1和τ2(τ1<τ2)，即回歸模型E[yi|xi]=β0+β1xi+δ1(xi-τ1)++δ2(xi-τ2)+。如果經過統(tǒng)計學檢驗拒絕H0，再假設H0為有1個轉折點τ1，即E[yi|xi]=β0+β1xi+δ1(xi-τ1)+；H1不變。如果經過檢驗不拒絕H0，H0不變；再假設H1為有1個轉折點τ1，即E[yi|xi]=β0+β1xi+δ1(xi-τ1)+。這樣經過多次置換檢驗，最后得出有統(tǒng)計學意義的轉折點個數(shù)為1，回歸模型為E[yi|xi]=β0+β1xi+δ1(xi-τ1)+。

進行置換檢驗時，如果對所有的擬合模型都直接進行全數(shù)據(jù)集計算，計算量大、耗時長。因此，采用隨機數(shù)值生成器抽取Monte Carlo置換數(shù)據(jù)集樣本進行簡化計算，利用Bonferroni校正法設定多重檢驗的顯著性水平[14]。

②貝葉斯信息準則和修正貝葉斯信息準則[14]：有k個轉折點的JPR模型BIC值計算公式如下：

BIC(k)=ln{SSE(k)/#Obs}+{#Parm(k)/#Obs}×ln(#Obs),其中，SSE是誤差平方和，#Parm(k)=2×(k+1)是參數(shù)個數(shù)，#Obs是觀察值個數(shù)。

Zhang等[15]認為傳統(tǒng)BIC計算公式存在一些問題，提出其修正公式如下：

(3)

式中n為觀察記錄數(shù)，Γ(z)為gamma函數(shù)。

(4)

BIC值或MBIC值最小時的模型即為最優(yōu)模型。

(4)年度變化百分比和平均年度變化百分比及其可信區(qū)間的計算[14]

JPR模型中，年度變化百分比(annual percent change，APC)的計算公式如下，其中β為回歸系數(shù)。

APC={exp(β)-1}×100

(5)

(6)

(7)

JPR模型中，平均年度變化百分比(average annual percent change，AAPC)的計算公式如式(8)，其中wi為每個分段包括的年度數(shù)，βi為回歸系數(shù)。

(8)

(9)

(10)

如果轉折點為0，則AAPC=APC。

(5)預測精度評價

采用均方差(MSE)和平均相對誤差絕對值(mean absolute percentage error,MAPE)衡量預測的精確度。

(11)

(12)

2.JPR模型的應用條件和注意事項

(1)適用條件[14]

JPR模型主要用于時間序列數(shù)據(jù)或其他有序數(shù)據(jù)的趨勢變化特征分析。分析變量包括因變量、自變量，還可有分組變量。因變量數(shù)據(jù)類型包括：例數(shù)或“事件”數(shù)，如癌癥病例數(shù)、死亡人數(shù)等；粗發(fā)病率或死亡率、年齡標準化發(fā)病率或死亡率；構成比、百分比；或其他類型的數(shù)值變量，如氣溫值等。因變量可以在數(shù)據(jù)文件中直接提供(如發(fā)病數(shù)、死亡數(shù)等絕對數(shù))，也可以由數(shù)據(jù)文件提供的數(shù)據(jù)計算獲得(如發(fā)病率、標準化發(fā)病率、百分比、構成比等相對數(shù))。自變量應是有序變量，最常見的是時間，如年度、月份等；也可以是其他有序觀察值，如年齡等?？梢詫ψ宰兞窟M行重新編碼，但編碼必需是有序數(shù)值。若有分組變量，建立數(shù)據(jù)文件時必需先按分組變量進行排序才能正常分析。

如果需要比較兩組系列數(shù)據(jù)的分段回歸方程是否等價或趨勢一致性，可選用JPR模型中的“成對比較分析”(pairwise comparison)方法；如果由于疾病/死因編碼系統(tǒng)的改變等原因，引起監(jiān)測數(shù)據(jù)發(fā)生系統(tǒng)性的“躍變”，但并不影響所研究“事件”的潛在變化趨勢，可選用其中的“躍變”模型/可比性比例模型分析(jump model/comparability ratio(CR)model)方法[14，16]。

(2)注意事項

如果因變量觀察值是百分比，就不能大于100%；當有因變量觀察值為0時，不能直接計算，一般用“0.5”替換后再分析[14]。對因變量進行正態(tài)性檢驗時，Kolmogorov-Smirnov(K-S)檢驗法適用于樣本量大于100的連續(xù)型計量資料的分布一致性檢驗；當樣本量較小(小于100)時,不能直接使用其檢驗結果，應用Lilliefors修正臨界值表加以判斷[17]。

使用GSM時，需設定轉折點數(shù)的最小值(0)和最大值(≤9)。GSM允許最大轉折點數(shù)為9，但計算很耗時，一般設定不宜超過5；轉折點不可以是起始觀察值或末尾觀察值的位置；兩個轉折點之間至少有1個觀察值，轉折點之間的觀察值平均個數(shù)要在2個以上[14]。如果想得到一個轉折點，序列數(shù)據(jù)至少應有7個觀察值。

表1 觀察值個數(shù)和默認最大轉折點數(shù)關系表[14]

JPR模型是分析序列數(shù)據(jù)有顯著性意義變化的轉折點、變化方向和速度的理想方法，但并不表明模型的擬合度最優(yōu)[13]。在應用JPR模型時，其分析長期趨勢結果較為可靠，分析短期趨勢則會很大程度上受所分析的觀察值個數(shù)影響。如果以年度為單位進行趨勢分析，就無法揭示年度內短時間的變化特征和季節(jié)混雜因素的作用[8]，但可以通過改變時間維度進行深入分析；粗發(fā)病率(或死亡率)的趨勢變化不能直接反應人口構成對趨勢的影響，因此仍需通過標準化法排除人口構成變化引起的趨勢變化，或者按照年齡組分別進行趨勢分析；同時，分析長期趨勢時，還要考慮不同時期來源資料間的可比性問題[6]。

3.JRP軟件和JPR分析步驟[14]

(1)JRP軟件(joinpoint regression program)

JRP是美國國立癌癥研究所(National Cancer Institute)開發(fā)的基于Windows 操作系統(tǒng)進行JPR模型分析的專用軟件包，于2000年開始使用，2018年4月發(fā)布最新4.6.0.0版，可在“https://surveillance.cancer.gov/joinpoint/”網站免費下載安裝使用。

(2)分析步驟

使用JRP軟件進行JPR分析有五個基本步驟：創(chuàng)建數(shù)據(jù)文件、設置模型參數(shù)、運行程序、查看和選擇輸出結果，以及整理分析和解釋結果。數(shù)據(jù)文件須是ASCⅡ文本文件或excel文件?？梢允褂肧AS、SPSS、Excel、Word或其他軟件(如SEER*Stat Software)創(chuàng)建的文件。Excel文件須轉換為.CSV文件才可以使用。模型參數(shù)設置主要在“輸入文件”(input file)模塊設置因變量、自變量和模型類別等，在“方法和參數(shù)”(method and parameters)模塊設置模型擬合方法(GSM)、最大轉折點數(shù)和優(yōu)選模型分析方法(BIC、MBIC、permutation test)等。在高級分析(advanced analysis tools)模塊，可選擇進行趨勢“成對比較分析”(pairwise comparison)或“躍變”模型/可比性比例模型分析(jump model/CR Model)。

實例分析

以2014年廣東省登革熱大流行疫情為例，進行周發(fā)病數(shù)的趨勢變化特征分析及預測。

1.數(shù)據(jù)來源

資料來源于國家“傳染病報告信息管理系統(tǒng)”，各周發(fā)病數(shù)見表2，全年發(fā)病數(shù)共計45188例。

2.分析過程及主要結果

采用excel 2013建立Dg2014.xlsx數(shù)據(jù)表和Dg2014.CSV數(shù)據(jù)文件，將其中“發(fā)病數(shù)”為“0”的值用“0.5”替代；采用SPSS 19.0軟件對Dg2014.xlsx周發(fā)病數(shù)進行Lilliefors修正K-S檢驗，統(tǒng)計量Z=0.377，df=52，P<0.01，表明周發(fā)病數(shù)不服從正態(tài)分布，故采用JRP軟件選擇模型(2)對Dg2014.CSV周發(fā)病數(shù)進行趨勢分析，最大轉折點數(shù)設置為5。采用permutation test和MBIC進行優(yōu)選模型分析時MSE分別為22.30(有2個轉折點)、17.20(有3個轉折點)。后者值較小，擬合精度優(yōu)于前者，因此，采納MBIC分析結果。擬合回歸方程為：E[yi|xi]=e-0.757+0.055xi+0.438(xi-22)+-0.971(xi-40)+-0.213(xi-43)+。yi為周發(fā)病數(shù)，xi為周數(shù)，i=1～52，xi=1～52，如果(xi-22)+>0，則(xi-22)+=xi-22,如果(xi-22)+≤0，則(xi-22)+=0；如果(xi-40)+>0，則(xi-40)+=xi-40,如果(xi-40)+≤0，則(xi-40)+=0；如果(xi-43)+>0，則(xi-43)+=xi-43,如果(xi-43)+≤0，則(xi-43)+=0；MBIC=1.468。有3個有統(tǒng)計學意義的轉折點將周發(fā)病數(shù)分成了4個趨勢區(qū)段。轉折點分別在第22周(95%CI:7～33)、第40周(95%CI:39～41)和第43周(95%CI:42～46)。在第1～22周期間，APC=5.64%(95%CI:-10.6%～24.8%，P>0.05)；在第22～40周期間，APC=63.76%(95%CI:61.0%～66.6%，P<0.01)；在第40～43周期間，APC=-37.96%(95%CI:-42.8%～-32.7%，P<0.01)；在第43～52周期間，APC=-49.88%(95%CI:-52.9%～-46.6%，P<0.01)。結果見圖1。

根據(jù)擬合方程進行周發(fā)病數(shù)預測，結果見表2，平均相對誤差絕對值MAPE為17.05%。

討 論

JPR模型主要應用于時間序列數(shù)據(jù)或其他有序數(shù)據(jù)的趨勢變化特征分析，主要采用網格搜索法、置換檢驗和貝葉斯信息準則等數(shù)理運算法則，進行疾病流行趨勢變化擬合模型篩選和參數(shù)估計，分析趨勢改變的轉折點數(shù)量、確定轉折點位置、明確變化方向、計算變化速度等特征[2,14]。因此，比僅根據(jù)常用數(shù)據(jù)列表和統(tǒng)計圖直觀描述疾病流行趨勢變化更加科學。JRP軟件為JPR模型運算提供了簡便適用的程序包。

圖1 2014年廣東省登革熱周發(fā)病數(shù)趨勢變化JPR分析圖

2014年廣東省發(fā)生了登革熱大流行疫情，總報告病例數(shù)超過4.5萬例[18-19]?？茖W分析了解傳染病大流行疫情的趨勢變化特征和規(guī)律，將有利于相關部門及時采取有效的應對策略和干預措施，更好預防控制疫情發(fā)生發(fā)展。目前，國內對于登革熱的流行病學研究，大多采用傳統(tǒng)流行病學方法描述其三間分布[18-19]；也有少數(shù)采用空間自相關和時空掃描聚類等分析方法進行病例時空聚集性研究，或采用SARIMA 模型及其與GRNN的組合模型預測登革熱的月發(fā)病數(shù)等[20-21]。因此，采用JPR模型分析2014年廣東省登革熱疫情發(fā)生發(fā)展的時間趨勢變化特征具有特殊意義。

結果顯示，2014年廣東省登革熱大流行時的周發(fā)病數(shù)擬合JPR模型的MBIC、MSE和MAPE值均較小，而且在爆發(fā)期的預測精度較高，說明擬合模型基本有效，此類疫情對JPR模型擬合有較好的適用性。登革熱周發(fā)病數(shù)的趨勢變化有3個非常關鍵的轉折點，即第22周(5月26日-6月1日)、第40周(9月29日-10月5日)和第43周(10月20-26日)。登革熱的潛伏期為1～15d[22]，分別將以上3個轉折點時間往前推移一個最長潛伏期(按2周計)，即是第20周(5月12日-18日)、第38周(9月15-21日)、第41周(10月6-12日)，此三個時段應當是2014年登革熱防控工作的關鍵時期。結果表明，如果不采取“早預防”措施，即使早期病例不多，也可能隱藏較大的疫情大流行風險；如果發(fā)展到疫情快速上升期，周發(fā)病數(shù)平均增長速度可達很高(APC超過60%)，疫情進展將十分迅猛，極易形成爆發(fā)，控制難度加大，時間延續(xù)較長；經過采取全面有效的控制措施以后，疫情下降速度也較快(APC超過-30%)。而且，隨著各項防控措施的持續(xù)全面深入開展，效果進一步顯現(xiàn)，疫情下降速度加快(APC超過-40%)?？梢姡琂PR模型對于揭示登革熱大流行時有顯著性意義的時間趨勢變化特征，包括關鍵控制期、疫情變化方向和速度等是較理想的分析方法。雖然該模型的擬合度并不一定最優(yōu)[13]，但可以用于疫情趨勢分析和預測，從而為防控工作提供重要參考。

表2 2014年廣東省登革熱各周發(fā)病數(shù)及擬合值比較(例)

*：實際發(fā)病數(shù)為0，不能計算相對誤差，也未列入平均相對誤差絕對值計算。

JPR模型是一種較新的趨勢分析方法，具有傳統(tǒng)回歸分析方法不具備的某些分析功能和特點[13，23]；它可以用于分析某些傳染病大流行期間的時間趨勢變化特征和預測，其分析結果能為防控工作提供決策參考。因此，它在各種傳染病流行趨勢變化特征分析和預測中的應用價值值得深入探討。在應用時，要充分了解其基本原理、分析方法、適用條件和注意事項，科學解釋并合理運用其分析結果為傳染病防控工作服務。