帶有有序變量的結(jié)構(gòu)方程模型中的模型選擇

李云仙，王學(xué)仁

（1.云南財(cái)經(jīng)大學(xué) 保險(xiǎn)系，昆明 650021；2.云南大學(xué) 統(tǒng)計(jì)系，昆明 650091）

0 引言

有序分類數(shù)據(jù)在社會(huì)學(xué)、心理學(xué)、醫(yī)學(xué)和行為學(xué)中非常常見，例如調(diào)查人們的疼痛不適感可由取值為1～5的有序變量來衡量 （其中1表示非常疼痛，2表示有些疼痛，3表示一般，4表示不疼痛，5表示感覺非常舒服）；人們?cè)谏钪械姆e極性同樣可由類似有序變量來衡量。近年，很多學(xué)者對(duì)此類問題進(jìn)行過研究，其中帶有有序變量的結(jié)構(gòu)方程模型[11]被廣泛用于分析該類問題。在該類模型的應(yīng)用過程中，一個(gè)主要問題就是如何尋找一個(gè)好的模型，使之能夠較好的反映顯變量、協(xié)變量和潛在變量之間的關(guān)系。因此，模型選擇在結(jié)構(gòu)方程模型的應(yīng)用過程中是一個(gè)非常重要的問題。最近，貝葉斯方法在對(duì)結(jié)構(gòu)方程模型的分析中受到了人們的極大關(guān)注[1]。其中，貝葉斯因子[2]是最常用的模型選擇方法之一。由于該類模型結(jié)構(gòu)復(fù)雜，貝葉斯因子的計(jì)算也比較困難[3]，為了解決這個(gè)問題，Gelman and Meng[4]提出了一種比較有效的路徑抽樣方法來計(jì)算一個(gè)概率密度函數(shù)的正態(tài)化因子。之后，這種方法就被廣泛用于計(jì)算結(jié)構(gòu)方程模型的貝葉斯因子[5][6]。然而，當(dāng)兩模型結(jié)構(gòu)相差較大時(shí)，應(yīng)用貝葉斯因子很難將這兩個(gè)模型連接起來，這時(shí)就需要構(gòu)造一些輔助模型。另外，眾所周知，未知參數(shù)先驗(yàn)分布的選取在貝葉斯因子的計(jì)算中很重要?；谶@些局限，本文將應(yīng)用另外一個(gè)貝葉斯統(tǒng)計(jì)量，Lv測(cè)度[7]，對(duì)帶有有序分類變量的結(jié)構(gòu)方程模型進(jìn)行選擇。測(cè)度是一種基于貝葉斯準(zhǔn)則的方法，但對(duì)參數(shù)的先驗(yàn)信息沒有嚴(yán)格要求；同時(shí)，Lv測(cè)度的計(jì)算非常簡(jiǎn)單，可克服貝葉斯因子的一些缺陷。

1 Lv測(cè)度在帶有有序分類變量的結(jié)構(gòu)方程模型中的應(yīng)用

1.1 Lv測(cè)度的定義

其中

v（0≤v≤1）為權(quán)重系數(shù)，當(dāng) v=1時(shí)，該準(zhǔn)則等價(jià)于平方殘差和。從式（1）給出的定義，Lv測(cè)度包括兩部分，第一部分為后驗(yàn)預(yù)測(cè)方差，可以看成是懲罰項(xiàng)；第二部分為預(yù)測(cè)偏差，可以看成是對(duì)模型的擬合優(yōu)度的測(cè)量。因此，具有最小Lv測(cè)度的模型被認(rèn)為是最優(yōu)模型。

1.2 模型定義

假設(shè) yi（i=1，…，n）是一個(gè) p×1 的隨機(jī)向量，同時(shí) yi滿足以下測(cè)量方程：

其中 u（p×1）表示均值向量，ωi（q×1）表示關(guān)于潛在變量的隨機(jī)向量，Λ（p×q）表示載荷矩陣，εi～N（0,Ψε）表示殘差向量，其中 Ψε=diag（ψε1，…，ψεp）是一個(gè)對(duì)角陣。同時(shí)，假設(shè) ωi和εi是相互獨(dú)立的。 潛在變量 ωi可以進(jìn)一步分解為由此可以定義測(cè)量潛在變量間關(guān)系的結(jié)構(gòu)方程為：

其中 ηi（ q1×1）為內(nèi)生潛在變量，ξi（q2×1）為外生潛在變量，q1+q2=q；F（ξi）=（ f1（ξi），…， fm（ξi））T為包含 m 個(gè)實(shí)值可微函數(shù)的向量，其中m≥q。另外，假設(shè)Π0=Iq1-Π是一個(gè)正定矩陣，同時(shí)該矩陣的行列式與 Π 中的元素?zé)o關(guān)；ξi～N（0，Φ）與δi～N（0，Φδ）相互獨(dú)立，其中 Φδ=diag（ψδ1，…，ψδq1）。 如果令 Λω=（Π，Γ），以及那么式（3）可改寫為 ηi=ΛωG（ωi）+δi。 另外，如果令 Λη和 Λξ分別為對(duì)應(yīng)于 ηi和 ξi的Λ的兩個(gè)子矩陣，由式（2）和（3）定義的結(jié)構(gòu)方程模型可以改寫為 yi=u+ΛηΠ0（ΓF（ξi）+δi）+Λξξi+εi。

為了處理有序分類變量，假設(shè) yi=（yo，i，yu，i），其中 yo，i（r×1）為可觀測(cè)連續(xù)隨機(jī)向量，yu，i（s×1）為不可觀測(cè)連續(xù)隨機(jī)向量，r+s=p。其中關(guān)于yu，i的信息可以通過可觀測(cè)有序分類變量zi（s×1）得到，它們之間的關(guān)系定義如下:

其中 zik（k=1，…s）為在{0，1，…，bk}中取值的整數(shù)，αk=（αk，0，…，αk，bk+1）（k=1，…，s）表示閥值。 因此，對(duì)于第 k 個(gè)變量，根據(jù)閥值所定義的類別就有bk+1類。在該模型中，一般假定 αk，0=-∞，αk，bk+1=∞。 根據(jù)很多學(xué)者的研究，為了確保模型的可識(shí)別性，αk，1和 αk，bk，以及Λ和Π中部分元素的值會(huì)預(yù)先給定。為了方便起見，記式（2）到（4）定義的模型為M。

1.3 帶有有序分類變量的結(jié)構(gòu)方程模型中的測(cè)度

為了定義上述出模型的 Lv測(cè)度， 令其中，表示可觀測(cè)連續(xù)變量的樣本矩陣Yu=（yu，1，…，yu，n），為不可觀測(cè)連續(xù)隨機(jī)變量的樣本矩陣同時(shí)令令為有序分類變量的樣本矩陣，其中

根據(jù)這個(gè)定義可得，p（zik，j=1）=p（zik=j-1）=p（αk，j-1＜yu，ik≤αk，j）?pik，j，以及 p（zik，j=0）=1-p（zik，j=1）=1-pik，j。 從上面的定義可知，zik，j～Bernoulli（pik，j），其分布函數(shù)為：

f（zik，j，pik，j）=（pik，j）zik，j（1-pik，j）1-zik，j因此，本模型中的觀測(cè)數(shù)據(jù)即為令和 Zrep分別為 Yobs和 Zobs的預(yù)測(cè)值。 此外，令 Ω =（ω1，…，ωn），Ω1=（η1，…，ηn），Ω2=（ξ1，…，ξn），α 為包括所有閥值的向量，θ為包括模型中所有未知參數(shù)的向量。

在上面介紹的Lv測(cè)度是基于指數(shù)分布族定義的，由于本模型包括兩類變量，連續(xù)變量和有序分類變量。根據(jù)模型的定義，連續(xù)變量的分布屬于指數(shù)分布族，可以直接應(yīng)用式（1）給出的Lv測(cè)度。然而，有序分類變量的分布并不屬于指數(shù)分布族，因此式（1）給出的Lv測(cè)度不能直接應(yīng)用。為了解決這個(gè)問題，本文將采用Chen,Dey and Ibrahim[8]給出的一種方法將有序分類數(shù)據(jù)轉(zhuǎn)化成二元數(shù)據(jù)。引入一個(gè)新的bk+1的向量，其中的元素定義為：

此分布屬于指數(shù)分布族，且該分布的期望和方差分別為

E（zik，j）=pik，j，Var（zik，j）=pik，j（1-pik，j）

根據(jù)式（5），可將Zobs和Zrep中的數(shù)據(jù)全部轉(zhuǎn)換成來自伯努利分布的數(shù)據(jù)，轉(zhuǎn)化后的數(shù)據(jù)分別用Zobs*和Zrep*表示。對(duì)于有序分類數(shù)據(jù)，定義如下二次損失Lv測(cè)度（quadratic loss Lvmeasure）[8]：

其中，ej是一個(gè)bk+1維列向量，除了第j個(gè)元素為1，其他元素全部為0。在上式中

通過簡(jiǎn)單變換，式（6）可寫為

M）+的第 j個(gè)對(duì)角元素是條件期望的第j個(gè)元素。故

根據(jù)模型的定義，在θ和ξ給定的條件下，yu，i具有正態(tài)分布，其均值為協(xié)方差陣為 Λu，η（1ΓF（ξi）T）+Ψεu，其中 uu，Λu，η，Λu，ξ和 Ψεu分別為 u，Λη，Λξ和Ψε的對(duì)應(yīng)于yu，i的子向量或子矩陣。式（8）中最后一個(gè)概率函數(shù)為：

uu，k為 uu的第 k 個(gè)元素，Λu，ηk和 Λu，ξk分別為 Λu，η和 Λu，ξ的第 k 行，ψεuk為 Ψεu的第 k 個(gè)對(duì)角元素，Φ（·）表示標(biāo)準(zhǔn)正態(tài)分布的累計(jì)分布函數(shù)。根據(jù)上面的公式便可以得到式（7）中的條件期望。根據(jù)這個(gè)條件期望，可以進(jìn)一步得到式（7）中的條件方差：

對(duì)于模型中的連續(xù)變量，由式（1）定義的Lv測(cè)度可以直接采用。因此，本文定義的模型M的Lv測(cè)度為：

其中 uo，k為 uo的第 k 個(gè)元素，Λo，ηk和 Λo，ξk分別為o，η和Λo，ξ的第k行。式（9）中第一個(gè)求和項(xiàng)中的條件方差為：

至此，在式（9）中需要的條件期望和方差都得到了，然而在求這些條件期望和條件方差的過程中需要求高階積分，而這些積分是沒有辦法求出來的，因此不能得到Lv測(cè)度的顯式表達(dá)式。本文將采用MCMC的方法來計(jì)算Lv測(cè)度。

根據(jù) Lv測(cè)度的定義，只要得到關(guān)于（θ，α，Ω，Yu）的一個(gè)足夠大的后驗(yàn)樣本，Lv測(cè)度就可以計(jì)算出來了。本文將采用Gibbs抽樣[9]從后驗(yàn)分布中抽取所需要的樣本。 具體為，給定當(dāng)前數(shù)值（θ（g），α（g），Ω（g），）：

在這個(gè)方法收斂后就可以得到關(guān)于未知參數(shù)、潛在變量、閥值、不可觀測(cè)變量的后驗(yàn)樣本{θ（g），Ω（g），α（g），Yu（g）；g=1，…，G}，G是一個(gè)足夠大的數(shù)據(jù)。參數(shù)和潛在變量的后驗(yàn)均值可以作為他們的貝葉斯估計(jì)值，而Lv測(cè)度也可以根據(jù)這個(gè)后驗(yàn)樣本得到。

根據(jù)前面的討論，要完成Gibbs抽樣，計(jì)算測(cè)度的值，就需要求出關(guān)于位置參數(shù)以及潛在變量，閥值等的后驗(yàn)分布。首先，對(duì)于Gibbs抽樣的第一步，需要給出未知參數(shù)的先驗(yàn)分布。本文將采用以下常用的共軛先驗(yàn)分布：

其中IWq2表示q2維數(shù)為逆Wishart分布，在上述先驗(yàn)分布中 ρ0，α0δk，β0δk，α0εj，β0εj，u0，Σ0，Λ0ωk，H0ωk，Λ0k和 H0k是根據(jù)先驗(yàn)信息預(yù)先給定的一些超參數(shù)。根據(jù)這些先驗(yàn)分布，可以很容易的求出關(guān)于它們的后驗(yàn)分布，這些后驗(yàn)分布都是比較熟悉的正態(tài)分布，伽馬分布，Wishart分布，從這些分布中抽取隨機(jī)數(shù)比較簡(jiǎn)單。為了節(jié)省空間，這里將不給出具體的推導(dǎo)過程。對(duì)于閥值α，本文將采用如下無信息先驗(yàn)分布：

其中 C 是一個(gè)常數(shù)。 令 yu（k）=（yu，1k，…，yu，nk），由此可得

根據(jù)前面的討論，在Gibbs抽樣中①中所需要的條件分布是比較常見的分布，可以直接從這些后驗(yàn)分布中抽取未知參數(shù)的隨機(jī)數(shù)。然而，②和③中所需要的后驗(yàn)分布均是不常見的分布，從這些分布中抽取隨機(jī)數(shù)就比較復(fù)雜。為了解決這個(gè)問題，Metropolis-Hastings（MH）方法[10][11]將被應(yīng)用與隨機(jī)抽樣。因此，本文所采用的MCMC算法是Gibbs抽樣和MH方法的一種混合算法。

2 實(shí)例分析

生活質(zhì)量（quality of life,QOL）的度量在臨床實(shí)驗(yàn)以及與健康相關(guān)的分析中非常重要，在這個(gè)問題中，有序分類變量比較常見。針對(duì)這個(gè)問題，Lee[12]將帶有有序變量的結(jié)構(gòu)方程模型應(yīng)用到該問題中進(jìn)行過分析，并應(yīng)用貝葉斯方法估計(jì)結(jié)構(gòu)方程模型的未知參數(shù)。本文將同樣采用帶有有序分類變量的結(jié)構(gòu)方程模型對(duì)QOL數(shù)據(jù)進(jìn)行分析，然而本文的目的在于應(yīng)用測(cè)度進(jìn)行模型選擇。WHOQOL－100量表（世界衛(wèi)生組織生存質(zhì)量測(cè)定量表）是世界衛(wèi)生組織20余個(gè)國(guó)家和地區(qū)共同研制的適用于一般人群的普適性量表。該量表主要用于評(píng)估四個(gè)潛在變量的狀況：身體健康ξ1，心理健康ξ2，社會(huì)關(guān)系ξ3及環(huán)境影響ξ4，以及一個(gè)關(guān)于總體健康狀況（z1）及精神支柱（z2）的綜合因素η。這四個(gè)潛在變量可以由24個(gè)小方面（z3～z26）來反映，它們之間的關(guān)系為：z3～z9反映人的身體健康狀況，z10～z15反映人的心理健康狀況，z16～z18反映社會(huì)關(guān)系，z19～z26反映環(huán)境因素。根據(jù)問題設(shè)計(jì)，這26個(gè)變量均為有序分類變量，其取值均為1～5（其中1表示完全不滿意，2表示有些不滿意，3表示適度，4表示滿意，5表示非常滿意）。整個(gè)樣本數(shù)據(jù)量非常大，這里只是為了說明本文所提出的方法在模型選擇中的應(yīng)用及效用，只分析樣本容量n=338的一組綜合數(shù)據(jù)。本例將比較兩個(gè)結(jié)構(gòu)方程模型：M1和M2，其中M1的定義如下：

M1：y=Λ1ω1+ε，η=γ1ξ1+γ2ξ2+γ3ξ3+γ4ξ4+ε

其中載荷矩陣為：

模型M2定義為：

M2：y= Λ2ω2+ε，η=γ1ξ1+γ2ξ2+γ3ξ3+γ4ξ4+ε

其中載荷矩陣為:

ω2=（η，ξ1，ξ2，ξ3），即將社會(huì)關(guān)系領(lǐng)域及環(huán)境領(lǐng)域合為一個(gè)因素。 在該模型中，假設(shè)，ε～N［0，Ψε2］，ξ＝（ξ1，ξ2，ξ3）～N［0，Φ2］，δ～N［0，］。

在上面兩個(gè)模型中，y=（y1，…，y26）T是隱含的連續(xù)變量，其信息可通過式（4）由觀測(cè)變量 z=（z1，…，z26）T得到。其中，閥值為 α =（α1，…，α26），αk=（αk1，αk2，…，αk5，αk6）（k=1，…26），αk1=-∞，αk6=∞。為了使模型可識(shí)別，需要固定αk2和αk5為一定值，本文采用標(biāo)準(zhǔn)正態(tài)分布來確定αk2和αk5的給定值。具體做法為，分別計(jì)算zk=1概率p1以及zk=5概率p5，則根據(jù)標(biāo)準(zhǔn)正態(tài)分布，可以設(shè)定 αk2=Φ-1（ p1），αk5=Φ-1（ p5），其中 Φ-1表示標(biāo)準(zhǔn)正態(tài)分布累計(jì)分布函數(shù)的反函數(shù)。

為了計(jì)算Lv測(cè)度，需要給定共軛先驗(yàn)分布中先驗(yàn)信息，即超參數(shù)的取值。本例給定如下先驗(yàn)信息：對(duì)于Λ1和Λ2中未知元素的先驗(yàn)分布的均值均設(shè)定為 0.8，對(duì)應(yīng)于（γ1，γ2，γ3，γ4）的先驗(yàn)分布均值設(shè)為（0.6，0.6，0.4，0.4），Φ1和 Φ2的先驗(yàn)分布，逆wishart分布中的超參數(shù)為8I3，其中 Id表示 d 維單位陣，關(guān)于 Ψε1和 Ψε2中的元素，先驗(yàn)分布為伽馬分布，其中的超參數(shù)為 α0εk1=α0εk1=10，β0εk1=β0εk1=8，而對(duì)于 σ1δ和 σ2δ的先驗(yàn)分布，超參數(shù)為 α0δ1=α0δ2=10 α0δ1=α0δ2=10。 本文將采用 R2WinBUGS 來計(jì)算 Lv測(cè)度，其中，所用樣本為Gibbs抽樣過程中刪除前4000次之后的2000次抽樣構(gòu)成。在時(shí)，計(jì)算結(jié)果為：

L0.5=（M1）=7273.01，L0.5=（M2）=7343.82

因此，本文認(rèn)為模型M1優(yōu)于M2。

為了比較不同模型選擇方法的結(jié)果，本文還針對(duì)此例計(jì)算了貝葉斯因子，得logB12=23.51，根據(jù)貝葉斯因子在模型選擇中的準(zhǔn)則,同樣選擇模型。

3 結(jié)論

本文將Lv測(cè)度應(yīng)用到帶有有序分類變量的結(jié)構(gòu)方程模型中進(jìn)行模型選擇，并通過一個(gè)實(shí)例分析說明了該方法的有效性。從實(shí)例分析可以看出，Lv測(cè)度與貝葉斯因子可以得出相同的結(jié)論。然而，Lv測(cè)度的計(jì)算比貝葉斯因子簡(jiǎn)單且節(jié)省時(shí)間，并對(duì)參數(shù)先驗(yàn)分布沒有嚴(yán)格要求。本文所考慮的Lv測(cè)度均是基于v=0.5得到的，當(dāng)然的取值可以在內(nèi)變化，針對(duì)此問題，作者也試過使取其它值，但結(jié)果變化不大，因此本文僅給出v=0.5的結(jié)果。

其次，Lv測(cè)度的應(yīng)用比較廣泛，除本文提出的模型及數(shù)據(jù)類型外，還可以將其應(yīng)用到其他一些諸如帶有無序分類變量的結(jié)構(gòu)方程模型的模型選擇以及其類型的模型。同時(shí)，v的最優(yōu)取值也是一個(gè)比較有意義的研究問題。

[1]R.Schines,H.Hoijtink，A.Boomsma.Bayesian Estimation and Testing of Structural Equation Models[J].Psychometrika，1999，（64）.

[2]R.E.Kass，A.E.Raftery.Bayes Factors[J].Journal of the American Statistical Association，1995，（90）.

[3]T.J.DiCiccio,R.E.Kass,et.al.ComputingBayesFactorsby Combining Simulation and Asymptotic Approximations[J].Journal of the American Statistical Association，1997，（92）.

[4]A.Gelman，X.L.Meng.Simulating Normalizing Constantsfrom Importance Sampling to Bridge Sampling to Path Sampling[J].Sta tistical Science，1998，（13）.

[5]S.Y.Lee，X.Y.Song.Bayesian Model Comparison of Nonlinear Structural Equation Models with Missing Continuous and Ordinal Categorical Data[J].British Journal of Mathematical and Statistical Psychology，2004，（57）.

[6]S.Y.Lee，N.S.Tang.Bayesian Analysis of Structural Equation Models with Mixed Exponential Family and Ordered Categorical Data[J].British Journal of Mathematical and Statistical Psychology，2006，（59）.

[7]J.G.Ibrahim,M.H.Chen，D.Sinha.Criterion-based Methods for Bayesian Model Assessment[J].Statistica Sinica，2001，（11）.

[8]M.H.Chen,D.K.Dey，J.G.Ibrahim.Bayesian Criterion Based Model Assessment for Categorical Data[J].Biometrika，2004，（91）.

[9]S.Geman,D.Geman,Stochastic Relaxation,Gibbs Distribution,and the Bayesian Restoration of Images[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1984，（6）.

[10]N.Metropolis,A.W.Rosenbluth,M.N.Rosenbluth,et.al.Equa tions of State Calculations by Fast Computing Machine[J].Journal of Chemical Physics，1953，（21）.

[11]W.K.Hastings. Monte Carlo Sampling Nethods Using Markov Chains and Their Application[J].Biometrika，1970，（57）.

[12]S.Y.Lee. Structural Equation Modeling:A Bayesian Approach[M].New York:Wiley，2007.