李云仙，王學(xué)仁

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

0 引言

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

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

1.1 Lv測度的定義

其中

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

1.2 模型定義

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

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

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

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

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

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

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

根據(jù)這個定義可得，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因此，本模型中的觀測數(shù)據(jù)即為令和 Zrep分別為 Yobs和 Zobs的預(yù)測值。 此外，令 Ω =（ω1，…，ωn），Ω1=（η1，…，ηn），Ω2=（ξ1，…，ξn），α 為包括所有閥值的向量，θ為包括模型中所有未知參數(shù)的向量。

在上面介紹的Lv測度是基于指數(shù)分布族定義的，由于本模型包括兩類變量，連續(xù)變量和有序分類變量。根據(jù)模型的定義，連續(xù)變量的分布屬于指數(shù)分布族，可以直接應(yīng)用式（1）給出的Lv測度。然而，有序分類變量的分布并不屬于指數(shù)分布族，因此式（1）給出的Lv測度不能直接應(yīng)用。為了解決這個問題，本文將采用Chen,Dey and Ibrahim[8]給出的一種方法將有序分類數(shù)據(jù)轉(zhuǎn)化成二元數(shù)據(jù)。引入一個新的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*表示。對于有序分類數(shù)據(jù)，定義如下二次損失Lv測度（quadratic loss Lvmeasure）[8]：

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

通過簡單變換，式（6）可寫為

M）+的第 j個對角元素是條件期望的第j個元素。故

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

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

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

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

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

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

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

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

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

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

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

2 實例分析

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

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

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

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

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

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

3 結(jié)論

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

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

[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.