基于證候要素的二階證實(shí)性因子分析在R語言lavaan包中的實(shí)現(xiàn)*

黃 灝 薛芳靜 許碧云 申春悌 陳啟光 陳炳為△

基于證候要素的二階證實(shí)性因子分析在R語言lavaan包中的實(shí)現(xiàn)*

黃 灝1薛芳靜1許碧云2申春悌3陳啟光1陳炳為1△

目的探索高血壓陰虛夾濕證的證候要素，并介紹R語言lavaan包在結(jié)構(gòu)方程模型中的應(yīng)用。方法介紹lavaan包的語法，然后利用lavaan包中的二階證實(shí)性因子分析對1280例高血壓的陰虛夾濕證進(jìn)行分析，以得到高血壓陰虛夾濕證的證候要素。結(jié)果陰虛夾濕證包含3個(gè)證候要素，病位主要在肝與腎、病性主要為虛與濕。結(jié)論利用R語言lavaan包實(shí)現(xiàn)結(jié)構(gòu)方程模型是十分方便的。二階證實(shí)性因子分析方法可用于中醫(yī)證候分類及證候要素提取和命名。

lavaan包 二階證實(shí)性因子分析 證候要素

辨證論治是中醫(yī)學(xué)的精髓，證候是指內(nèi)外環(huán)境變化所導(dǎo)致的機(jī)體病理生理整體反應(yīng)狀態(tài)，是疾病發(fā)生和演變過程中某一階段病理本質(zhì)的反映，它以一組相關(guān)的癥狀和體征為依據(jù)，從不同程度揭示出患者當(dāng)前的病機(jī)綜合而成。證候要素是指組成“證候”本質(zhì)的基本構(gòu)成單位，是指組成證候的主要元素。本文利用二階證實(shí)性因子分析對中醫(yī)證候中的陰虛夾濕證進(jìn)行證候要素分解［1］。

二階證實(shí)性因子分析屬于結(jié)構(gòu)方程模型（structural equation model，SEM）的一種［2］。目前市場有很多軟件可用于擬合SEM模型，如R、LISREL、Mplus、EQS、Stata和AMOS等。其中R是一款免費(fèi)的、開源的統(tǒng)計(jì)編程語言和計(jì)算環(huán)境［3］。R中有多個(gè)SEM的包，包括：lavaan、sem、lava和OpenMx［4］。目前使用較多的有l(wèi)avaan包和sem包，本文以中醫(yī)證候要素的提取為實(shí)例介紹R語言lavaan包中結(jié)構(gòu)方程模型的實(shí)現(xiàn)。借助R語言中l(wèi)avaan包中的二階證實(shí)性因子分析對中醫(yī)證候中的陰虛夾濕證進(jìn)行證候要素的分解。

結(jié)構(gòu)方程模型及l(fā)avaan包簡介

1．結(jié)構(gòu)方程模型簡介

在SEM中包括了兩種子模型，第一個(gè)為測量模型（measurementmodel），它是確定觀察變量與潛在變量之間關(guān)系；另一為結(jié)構(gòu)模型（structuralmodel），是表示潛變量與潛變量之間關(guān)系的模型。它用公式表達(dá)如下：

測量模型：X＝ΛXξ＋δ與Y＝Λyη＋ε

結(jié)構(gòu)模型：η＝Bη＋Γξ＋ζ

其中ΛX是外生顯在變量（可觀測變量）X在外生潛在變量ξ上的因子載荷矩陣，ΛY是內(nèi)生顯在變量Y在內(nèi)生潛在變量η上因子載荷矩陣。δ、ε、ζ分別表示誤差［1］。下面以4個(gè)內(nèi)生顯在變量及4個(gè)外生顯在變量，2個(gè)內(nèi)生潛在變量及2個(gè)外生潛在變量來說明，其通徑關(guān)系如圖1。

圖1 一個(gè)假設(shè)的結(jié)構(gòu)方程的路徑圖

2．lavaan包簡介

對于R語言的lavaan包中需先對模型進(jìn)行定義，在模型定義中主要用到一些符號(hào)與關(guān)鍵字，將其列出如表1。

在R語言中，首先需要安裝lavaan包，安裝程序如下：

install．packages（“l(fā)avaan”）＃安裝lavaan包

library（“l(fā)avaan”）＃調(diào)用lavaan包

lavaan包中結(jié)構(gòu)方程模型的分析主要從以下幾步考慮：

（1）模型指定：以圖1的模型來說明，生成模型代碼（mymodel）如下：

mymodel＜－′ks1＝～x1＋x2＃ks1由變量x1與x2度量

ks2＝～x3＋x4＃ks2由變量x3與x4度量

et1＝～y1＋y2＃et1由變量y1與y2度量

et2＝～y3＋y4＃et2由變量y3與y4度量

et1～ks1＋ks2＃ks1與ks2預(yù)測et1

et2～ks1＋ks2＋et1＃ks1、ks2與et1預(yù)測et2

ks1～～ks2′＃ks1與ks2相關(guān)

表1 R語言中l(wèi)avaan包中主要命令

（2）模型運(yùn)行：當(dāng)完成結(jié)構(gòu)方程模型設(shè)置后，可運(yùn)行如下程序。

fit＜－sem（model＝，data＝，meanstructure＝，sample．nobs＝，sample．mean＝，sample．cov＝，group＝，bootstrap＝，estimator＝）

當(dāng)模型為結(jié)構(gòu)方程模型時(shí)，采用sem命令；當(dāng)模型為證實(shí)性因子分析時(shí)，采用cfa或sem命令；當(dāng)模型為潛在生長曲線模型時(shí)，采用grow th或sem命令。model用于指定模型名稱，如前面自定義的模型名稱為mymodel。對于結(jié)構(gòu)方程模型一般考慮協(xié)方差矩陣來估計(jì)參數(shù)，但有些模型需要指定常數(shù)項(xiàng)，如潛在生長曲線模型，這時(shí)除了協(xié)方差矩陣外，還需要均數(shù)向量，這時(shí)模型的meanstructure指定TRUE，默認(rèn)為不考慮均數(shù)向量。對于分析的數(shù)據(jù)有兩種格式，第一種是原始數(shù)據(jù)，則可用data選項(xiàng)，其中數(shù)據(jù)集為R語言的frame格式。如果數(shù)據(jù)形式為協(xié)方差矩陣時(shí)，則可用sample．nobs＝樣本例數(shù)，sample．cov＝樣本協(xié)方差矩陣的名稱（協(xié)方差的格式為下三角陣形式），sample．mean＝樣本均數(shù)向量的名稱。如果模型中需要考慮多個(gè)組的結(jié)構(gòu)方程模型，則選項(xiàng)group＝組別變量名；如果需要進(jìn)行bootstrap抽樣，則bootstrap＝抽取的次數(shù)；estimator選項(xiàng)是指定模型估計(jì)的方法，對應(yīng)的選項(xiàng)有ML、GLS、WLS、ULS、DWLS等，即采用最大似然估計(jì)、廣義最少二乘估計(jì)、加權(quán)最二乘估計(jì)、不加權(quán)最少二乘估計(jì)及對角線加權(quán)最少二乘估計(jì)等。R語言統(tǒng)計(jì)分析將結(jié)果都存在一個(gè)fit的對象（object）里面。

（3）結(jié)果輸出：

R語言在分析執(zhí)行結(jié)束后并不顯示任何結(jié)果，但可以有選擇地從結(jié)果中提取感興趣的部分［5］?？梢酝ㄟ^以下的命令來提取結(jié)果。

①summary（fit，standardized＝，fit．measures＝，rsquare＝）

summary對模型擬合情況進(jìn)行描述，給出模型擬合的卡方值、自由度和P值，同時(shí)給出統(tǒng)計(jì)量的方差。Standardized選項(xiàng)值為TRUE時(shí)輸出標(biāo)準(zhǔn)化回歸系數(shù)，否則默認(rèn)不輸出。fit．measures＝TRUE時(shí)輸出模型擬合的CFI、TLI、AIC、BIC、RMSEA及其95%可信區(qū)間、SRMR等擬合度的指數(shù)，默認(rèn)時(shí)只輸出卡方、自由度及P值。rsquare＝TRUE給出每個(gè)變量（可觀測變量及潛在變量）的復(fù)相關(guān)系數(shù)。

②modindices（fit，sort．＝，minimum．value＝，maximum．number＝）

modindices給出模型的修正指數(shù)MI值，為模型修正提供依據(jù)。sort．＝TRUE，按MI值從大到小排序，默認(rèn)不進(jìn)行排序。minimum．value＝m，只保留MI≥m的修正結(jié)果，maximum．number＝n，列出前n個(gè)較大MI值。

③fitMeasures（fit，“擬合指數(shù)名稱”）

fitMeasures函數(shù)計(jì)算了各種模型擬合的指數(shù)。默認(rèn)給出所有的衡量指標(biāo)，包括：CFI、TLI、GFI、AGFI、SRMR、RMSEA等。若只需要單個(gè)指標(biāo)，只需要單個(gè)指標(biāo)加入函數(shù)中即可，如fitMeasures（fit，“cfi”）。

④standardizedSolution（fit，type＝）

standardizedSolution描述潛變量模型的標(biāo)準(zhǔn)化系數(shù)。type默認(rèn)為“std．a(chǎn)ll”，即對顯變量和潛變量的分別進(jìn)行標(biāo)準(zhǔn)化，若選項(xiàng)為“std．lv”，則只對模型中潛變量的進(jìn)行標(biāo)準(zhǔn)化。因此，type大部分情況采用“std．a(chǎn)ll”的選項(xiàng)。

實(shí)例分析

資料來自2006年7月至2008年12月在常州、南京、沈陽和珠海四個(gè)地區(qū)五個(gè)三級(jí)甲等中醫(yī)院收集到的高血壓病例共計(jì)1280例。證候?yàn)殛幪搳A濕的中醫(yī)四診信息有：目澀（x1）、耳鳴（x2）、目癢（x3）、目脹（x4）、耳聾（x5）、健忘（x6）、視物模糊（x7）、迎風(fēng)流淚（x8）、足痛（x9）、尿后余瀝（x10）、小便不暢（x11）、便秘（x12）、小便清長（x13）、身重（x14）、口粘膩（x15）、夜間多尿（x16）、白苔（x17）、口咽干燥（x18）、浮腫（x19）、舌紫（x20）、四肢麻木（x21）指標(biāo)，數(shù)據(jù)集稱為GHKX。

對于結(jié)構(gòu)方程模型初始模型的建立有兩種方法，第一種方法是采用探索性因子分析來探索包含幾個(gè)因子；第二種方法是根據(jù)專業(yè)的知識(shí)構(gòu)造初始模型。先采用因子分析得到陰虛夾濕包含三因子，并根據(jù)因子載荷系數(shù)構(gòu)建二階證實(shí)性因子分析初始模型（mymodel）：

構(gòu)建初始模型后，采用cfa或sem命令對數(shù)據(jù)進(jìn)行擬合：

fit＜-cfa（mymodel，data＝GHKX）

summary（fit，standardized＝TRUE，fit．measures＝TRUE）

modindices（fit，sort．＝TRUE，minimum．value＝10）

采用summary對結(jié)果進(jìn)行輸出，得到部分?jǐn)M合度指數(shù)、載荷系數(shù)、方差及協(xié)方差等。利用modindices列出模型的修正指數(shù)MI值，根據(jù)每一步的MI值大小，對模型進(jìn)行調(diào)整。模型修正中x10與x11間的相關(guān)系數(shù)MI最大，其值為78．114，它是模型修正值。因此，在初始模型中加入“x10～～x11”重新擬合模型。如此循環(huán)反復(fù)，直到所有的MI值在5～10之間則停止修正，從而得到最終的模型。最終模型利用summary、fit-Measures與standardizedSolution得到擬合度、載荷系數(shù)及標(biāo)準(zhǔn)化系數(shù)。

表2 模型整體擬合度評價(jià)指標(biāo)

模型擬合度中CFI、TLI、GFI、AGFI越大說明模型擬合越好［6］。Hu和Bentler（1998）推薦SRMR小于0．08認(rèn)為模型可以接受［7］。RMSEA（近似誤差均方根）＜0．05擬合好；0．08～0．10擬合一般；＞0．10擬合不好。表2中x/df＝1．596＜3，CFI、TLI、GFI、AGFI均大于＞0．90，SRMR與RMSEA均小于0．05，可見最終模型的擬合較好。

表3 二階證實(shí)性因子分析最終模型的標(biāo)準(zhǔn)化系數(shù)

結(jié)果分析可見，陰虛夾濕證包含3個(gè)證候要素，證候要素1包含6個(gè)指標(biāo)，其中5個(gè)指標(biāo)（目澀、目癢、健忘、耳鳴、迎風(fēng)流淚）標(biāo)準(zhǔn)化系數(shù)高于0．4，證候要素1可以解釋為病位在肝，病性為氣虛；證候要素2包含3個(gè)指標(biāo)，其中2個(gè)指標(biāo)（尿后余瀝、小便不暢）標(biāo)準(zhǔn)化系數(shù)高于0．4，證候要素2可以解釋為病位在腎，病性為陰虛；證候要素3包含5個(gè)指標(biāo)，標(biāo)準(zhǔn)化系數(shù)都高于0．4，證候要素3可以解釋為病性為濕。由此可見，陰虛夾濕證與篩選出的14個(gè)指標(biāo)間有高度相關(guān)性，成功實(shí)現(xiàn)了中醫(yī)的證候要素的提取。

討 論

證候要素的提取便于中醫(yī)證候的命名，及中醫(yī)上的診療，因此將證候分解成中醫(yī)的證候要素有重要臨床意義。目前應(yīng)用較多的實(shí)現(xiàn)潛在變量模型的軟件有LISREL、Mplus及AMOS，以上三個(gè)軟件均為商業(yè)版，R語言是一款免費(fèi)的統(tǒng)計(jì)軟件包，在醫(yī)學(xué)研究中目前應(yīng)用越來廣泛，每個(gè)軟件都有各自的優(yōu)缺點(diǎn)。AMOS可以通過通徑圖來構(gòu)建模型，也可通過程序，而R語言、LISREL、M plus是利用編程得到的。然而在模型修正時(shí)，選擇基于編程語言的軟件如R和LISREL比點(diǎn)擊式操作的軟件更具有系統(tǒng)性和便利性。

［1］陳炳為，陳啟光，許碧云．潛在變量模型及其在中醫(yī)證候中的應(yīng)用概述．中國衛(wèi)生統(tǒng)計(jì)，2009，26（5）：535-538．

［2］Oberski DL．Lavaan survey：an R package for complex survey analysis of structural equation models．Journal of Statistical Software，2014，57（1）：1-27．

［3］Fox J．Teacher′scorner：Structural equation modeling with the sempackage in R．Structural equation modeling：A Multidisciplinary Journal，2006，13（3）：465-486．

［4］Song YE，Stein CM，Morris NJ．Strum：an R package for structural modeling of latent variables for general pedigrees．BMCgenetics，2015，16（1）：1-2．

［5］湯銀才．R語言與統(tǒng)計(jì)分析．高等教育出版社，2008．

［6］秦浩，陳景武．結(jié)構(gòu)方程模型原理及其應(yīng)用注意事項(xiàng)．中國衛(wèi)生統(tǒng)計(jì)，2009，23（4）：367-369．

［7］侯杰泰，溫忠麟，成子娟．結(jié)構(gòu)方程模型及其應(yīng)用．北京：教育科學(xué)出版社，2004．

（責(zé)任編輯：鄧 妍）

國家自然科學(xué)基金（81273190）

1．東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)系（210009）

2．南京大學(xué)醫(yī)學(xué)院附屬鼓樓醫(yī)院

3．南京中醫(yī)藥大學(xué)附屬常州中醫(yī)院

△通信作者：陳炳為，E-mail：drchenbw＠126．com