兩組序列數(shù)據(jù)趨勢(shì)變化特征對(duì)比分析的Joinpoint回歸模型方法及應(yīng)用

曾四清 李 艷 劉 珺 謝仕蘭 付笑冰 龍其穗 林 鵬

【提 要】 目的 介紹兩組序列數(shù)據(jù)趨勢(shì)變化特征對(duì)比分析的Joinpoint回歸(JPR)模型的基本原理和方法，并以2006-2018年廣東省HIV/AIDS病例數(shù)為實(shí)例進(jìn)行分析。方法 該模型原理和方法主要來(lái)源于相關(guān)文獻(xiàn)和應(yīng)用軟件，實(shí)例數(shù)據(jù)來(lái)源于國(guó)家艾滋病綜合防治信息系統(tǒng)，采用Excel 2013建立數(shù)據(jù)文件，采用JRP軟件進(jìn)行分析。結(jié)果 2006-2018年廣東省新發(fā)現(xiàn)的經(jīng)注射毒品感染和經(jīng)性傳播的HIV/AIDS病例數(shù)分別擬合對(duì)數(shù)線性JPR模型的MSE為21.44和16.63，MAPE為9.51%和4.57%；二者的回歸方程不平行(P<0.001)。其中，經(jīng)注射毒品感染的HIV/AIDS病例數(shù)在2008年出現(xiàn)1個(gè)轉(zhuǎn)折點(diǎn)，其前后分段的APC分別為12.64%和-17.75%，AAPC為-13.31%，呈先升后降和總體下降趨勢(shì)；而經(jīng)性傳播的HIV/AIDS病例數(shù)在2008年和2015年各出現(xiàn)1個(gè)轉(zhuǎn)折點(diǎn)，其前后分段APC分別為101.47%、22.26%和1.91%，AAPC為27.02%，呈持續(xù)上升趨勢(shì)，但上升速度逐步下降，2015年后趨于平緩。兩組的AAPC差異有統(tǒng)計(jì)學(xué)意義(z=-9.50，P=0.00)，總體趨勢(shì)變化方向及速度不一致。結(jié)論 該模型的應(yīng)用擬合效果較好，兩組數(shù)據(jù)擬合的回歸方程及其趨勢(shì)變化特征指標(biāo)均存在顯著統(tǒng)計(jì)學(xué)差異。研究表明，該模型適用于某些傳染病序列數(shù)據(jù)趨勢(shì)變化特征量化對(duì)比分析。

非參數(shù)趨勢(shì)檢驗(yàn)適用于計(jì)數(shù)資料，只能得出趨勢(shì)結(jié)論，無(wú)法得到量化趨勢(shì)指標(biāo)，不同組別數(shù)據(jù)之間無(wú)法比較趨勢(shì)變化大??；線性回歸適用于計(jì)數(shù)資料和計(jì)量資料，有豐富的衍生模型，還可控制混雜因素，可計(jì)算趨勢(shì)變化量化指標(biāo)，不同組別數(shù)據(jù)之間可以比較趨勢(shì)變化大小[1-4]。Joinpoint回歸(joinpoint regression，JPR)是一種特殊的線性回歸衍生模型，除了可以進(jìn)行趨勢(shì)變化參數(shù)估計(jì)量化分析外，還可以分析趨勢(shì)變化的轉(zhuǎn)折點(diǎn)，從而進(jìn)行分段趨勢(shì)變化特征分析；還可以進(jìn)行兩組序列數(shù)據(jù)趨勢(shì)變化特征對(duì)比分析[5-6]。然而，目前有關(guān)該模型在傳染病趨勢(shì)分析中的應(yīng)用研究文獻(xiàn)較少。因此，本文介紹此模型的基本原理和方法，并以2006-2018年廣東省HIV/AIDS病例數(shù)為實(shí)例進(jìn)行分析，探討其應(yīng)用價(jià)值。

原理和方法

1.模型及計(jì)算原理

JPR模型分析的基本思想是通過(guò)模型擬合將一個(gè)長(zhǎng)期趨勢(shì)線分成若干有統(tǒng)計(jì)學(xué)意義的趨勢(shì)區(qū)段，各段用連續(xù)的線性進(jìn)行描述[5]。JPR模型分為線性數(shù)據(jù)模型和對(duì)數(shù)線性數(shù)據(jù)模型兩類。假設(shè)有一序列觀察值(x1，y1)，…，(xn，yn)，其中，yi(i=1，2，…，n)為因變量，xi(i=1，2，…，n)為自變量，x1≤…≤xn。線性數(shù)據(jù)模型表達(dá)式為：

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

(1)

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

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

(2)

式中，β0表示不變參數(shù)，β1表示斜率參數(shù)(回歸系數(shù))；δk=βn+1-βn表示分段函數(shù)的回歸系數(shù)；τk為未知轉(zhuǎn)折點(diǎn)，k為轉(zhuǎn)折點(diǎn)個(gè)數(shù)，當(dāng)(xi-τk)>0時(shí)，(xi-τk)+=(xi-τk)，否則，(xi-τk)+=0。

(3)

(4)

(1)兩組序列數(shù)據(jù)的JPR方程(回歸參數(shù))比較[6-7]

假設(shè)(xij，yij)是i組序列觀察值(i=1，2；j=1，2，…，n)，yij為因變量，xij為自變量。假設(shè)xij=x，yij的回歸方程為：

μi(x)=βi，0+βi，1x+δi，1(x-τi，1)++…+δi，ki(x-τi，ki)+

(5)

根據(jù)公式(1)和公式(2)，對(duì)于線性數(shù)據(jù)模型，μi(x)=yi(x)；對(duì)于對(duì)數(shù)線性數(shù)據(jù)模型，μi(x)=lnyi(x)。

檢驗(yàn)1：比較兩個(gè)回歸方程是否一致(包括模型參數(shù)截距、轉(zhuǎn)折點(diǎn)數(shù)、轉(zhuǎn)折點(diǎn)位置和回歸系數(shù))。

設(shè)模型參數(shù)(ki，θi)=(ki，τi，1，…，τi，ki，βi，0，βi，1，δi，1，…，δi，ki)

無(wú)效假設(shè)H0：(k1，θ1)=(k2，θ2)

備選假設(shè)H1：(k1，θ1)≠(k2，θ2)

檢驗(yàn)2：比較兩個(gè)回歸方程是否是不同截距的平行方程(包括模型參數(shù)轉(zhuǎn)折點(diǎn)數(shù)、轉(zhuǎn)折點(diǎn)位置和回歸系數(shù)，但不包括截距)。

檢驗(yàn)1和檢驗(yàn)2的統(tǒng)計(jì)量計(jì)算公式為：

(6)

(2)兩組序列數(shù)據(jù)JPR模型的APCi和AAPC比較[6，8]

APCi(annual percent change)是指序列數(shù)據(jù)的年度變化百分比，反映階段性變化速度；AAPC是指序列數(shù)據(jù)全局或選定區(qū)段的年度變化百分比的加權(quán)平均值，反映平均變化速度。對(duì)數(shù)線性模型才可計(jì)算APC和AAPC并進(jìn)行比較。

APCi及AAPC計(jì)算公式如下：

APCi={exp(βi)-1}×100

(7)

基于t分布假設(shè)計(jì)算APCi的100(1-α)%可信區(qū)間，公式為：

APCiL(α)=100{exp(βi-si*t-1di(1-α/2))-1}

(8)

APCiU(α)=100{exp(βi+si*t-1di(1-α/2))-1}

(9)

(10)

基于標(biāo)準(zhǔn)正態(tài)分布假設(shè)計(jì)算AAPC的100(1-α)%可信區(qū)間，公式為：

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

2.分析方法和步驟

采用Excel 2013建立序列數(shù)據(jù)的.CSV格式文件，采用JRP(Joinpoint regression program，JRP)軟件4.6.0.0版[6]進(jìn)行分析。采用網(wǎng)絡(luò)搜索法(grid search method，GSM)分析轉(zhuǎn)折點(diǎn)數(shù)量、位置及模型參數(shù)，采用Monte Carlo置換檢驗(yàn)(permutation test，PT)和修正貝葉斯信息準(zhǔn)則(modified Bayesian information criterion，MBIC)優(yōu)選模型，選用其中MSE較小者的分析結(jié)果；采用置換檢驗(yàn)方法進(jìn)行兩組模型回歸方程比較，采用標(biāo)準(zhǔn)正態(tài)分布原理進(jìn)行兩組AAPC差異顯著性檢驗(yàn)。

分析時(shí)，在JRP 軟件的“輸入文件”(Input file)窗口設(shè)置模型參數(shù)，設(shè)HIV/AIDS病例數(shù)為因變量(dependent variable)、年度為自變量(independent variable)、感染途徑為分組變量(By variable)，模型類別選擇為對(duì)數(shù)線性數(shù)據(jù)模型[ln(y)=xb]；假設(shè)HIV/AIDS病例數(shù)服從Poisson分布，將異方差選項(xiàng)(Heteroscedastic errors option)設(shè)置為“Poisson variance”；在“方法和參數(shù)”(Method and parameters)窗口設(shè)置模型擬合方法為網(wǎng)格搜索法(Grid search)、置換檢驗(yàn)(Permutation test)，或者為修正貝葉斯信息準(zhǔn)則(Modified BIC)，經(jīng)試分析本研究設(shè)最大轉(zhuǎn)折點(diǎn)數(shù)(Number of joinpoints)kmax值為2，誤差自相關(guān)選項(xiàng)(Autocorrelated errors options)設(shè)置為誤差不相關(guān)模型(Fit an uncorrelated errors model)，選用參數(shù)法(Parametric)計(jì)算AAPC；在“高級(jí)分析”(Advanced analysis tools)窗口，選擇對(duì)比分析(Pairwise comparison)方法[6]。

3.應(yīng)用條件和注意事項(xiàng)

(1)適用條件

經(jīng)典線性回歸模型基本適用條件包括：因變量y與自變量x呈線性關(guān)系；觀察值相互獨(dú)立；在一定范圍內(nèi)，任意給定x值，因變量y值都服從正態(tài)分布；不同x值所對(duì)應(yīng)y值的方差相等。JPR模型是線性回歸分析的一種衍生模型。對(duì)比分析JPR模型的因變量類型包括例數(shù)或“事件”數(shù)、發(fā)病率或死亡率、構(gòu)成比、百分比，或其他數(shù)值變量等；自變量應(yīng)是有序變量；同時(shí)，還應(yīng)有分組變量，建立數(shù)據(jù)文件時(shí)需按分組變量排序。若因變量y不完全符合以上條件，則需對(duì)其進(jìn)行數(shù)據(jù)轉(zhuǎn)換(如對(duì)數(shù)轉(zhuǎn)換)以達(dá)到適用條件，此時(shí)應(yīng)采用對(duì)數(shù)線性模型[2]。當(dāng)有因變量y值為“0”時(shí)，需用“0.5”替換后才能分析[6]。

(2)注意事項(xiàng)

對(duì)比分析JPR模型既可以是線性方程比較，也可以是對(duì)數(shù)線性方程比較，但如果要計(jì)算并比較APCi和AAPC，則應(yīng)選擇對(duì)數(shù)線性模型。

在JRP軟件中，單一序列數(shù)據(jù)的JPR趨勢(shì)分析既可選擇非相關(guān)誤差模型(模型中的隨機(jī)誤差具有方差齊同性或固定方差，其符合均值為0、方差為δ2的正態(tài)分布)，也可選用自相關(guān)誤差模型(即隨機(jī)誤差的方差隨時(shí)間而變化)；而進(jìn)行兩組對(duì)比分析時(shí)，只能采用非相關(guān)誤差模型[6-7]。

進(jìn)行置換檢驗(yàn)時(shí)，需設(shè)定檢驗(yàn)次數(shù)，取值范圍為1000～10000。為了計(jì)算結(jié)果(P值)穩(wěn)定，同時(shí)控制計(jì)算時(shí)間，一般設(shè)定至少為4499，JRP中的默認(rèn)值為4499；還需設(shè)定最大轉(zhuǎn)折點(diǎn)數(shù)kmax值，由于置換檢驗(yàn)較耗時(shí)，且kmax過(guò)大反而降低檢驗(yàn)效能，一般設(shè)定kmax≤4[7]。

通常AAPC不能代表APCi，但當(dāng)分段數(shù)為0時(shí)，APCi即為AAPC。此時(shí)，AAPC的可信區(qū)間采用t分布修正，以保證其與APCi的可信區(qū)間一致[6]；一般計(jì)算AAPC的區(qū)段是全局的，但其起始點(diǎn)或截止點(diǎn)也可以在全局范圍內(nèi)自行設(shè)定，應(yīng)用靈活。

對(duì)比分析JPR模型能根據(jù)分組變量對(duì)多組序列數(shù)據(jù)進(jìn)行兩兩比較分析，但不能同時(shí)進(jìn)行3組及以上組的比較分析；在進(jìn)行多組別之間的比較時(shí)，是否需要將相鄰兩組的序列數(shù)據(jù)進(jìn)行合并分析，可根據(jù)兩者比較的結(jié)果來(lái)判定。該對(duì)比分析方法不能用于只比較兩組序列數(shù)據(jù)的轉(zhuǎn)折點(diǎn)位置而不考慮其他回歸參數(shù)的分析[7]。

實(shí)例分析

實(shí)例數(shù)據(jù)來(lái)源于國(guó)家艾滋病綜合防治信息系統(tǒng)。2006-2018年廣東省新發(fā)現(xiàn)經(jīng)注射毒品感染與經(jīng)性傳播的HIV/AIDS病例數(shù)見(jiàn)表1，在此擬合對(duì)比分析JPR模型，進(jìn)行兩組病例數(shù)的年度趨勢(shì)變化特征對(duì)比分析。

表1 2006-2018年廣東省新發(fā)現(xiàn)經(jīng)注射毒品感染與經(jīng)性傳播的HIV/AIDS病例數(shù)

1.兩組病例數(shù)的年度趨勢(shì)變化回歸方程比較

2006-2018年新發(fā)現(xiàn)的經(jīng)注射毒品感染和經(jīng)性傳播HIV/AIDS病例數(shù)擬合線性數(shù)據(jù)模型的MSE分別為47.37、25.63，擬合對(duì)數(shù)線性數(shù)據(jù)模型的MSE分別為21.44、16.63，均是后者小于前者，故采用對(duì)數(shù)模型進(jìn)行對(duì)比JPR模型分析；置換檢驗(yàn)、修正貝葉斯信息準(zhǔn)則優(yōu)選模型方程時(shí)，經(jīng)注射毒品感染HIV/AIDS病例數(shù)擬合模型的MSE均為21.44，經(jīng)性傳播病例數(shù)擬合模型的MSE均為16.63，二者結(jié)果一致。兩組的趨勢(shì)變化特征比較見(jiàn)圖1，有關(guān)參數(shù)分析結(jié)果比較見(jiàn)表2。經(jīng)注射毒品感染的HIV/AIDS病例數(shù)回歸方程為：E[yi|xi]=e-231.27+0.12xi-0.31(xi-2008)+，經(jīng)性傳播的HIV/AIDS病例數(shù)回歸方程為：E[yi|xi]=e-1398.78+0.70xi-0.50(xi-2008)+-0.18(xi-2015)+。其中，xi為年度，取值為2006，2007，…，2018；E[yi|xi]為HIV/AIDS病例數(shù)。二者比較差異有統(tǒng)計(jì)學(xué)意義(P<0.001)。說(shuō)明兩組擬合的JPR回歸方程不一致、不平行。

表2 2006-2018年廣東省新發(fā)現(xiàn)經(jīng)注射毒品感染與經(jīng)性傳播HIV/AIDS病例數(shù)回歸方程參數(shù)

圖1 2006-2018年廣東省新發(fā)現(xiàn)經(jīng)注射毒品感染與經(jīng)性傳播HIV/AIDS病例數(shù)趨勢(shì)比較

2.APCi及AAPC比較

經(jīng)注射毒品感染與經(jīng)性傳播的HIV/AIDS病例數(shù)擬合對(duì)比分析JPR模型的趨勢(shì)區(qū)段、APCi及AAPC結(jié)果比較見(jiàn)表3。AAPC(1)-AAPC(2)=-40.30%(95%CI：-48.62%～-32.01%)；z=-9.50，P<0.001。說(shuō)明兩組的AAPC差異有統(tǒng)計(jì)學(xué)意義，即兩組的年平均變化速度差異有統(tǒng)計(jì)學(xué)意義。

表3 2006-2018年廣東省經(jīng)注射毒品感染與經(jīng)性傳播HIV/AIDS病例數(shù)趨勢(shì)對(duì)比分析APCi(%)及AAPC(%)結(jié)果

討 論

本文主要介紹了兩組序列數(shù)據(jù)趨勢(shì)變化特征對(duì)比分析的JPR模型的基本原理、計(jì)算方法、應(yīng)用條件及注意事項(xiàng)。同時(shí)，以2006-2018年廣東省新發(fā)現(xiàn)的經(jīng)注射毒品感染與經(jīng)性傳播的HIV/AIDS病例數(shù)為實(shí)例進(jìn)行了趨勢(shì)對(duì)比分析。

JPR回歸是分析某些疾病發(fā)病數(shù)或發(fā)病率序列數(shù)據(jù)趨勢(shì)變化特征的實(shí)用方法[5，9-10]。JPR回歸分析擬合的模型不一定最佳，但其優(yōu)勢(shì)在于能計(jì)算AAPC，反映整體的平均趨勢(shì)變化方向和速度，也可應(yīng)用于同一種類序列數(shù)據(jù)的不同組別之間進(jìn)行趨勢(shì)變化平均速度差異的量化比較，這是非參數(shù)趨勢(shì)檢驗(yàn)所不具備的功能[1-4]。而且，該模型能檢測(cè)出有顯著性意義的趨勢(shì)變化轉(zhuǎn)折點(diǎn)，并計(jì)算轉(zhuǎn)折點(diǎn)前后各區(qū)段的趨勢(shì)變化方向和速度APCi，從而揭示趨勢(shì)變化轉(zhuǎn)折點(diǎn)數(shù)量、位置和階段性趨勢(shì)變化指標(biāo)特征[11]。同時(shí)APCi由于與AAPC一樣，是一個(gè)無(wú)量綱的相對(duì)數(shù)，是反映趨勢(shì)變化方向和速度的指標(biāo)，因此，采用該模型進(jìn)行時(shí)間序列發(fā)病數(shù)據(jù)趨勢(shì)分析時(shí)，相同類別的不同病種之間、相同病種的不同時(shí)段之間或者不同組別人群之間等發(fā)病數(shù)據(jù)的APCi就可以進(jìn)行比較分析。這些都是經(jīng)典回歸分析所不具備的分析功能[2]。正因如此，本文探討了其在不同組別人群的傳染病時(shí)間趨勢(shì)對(duì)比分析中的應(yīng)用。

兩組的JPR模型擬合精確度均較高，兩組的JPR回歸方程比較差異有統(tǒng)計(jì)學(xué)意義，不一致、不平行。說(shuō)明需要用不同的回歸方程才能描述其各自的年度趨勢(shì)變化特征。同時(shí)，兩組的年平均變化方向和速度差異有統(tǒng)計(jì)學(xué)意義，說(shuō)明二者的總體趨勢(shì)變化也存在著差異，前者呈總體下降趨勢(shì)，而后者呈總體上升趨勢(shì)。其中，經(jīng)注射毒品感染的HIV/AIDS病例數(shù)在2008年開(kāi)始發(fā)生了顯著性轉(zhuǎn)折變化，從此前的上升趨勢(shì)轉(zhuǎn)變?yōu)榇撕蟮某掷m(xù)穩(wěn)定下降趨勢(shì)；提示廣東省2005年開(kāi)始開(kāi)展的監(jiān)管場(chǎng)所吸毒人群艾滋病篩查和2006年開(kāi)始開(kāi)展的社區(qū)戒毒藥物維持治療等防控措施取得了較好的成效。經(jīng)性傳播的HIV/AIDS病例數(shù)從2006年到2018年總體呈顯著上升趨勢(shì)，但在2008年和2015年開(kāi)始出現(xiàn)了兩次上升速度的統(tǒng)計(jì)學(xué)顯著性變化，上升速度呈階段性減低的特點(diǎn)，至2015年以后漸趨平穩(wěn)，也反映了廣東省采取的HIV/AIDS經(jīng)性傳播防控策略與措施產(chǎn)生了一定的積極效果。研究表明，該對(duì)比分析JPR模型適用于不同感染途徑的HIV/AIDS病例數(shù)時(shí)間序列趨勢(shì)變化特征比較研究，對(duì)比分析其趨勢(shì)變化的轉(zhuǎn)折點(diǎn)位置、方向和速度。參照上述方法，就可以對(duì)各種感染途徑或者按其他因素分組的HIV/AIDS新發(fā)病例數(shù)、發(fā)病率或構(gòu)成比等的趨勢(shì)變化特征進(jìn)行對(duì)比分析，以全面揭示其變化特征之規(guī)律及異同，進(jìn)而對(duì)其背后的原因做全面深入分析，可以為科學(xué)評(píng)估HIV/AIDS防控措施的效果、預(yù)測(cè)未來(lái)的發(fā)展趨勢(shì)和采取更有針對(duì)性的防控措施提供參考。

綜上，JPR模型是一種分段線性回歸分析方法，可以采用該模型特有的對(duì)比分析功能對(duì)某些傳染病的分組序列數(shù)據(jù)趨勢(shì)變化特征進(jìn)行量化對(duì)比分析，為深入揭示疾病的流行規(guī)律和開(kāi)展防控效果評(píng)估等提供特定方法，為實(shí)施精準(zhǔn)防控等提供決策參考。當(dāng)然，本研究也存在以下主要不足：一是只選用了一個(gè)案例進(jìn)行了分析，代表性有限；二是沒(méi)有同時(shí)采用其他的趨勢(shì)分析方法進(jìn)行對(duì)比分析。這些有待后續(xù)進(jìn)一步研究。