周廣帥 范冰冰 王春霞 游頂云 劉言訓(xùn) 薛付忠 陳 偉 張 濤△

【提 要】 目的 介紹交叉滯后路徑分析原理及其在變量間因果時序關(guān)系研究中的應(yīng)用。方法 交叉滯后路徑分析模型基于交叉滯后面板設(shè)計，估計的路徑系數(shù)具有明確的時間順序關(guān)系，滿足因果推斷中“因在前果在后”的時序性要求。利用健康隨訪數(shù)據(jù)構(gòu)建體重指數(shù)(BMI)與血尿酸(UA)的交叉滯后路徑分析模型，探索BMI和UA的因果時序關(guān)系。結(jié)果 調(diào)整混雜因素后，基線BMI到隨訪時UA的路徑系數(shù)(ρ2=0.060，P<0.001)明顯大于基線UA到隨訪時BMI的路徑系數(shù)(ρ1=-0.009，P=0.056)，且兩系數(shù)間的差異具有統(tǒng)計學(xué)意義(ρ2>ρ1，P<0.001)。在時間順序上BMI增加先于UA升高發(fā)生。結(jié)論 BMI增加可能是高尿酸血癥的原因，交叉滯后路徑分析模型可以有效的識別變量間的因果時序關(guān)系。

實際研究中，通常利用因果圖(casual diagram)的方式直觀地標(biāo)識危險因素之間的相關(guān)或因果關(guān)系，從而清晰地表達(dá)混雜效應(yīng)、中介效應(yīng)等多種因果推斷的關(guān)鍵概念[1]。如圖1所示的因果圖模型，X為解釋變量，Y為反應(yīng)變量，Z為另一個解釋變量。虛線雙箭頭表示存在關(guān)聯(lián)，實線箭頭表示有因果關(guān)系[2]。若X→Z則為中介效應(yīng)模型，反之為混雜效應(yīng)。因此，確定X和Z之間的因果方向是區(qū)分混雜和中介的關(guān)鍵。

圖1 混雜和中介的因果圖模型示意圖

交叉滯后路徑分析(cross-lagged path analysis)是利用縱向隊列中兩個及以上的面板數(shù)據(jù)，探索兩個變量間因果時序關(guān)系(temporal sequence)的一種統(tǒng)計分析方法[3]。貝葉斯網(wǎng)絡(luò)可采用有向無環(huán)圖表示因果網(wǎng)絡(luò)，被廣泛應(yīng)用于因果機(jī)制的發(fā)現(xiàn)和因果推斷中，但其通常是基于橫斷面數(shù)據(jù)。交叉滯后路徑分析的優(yōu)勢是結(jié)合了交叉滯后面板設(shè)計(cross-lagged panel design)，估計的路徑系數(shù)具有明確的時間順序關(guān)系，符合“因在前果在后”的因果推斷原則，常用于變量間因果時序關(guān)系的研究[4-7]。

本文將介紹交叉滯后路徑分析的基本原理和分析步驟，并通過體重指數(shù)(BMI)與血尿酸(UA)因果時序關(guān)系的研究實例介紹該方法的實際應(yīng)用。

交叉滯后路徑分析

交叉滯后路徑分析在交叉滯后面板設(shè)計基礎(chǔ)上，研究兩變量間的因果時序關(guān)系，在數(shù)據(jù)上要求兩個及以上的面板數(shù)據(jù)，其模型示意圖如圖2所示。假設(shè)任意兩變量分別為X和Y，其基線水平分別為X1和Y1，隨訪水平分別為Xt和Yt。其中β1和β2為自相關(guān)系數(shù)(autoregressive coefficient)，表示隨訪測量(Xt或Yt)與基線測量(X1或Y1)的自相關(guān)關(guān)系；r1為同步相關(guān)系數(shù)(synchronous coefficient)，是基線X1與Y1的相關(guān)系數(shù)；ρ1和ρ2為cross-lagged路徑系數(shù)(cross-lagged path coefficient)。需注意的是，交叉滯后路徑分析模型中并不包含Xt→Yt和Yt→Xt的因果路徑，即Xt和Yt之間不存在同步效應(yīng)。否則，模型中需要估計的參數(shù)個數(shù)將超過已有的方差和協(xié)方差數(shù)量，模型將無法識別[8-9]。

圖2 交叉滯后路徑分析模型示意圖

在路徑分析框架下，交叉滯后路徑分析模型的形式如下：

Xt=α1+β1X1+ρ1Y1+e1

Yt=α2+β2Y1+ρ2X1+e2

其中α1和α2為模型截距，e1和e2為模型殘差。

路徑系數(shù)ρ是交叉滯后路徑分析進(jìn)行因果時序關(guān)系推斷的依據(jù)，其含義為：在控制了因變量基線狀態(tài)(如Y1)后，自變量(如X1)基線測量結(jié)果對因變量(如Yt)的影響，即Yt=α2+β2Y1+ρ2X1+e2。該模型經(jīng)轉(zhuǎn)換可表示為Y的動態(tài)變化同基線X1和Y1的回歸形式，即ΔY=(Yt-Y1)=α2+(β2-1)Y1+ρ2X1+e2。計算的路徑系數(shù)ρ在時間上具有明確的先后順序(X1→Yt，Y1→Xt)。

通過路徑系數(shù)ρ1和ρ2的比較可以確定變量間的因果時序關(guān)系，包含以下四種情況：若ρ1=0，ρ2=0，則X，Y間無因果時序關(guān)系；若ρ1≠0，ρ2=0，且兩系數(shù)差異顯著，則兩變量為Y→X的單向因果時序關(guān)系；若ρ1=0，ρ2≠0，且兩系數(shù)差異顯著，則兩變量為X→Y的單向因果時序關(guān)系；若ρ1≠0，ρ2≠0，則X，Y為相互調(diào)控的雙向因果時序關(guān)系，此時若ρ1>ρ2且兩者差異顯著，可進(jìn)一步確定主要因果時序效應(yīng)為Y→X，反之則為X→Y。

研究實例

肥胖和高尿酸血癥是心腦血管疾病的危險因素，大量流行病學(xué)和臨床研究發(fā)現(xiàn)肥胖和高尿酸血癥存在顯著的相關(guān)性，然而其因果時序關(guān)系仍然存在爭議。通常認(rèn)為體重指數(shù)(BMI)增加會導(dǎo)致高尿酸血癥，但也有研究發(fā)現(xiàn)血尿酸(UA)水平高的人群更易患有肥胖[10-12]。此外，基礎(chǔ)研究表明尿酸可引起細(xì)胞內(nèi)脂肪累積，導(dǎo)致體重增加[13]。本文將通過交叉滯后路徑分析探討B(tài)MI和UA之間的因果時序關(guān)系。

本研究基于山東某三甲醫(yī)院2004到2015年的健康隨訪數(shù)據(jù)，選取基線年齡≥20歲，體檢記錄中包含兩次及以上BMI和UA測量的體檢參與者。由于數(shù)據(jù)中不包含高尿酸血癥的相關(guān)治療信息，為防止服藥對尿酸水平的影響，排除了患有高尿酸血癥的人群。本研究最終共納入15348名對象，選取其基線和最后一次隨訪數(shù)據(jù)進(jìn)行交叉滯后路徑分析。研究對象的基線與隨訪特征如表1所示，其中男性9578人，女性5770人，基線平均年齡為39.68歲，平均隨訪時間為3.36年。

表1 研究隊列特征

交叉滯后路徑分析的建模步驟

交叉滯后路徑分析作為路徑分析的一種，可使用結(jié)構(gòu)方程模型的分析軟件建立模型，如：LISREL、AMOS、SAS中CALIS過程步、R中sem包、lavaan包等。

1.數(shù)據(jù)預(yù)處理

交叉滯后路徑分析前，為控制潛在的混雜效應(yīng)和消除量綱影響，需對研究變量進(jìn)行回歸殘差分析，具體步驟為：(1)構(gòu)建研究變量與混雜因素(如：年齡、性別、吸煙、飲酒、隨訪時間等)的多元回歸方程；(2)對上述方程取殘差，并進(jìn)行Z轉(zhuǎn)化的標(biāo)準(zhǔn)化處理。

2.計算協(xié)方差矩陣

結(jié)構(gòu)方程模型又稱協(xié)方差結(jié)構(gòu)模型，其分析軟件多以協(xié)方差矩陣為基礎(chǔ)開始模型擬合，如LISREL，sem包等。此外一些結(jié)構(gòu)方程模型分析軟件也可以直接使用原始數(shù)據(jù)進(jìn)行模型分析，如：AMOS、CALIS過程步、lavaan包等。交叉滯后路徑分析模型在預(yù)處理時已將研究變量標(biāo)準(zhǔn)化，其協(xié)方差矩陣與相關(guān)系數(shù)矩陣相同。

3.模型構(gòu)建

(1)根據(jù)現(xiàn)有研究確定各時點BMI與UA間的路徑方向，建立如圖3所示的交叉滯后路徑分析模型，其結(jié)構(gòu)方程組如下：

BMIt=β1BMI1+ρ1UA1+e1

UAt=β2UA1+ρ2BMI1+e2

上述方程組中基線BMI1、基線UA1、隨訪BMIt、隨訪UAt的協(xié)方差矩陣為[14]：

∑=

其中r1為基線BMI、基線UA的同步相關(guān)系數(shù)；re為e1、e2的相關(guān)系數(shù)。

(2)通過結(jié)構(gòu)方程模型分析軟件對上述模型進(jìn)行參數(shù)估計。本研究采用lavaan包直接對原始數(shù)據(jù)進(jìn)行模型分析。具體代碼如下，#號后文字是對代碼的注釋：

install.packages(“l(fā)avaan”)#安裝lavaan包

library(lavaan)#載入lavaan包

model1<-"BMIt～a1*BMI1+b1*UA1，

UAt～a2*UA1+b2*BMI1，

BMI1～～r1*UA1，

BMIt～～0*UAt"#輸入模型公式

result<-sem(model=model1，data=usedata)#模型求解

summary(result)#查看模型結(jié)果

fitMeasures(result，c(“rmr”，“cfi”))#查看模型擬合情況

其中r1對應(yīng)同步相關(guān)系數(shù)r1，a1、a2分別對應(yīng)自相關(guān)系數(shù)β1和β2，b1、b2分別對應(yīng)路徑系數(shù)ρ1和ρ2；使用sem函數(shù)進(jìn)行求解，model為需要輸入的模型結(jié)構(gòu)，data為分析所用的數(shù)據(jù)，此處usedata已經(jīng)過預(yù)處理控制了混雜因素的影響；使用summary查看上述分析結(jié)果，使用fitMeasures查看模型擬合情況。

Lavaan包中默認(rèn)使用極大似然法(maximum likelihood，ML)進(jìn)行參數(shù)估計。其原理為通過迭代計算使模型隱含的協(xié)方差矩陣即“再生矩陣”與樣本協(xié)方差矩陣間“距離”最小，得到模型最優(yōu)解。ML屬于完全信息方法，具有無偏、一致、尺度不變和漸進(jìn)有效性的優(yōu)點，是當(dāng)前結(jié)構(gòu)方程模型最常用的參數(shù)估計方法，其假設(shè)為變量服從多元正態(tài)分布[14]。

(3)判斷模型參數(shù)估計是否合理，并通過擬合指數(shù)評價模型擬合情況。

常用的模型擬合指數(shù)有CFI(comparative fit index)和RMR(root mean-square residual)。其中CFI反映當(dāng)前假定模型與變量間相互獨立的虛模型(最不理想模型)的差異程度，具有不易受樣本量影響，小樣本時表現(xiàn)良好的優(yōu)點，其值越大擬合越好；RMR則是基于假定模型整體殘差的模型擬合指數(shù)，具有對誤設(shè)模型較為敏感的優(yōu)點，其值越小擬合越好。當(dāng)CFI>0.90且RMR<0.05時表示模型擬合良好[4，15]。

(4)通過Fisher Z′ test對路徑系數(shù)進(jìn)行比較，確定主要因果時序效應(yīng)。

實例分析結(jié)果

調(diào)整年齡、性別、吸煙、飲酒、隨訪時間后，BMI和UA在基線與隨訪時的相關(guān)系數(shù)和散點圖如圖4所示，其中基線UA與隨訪UA的相關(guān)系數(shù)為0.615，基線BMI與基線UA、隨訪BMI的相關(guān)系數(shù)分別為0.246、0.831，圖中所有的相關(guān)系數(shù)均具有統(tǒng)計學(xué)意義。

圖3 BMI與UA的交叉滯后路徑分析模型

BMI與UA的交叉滯后路徑分析模型結(jié)果如圖3所示，調(diào)整年齡、性別、吸煙、飲酒、隨訪時間后，由基線BMI到隨訪UA的路徑系數(shù)ρ2(ρ2=0.060，P<0.001)明顯大于由基線UA到隨訪BMI的路徑系數(shù)ρ1(ρ1=-0.009，P=0.056)，且兩系數(shù)間差異具有統(tǒng)計學(xué)意義(ρ2>ρ1，P<0.001)。BMI和UA兩次測量間的自相關(guān)系數(shù)分別為0.834和0.600。模型擬合指數(shù)CFI和RMR分別為0.979和0.027，提示當(dāng)前模型擬合良好。

圖4 基線與隨訪的BMI、UA相關(guān)系數(shù)和散點圖

如圖5所示，調(diào)整年齡、性別、吸煙、飲酒后，隨著基線BMI四分位數(shù)的增加UA的年變化率逐漸升高，且增長趨勢具有統(tǒng)計學(xué)意義；而隨著基線UA四分位數(shù)的增加BMI的年變化率基本保持不變，增長趨勢沒有統(tǒng)計學(xué)意義。圖5結(jié)果與圖3中交叉滯后路徑分析結(jié)果一致。

圖5 不同基線四分位數(shù)下UA、BMI的年變化率

此外，本文選取隨訪隊列中體檢三次及以上的4797名參與者的隨訪信息，構(gòu)建了三斷面交叉滯后路徑分析模型，進(jìn)行了敏感性分析。模型結(jié)果如圖6所示，調(diào)整混雜因素后，在不同斷面間均存在的單向因果時序關(guān)系(ρ2=0.046，P<0.001；ρ4=0.061，P<0.001)，與圖3兩斷面模型結(jié)果一致。CFI=0.931，RMR=0.052提示當(dāng)前模型擬合良好。綜上，可以認(rèn)為BMI與UA間因果時序方向為BMI→UA，即在時間順序上BMI增加先于UA水平升高發(fā)生，BMI增加可能是導(dǎo)致高尿酸血癥的原因。

討 論

變量間的時序關(guān)系是流行病學(xué)中進(jìn)行因果推斷的重要原則之一，也是因果圖模型中確定節(jié)點間箭頭方向的基礎(chǔ)。交叉滯后路徑分析是利用縱向隊列中兩個及以上面板數(shù)據(jù)，探索兩變量間因果時序關(guān)系的方法，其估計的路徑系數(shù)具有明確的時間順序關(guān)系，符合“因在前果在后”的因果推斷原則，為進(jìn)一步的因果關(guān)系探索提供了依據(jù)。

交叉滯后路徑分析需滿足以下假設(shè)：(1)變量服從多元正態(tài)分布，當(dāng)違背正態(tài)性假定時，使用ML法的模型參數(shù)標(biāo)準(zhǔn)誤將被低估，增大I型錯誤風(fēng)險，此時可通過變量轉(zhuǎn)換，或者通過S-B調(diào)整的ML法、漸進(jìn)自由分布法、自助抽樣法等進(jìn)行模型擬合[14]；(2)人群內(nèi)部一致性，交叉滯后路徑分析關(guān)注總體層面的平均因果效應(yīng)，未充分考慮個體間差異，當(dāng)研究變量個體變異較大時可能影響模型結(jié)果[16]；(3)隨訪間隔接近真實因果效應(yīng)時間，變量間的相互作用常需一定發(fā)生時間，如果隨訪間隔太短或太長，因果效應(yīng)可能尚未發(fā)生或已經(jīng)消失，影響時序關(guān)系判斷，模型選取的隨訪間隔應(yīng)盡量接近真實的因果效應(yīng)時間，但實際應(yīng)用中受限于現(xiàn)有數(shù)據(jù)，該假設(shè)往往容易被違背[17]。此外需注意的是交叉滯后路徑分析仍借助常規(guī)回歸方法控制混雜效應(yīng)，模型結(jié)果可能受到未知混雜影響。另一方面，交叉滯后路徑分析屬于驗證性分析，模型假定的因果方向均應(yīng)來自現(xiàn)有人群或理論研究。

兩變量兩斷面交叉滯后路徑分析是該模型最基本的形式，當(dāng)包含三個或者更多面板時，可構(gòu)建包含多個隨訪間隔的交叉滯后路徑分析模型，通過比較各隨訪期內(nèi)的路徑系數(shù)，探索不同時點下變量間因果時序方向和強(qiáng)度的變化規(guī)律[18]。多斷面模型中，每個面板的隨訪人群均應(yīng)處于相同年齡段或生命歷程階段(如：兒童發(fā)育階段、女性絕經(jīng)狀態(tài)等)。其中三斷面模型要求三個隨訪斷面均處于相同生命歷程階段(如：均為成年期)，常用于兩斷面模型的敏感性分析，以檢驗隨訪間隔改變對模型結(jié)果的影響。本文實例分析部分BMI與UA的三斷面交叉滯后路徑分析模型與兩斷面模型結(jié)果一致，均提示BMI增加是導(dǎo)致高尿酸血癥的原因。多斷面模型的模型構(gòu)建和求解過程與兩斷面模型類似，鑒于文章篇幅不再展開。此外，隨著模型發(fā)展交叉滯后路徑分析已衍生多種模型形式，關(guān)注較多的有RI-CLPM(random-intercept cross-lagged panel analysis model)、ALT-SD(autoregressive latent trajectory model with structured residuals)等，前者通過在模型中引入截距因子(潛變量)，為控制個體間差異對變量間真實時序關(guān)系的影響提供了可能[16]；后者則將交叉滯后路徑分析與潛變量增長曲線模型結(jié)合，使模型在估計變量間因果時序關(guān)系的同時也可獲得變量的總體變化軌跡，擴(kuò)展了交叉滯后路徑分析模型的應(yīng)用范圍[19]。