胡良平

（1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專(zhuān)業(yè)委員會(huì)，北京 100029 *通信作者：胡良平，E-mail：lphu812＠sina.com）

1 概 述

1.1 有限混合模型

有限混合模型［1］見(jiàn)式（1）：

在式（1）中，假定y是一個(gè)可觀測(cè)的隨機(jī)變量，其分布取決于一個(gè)不可觀測(cè)的隱變量S，它有k種可能的狀態(tài)；k的取值是未知的，但至少是有限的。

令πj表示S取值j的概率，在S=j的條件下，假定反應(yīng)變量 Y的分布是 fj（y；αj，βj|S=j）（注：它是一個(gè)概率密度函數(shù)）。

1.2 有限混合模型回歸分析應(yīng)用的場(chǎng)合

根據(jù)研究目的確定了一個(gè)具有同質(zhì)性的研究總體，若從該總體中隨機(jī)抽取樣本含量為n的個(gè)體，從每個(gè)個(gè)體身上測(cè)量某計(jì)量指標(biāo)的數(shù)值，依此法收集到的n個(gè)計(jì)量數(shù)據(jù)被稱(chēng)為“單組設(shè)計(jì)一元計(jì)量資料”。描述這組計(jì)量數(shù)據(jù)的方法有如下幾種：第一，編制頻數(shù)分布表；第二，繪制直方圖；第三，擬合頻數(shù)分布曲線。

當(dāng)直方圖顯示只有一個(gè)高峰時(shí)，就稱(chēng)其為“單峰分布”。此時(shí)，可以根據(jù)高峰所處的位置，分為“正偏態(tài)分布（高峰位于左側(cè)）”“對(duì)稱(chēng)分布（高峰居中，其特例為正態(tài)分布）”和“負(fù)偏態(tài)分布（高峰位于右側(cè)）”。此時(shí)，若希望擬合頻數(shù)分布曲線，可以利用SAS中CAPABILITY過(guò)程，見(jiàn)文獻(xiàn)［2］。

當(dāng)直方圖顯示有兩個(gè)或多個(gè)峰時(shí)，就稱(chēng)其為“多峰分布”。它們很可能是由多個(gè)“單峰分布”疊加而形成，被稱(chēng)為“有限混合分布樣本”。在實(shí)際問(wèn)題中，這種混合分布樣本很可能來(lái)自一個(gè)“不同質(zhì)”的總體，例如正常人樣本、某種疾病不同嚴(yán)重程度（輕、中、重度）的患者樣本。若希望擬合頻數(shù)分布曲線，可以利用SAS中FMM過(guò)程，見(jiàn)文獻(xiàn)［1］。

1.3 有限混合模型回歸分析的計(jì)算原理

1.3.1 概率分布

對(duì)于樣本資料而言，描述單組設(shè)計(jì)一元計(jì)量資料的頻數(shù)分布情況，所采用的方法被稱(chēng)為“擬合頻數(shù)分布曲線”；而對(duì)于總體資料而言，描述其一元連續(xù)性變量的變化規(guī)律，所采用的方法被稱(chēng)為“呈現(xiàn)概率分布密度函數(shù)”。在運(yùn)用統(tǒng)計(jì)學(xué)解決實(shí)際問(wèn)題時(shí)，通常都是通過(guò)樣本信息去推論總體的規(guī)律，故通常都是以“擬合頻數(shù)分布曲線”取代“呈現(xiàn)概率分布密度函數(shù)”。換言之，就是以“精確的概率分布密度函數(shù)”作為理論依據(jù)，描述“頻數(shù)分布曲線”的變化規(guī)律。

1.3.2 計(jì)算原理

式（1）描述的“頻數(shù)分布曲線”僅適用于“單組設(shè)計(jì)一元計(jì)量資料”，然而，當(dāng)資料中還有影響計(jì)量結(jié)果變量的自變量z和x時(shí)，需要將式（1）修改成式（2）的形式：

在式（2）中，有k個(gè)樣本混合在一起，k值需要結(jié)合專(zhuān)業(yè)知識(shí)事先給定；求和符號(hào)之后的第1項(xiàng)為樣本來(lái)自第j個(gè)總體的概率，其自變量為z；而其后的那一項(xiàng)為第 j個(gè)樣本的“頻數(shù)分布（對(duì)樣本而言）”或“概率密度函數(shù)（對(duì)總體而言）”，其自變量為x。式（2）中求和符號(hào)之后的第1項(xiàng)還應(yīng)滿足下式要求，見(jiàn)式（3）：在式（3）中，πj≥0，對(duì)于所有的 j（j=1到 k）都成立。

要想在式（3）的約束條件下求解式（2），涉及到“一般混合模型”的求解理論和技術(shù)方法，涉及到基于多種常用概率密度函數(shù)構(gòu)造“對(duì)數(shù)似然函數(shù)”并求解，還涉及到貝葉斯統(tǒng)計(jì)模擬算法［1，3-4］等深入的內(nèi)容，因篇幅所限，此處從略。

2 基于有限混合模型回歸分析解決實(shí)際問(wèn)題

2.1 問(wèn)題與數(shù)據(jù)結(jié)構(gòu)

【例1】下面是一個(gè)關(guān)于牛的喂食時(shí)間間隔的數(shù)據(jù)（計(jì)量數(shù)據(jù)）。按兩種情形對(duì)時(shí)間進(jìn)行劃分：第1種情形：分為3類(lèi)不同的時(shí)間段（相當(dāng)于前面所說(shuō)的混合樣本個(gè)數(shù)k=3），即：①各次進(jìn)食之間的時(shí)間間隔很短；②各次進(jìn)食之間的時(shí)間間隔稍長(zhǎng)一點(diǎn)，間隙讓牛飲水；③每?jī)纱芜M(jìn)食之間的時(shí)間間隔很長(zhǎng)。第2種情形：在上面的3類(lèi)時(shí)間間隔組中的每一種時(shí)間間隔內(nèi)，不同牛的進(jìn)食時(shí)間是不盡相同的。

測(cè)量141 414頭牛中每頭牛每?jī)纱芜M(jìn)食之間的“時(shí)間間隔數(shù)據(jù)（int）”，再取其對(duì)數(shù)，記為“l(fā)ogint”。以其為“計(jì)量結(jié)果變量（特別說(shuō)明：選取這樣的指標(biāo)作為結(jié)果變量，只有在特定專(zhuān)業(yè)領(lǐng)域內(nèi)才有意義，在一般的實(shí)際問(wèn)題中，‘時(shí)間間隔’是不可能用作結(jié)果變量的）”，由于原始數(shù)據(jù)的精確度很高，保留其精確度為“0.05”，這就導(dǎo)致了相同的數(shù)據(jù)很多，于是，可用類(lèi)似“頻數(shù)表”簡(jiǎn)化地呈現(xiàn)原始數(shù)據(jù)，其數(shù)據(jù)格式見(jiàn)表1。

表1 141 414頭牛中每頭牛每?jī)纱芜M(jìn)食之間的“時(shí)間間隔數(shù)據(jù)（logint）”

【對(duì)數(shù)據(jù)結(jié)構(gòu)的分析】以上采用“頻數(shù)分布表”形式呈現(xiàn)了資料，而資料的原始形式很簡(jiǎn)單，即一個(gè)計(jì)量變量“l(fā)ogint”，它有141 414個(gè)取值。通常這樣的數(shù)據(jù)被稱(chēng)為“單組設(shè)計(jì)一元計(jì)量資料”；而在本例中，因所有數(shù)據(jù)分別屬于“3類(lèi)不同的時(shí)間段”。也就是說(shuō)，若用一個(gè)變量k代表“3類(lèi)不同的時(shí)間段”，將此變量“k”及其具體取值也體現(xiàn)在每個(gè)“l(fā)ogint”數(shù)據(jù)之前，相當(dāng)于多了一個(gè)“分組變量”。此時(shí)，全部數(shù)據(jù)就可被視為“單因素三水平設(shè)計(jì)一元計(jì)量資料”了。

【統(tǒng)計(jì)分析方法的選擇】若希望采用 “單因素三水平設(shè)計(jì)一元計(jì)量資料”方差分析處理此資料，就必須在資料中全面反映出變量k及其取值；而在本例中，統(tǒng)計(jì)分析目的是 “擬合頻數(shù)分布曲線”，就只需要告知有 “k個(gè)樣本” （注意：k必須是一個(gè)具體的正整數(shù)），并且還需要告知這k個(gè)樣本所代表的 “分布”分別是什么。與特定分布對(duì)應(yīng)的“參數(shù)”可以告知，也可以不告知。例如 “dist=normal k=2 parms（3 1，5 1）”，這就是告知：有2個(gè)正態(tài)分布的樣本，它們的均值分別為3和5、方差分別為1和1；又例如 “dist=normal k=2”，這就是告知：有2個(gè)正態(tài)分布的樣本，它們的均值和方差都沒(méi)有指定，由統(tǒng)計(jì)軟件根據(jù)實(shí)際數(shù)據(jù)去估算。

2.2 創(chuàng)建SAS數(shù)據(jù)集

創(chuàng)建一個(gè)名為“cattle”的臨時(shí)SAS數(shù)據(jù)集的SAS數(shù)據(jù)步程序：

data cattle；

input LogInt Count＠＠；

datalines；

0.70 195 1.10 233 1.40 355 1.60 563

1.80 822 1.95 926 2.10 1018 2.20 1712

2.30 3190 2.40 2212 2.50 1692 2.55 1558

2.65 1622 2.70 1637 2.75 1568 2.85 1599

2.90 1575 2.95 1526 3.00 1537 3.05 1561

3.10 1555 3.15 1427 3.20 2852 3.25 1396

3.30 1343 3.35 2473 3.40 1310 3.45 2453

3.50 1168 3.55 2300 3.60 2174 3.65 2050

3.70 1926 3.75 1849 3.80 1687 3.85 2416

3.90 1449 3.95 2095 4.00 1278 4.05 1864

4.10 1672 4.15 2104 4.20 1443 4.25 1341

4.30 1685 4.35 1445 4.40 1369 4.45 1284

4.50 1523 4.55 1367 4.60 1027 4.65 1491

4.70 1057 4.75 1155 4.80 1095 4.85 1019

4.90 1158 4.95 1088 5.00 1075 5.05 912

5.10 1073 5.15 803 5.20 924 5.25 916

5.30 784 5.35 751 5.40 766 5.45 833

5.50 748 5.55 725 5.60 674 5.65 690

5.70 659 5.75 695 5.80 529 5.85 639

5.90 580 5.95 557 6.00 524 6.05 473

6.10 538 6.15 444 6.20 456 6.25 453

6.30 374 6.35 406 6.40 409 6.45 371

6.50 320 6.55 334 6.60 353 6.65 305

6.70 302 6.75 301 6.80 263 6.85 218

6.90 255 6.95 240 7.00 219 7.05 202

7.10 192 7.15 180 7.20 162 7.25 126

7.30 148 7.35 173 7.40 142 7.45 163

7.50 152 7.55 149 7.60 139 7.65 161

7.70 174 7.75 179 7.80 188 7.85 239

7.90 225 7.95 213 8.00 235 8.05 256

8.10 272 8.15 290 8.20 320 8.25 355

8.30 307 8.35 311 8.40 317 8.45 335

8.50 369 8.55 365 8.60 365 8.65 396

8.70 419 8.75 467 8.80 468 8.85 515

8.90 558 8.95 623 9.00 712 9.05 716

9.10 829 9.15 803 9.20 834 9.25 856

9.30 838 9.35 842 9.40 826 9.45 834

9.50 798 9.55 801 9.60 780 9.65 849

9.70 779 9.75 737 9.80 683 9.85 686

9.90 626 9.95 582 10.00 522 10.05 450

10.10 443 10.15 375 10.20 342 10.25 285

10.30 254 10.35 231 10.40 195 10.45 186

10.50 143 10.55 100 10.60 73 10.65 49

10.70 28 10.75 36 10.80 16 10.85 9

10.90 5 10.95 6 11.00 4 11.05 1

11.15 1 11.25 4 11.30 2 11.35 5

11.40 4 11.45 3 11.50 1

；

run；

2.3 利用SAS／STAT中KDE過(guò)程繪制資料的頻數(shù)分布曲線

利用下面的SAS過(guò)程步程序，可以繪制出反映計(jì)量變量logint的頻數(shù)分布直方圖和頻數(shù)分布曲線圖。

ods graphics on；

proc kde data=cattle；

univar LogInt／bwm=4；

freq count；

run；

【SAS主要輸出結(jié)果】

圖1 本例資料的頻數(shù)分布直方圖與頻數(shù)分布曲線圖

圖1 給人的印象是由兩個(gè)頻數(shù)分布混合而成的，然而，由專(zhuān)業(yè)知識(shí)可知，它實(shí)際上是由三個(gè)頻數(shù)分布混合而成的。

2.4 利用SAS／STAT中FMM過(guò)程擬合“三分量頻數(shù)分布曲線”

“三分量頻數(shù)分布曲線”就是擬合“由三個(gè)不同分布樣本混合而成的一個(gè)總樣本”的頻數(shù)分布曲線。程序?qū)⒔o出三條各自的頻數(shù)分布曲線以及一條混合的頻數(shù)分布曲線。

領(lǐng)域?qū)＜医o出的三個(gè)分布分別為：①正態(tài)分布，N（3，12）；②正態(tài)分布，N（5，12）；③威布爾分布。利用下面的SAS過(guò)程步程序，可以擬合并繪制出反映計(jì)量變量logint的頻數(shù)分布曲線圖。

proc fmm data=cattle gconv=0；

model LogInt=／dist=normal k=2 parms（3 1，5 1）；

model+ ／dist=weibull；

freq count；

run；

【SAS主要輸出結(jié)果及解釋】

Fit Statistics

以上是擬合統(tǒng)計(jì)量的有關(guān)計(jì)算結(jié)果：前5行都是關(guān)于擬合效果的評(píng)價(jià)指標(biāo)及其取值，這些數(shù)值只有在兩個(gè)或多個(gè)模型比較時(shí)才有參考的價(jià)值。

Parameter Estimates for Normal Model

以上是基于“正態(tài)分布”假定條件下，計(jì)算出兩個(gè)正態(tài)分布對(duì)應(yīng)的“參數(shù)估計(jì)”結(jié)果。因?yàn)閷?duì)于每一個(gè)特定的正態(tài)分布而言，只要給定了“均值”與“標(biāo)準(zhǔn)差（或方差）”，該正態(tài)分布也就唯一被確定了。

第1個(gè)正態(tài)分布為：N（3.3415，0.6718）=N（3.3415，0.81932）；

第 2個(gè)正態(tài)分布為：N（4.8940，1.4497）=N（4.8940，1.20402）。

Parameter Estimates for Weibull Model

以上是基于“威布爾分布”假定條件下，計(jì)算出對(duì)應(yīng)的“參數(shù)估計(jì)”結(jié)果。

第3個(gè)威布爾分布為：W（α，β，δ），其中，α＞0為形狀參數(shù)，β＞0為尺度參數(shù)，δ≥0為位置參數(shù)。上面計(jì)算的結(jié)果為：α=exp（2.2531）=9.5174、β=0.0685、δ=0。

Parameter Estimates for Mixing Probabilities

以上輸出的是各分布在混合分布中出現(xiàn)的概率，第1個(gè)正態(tài)分布出現(xiàn)的概率為0.4545，第2個(gè)正態(tài)分布出現(xiàn)的概率為0.3435，而第3個(gè)威布爾分布出現(xiàn)的概率為1-（0.4545+0.3435）=0.2020。于是，就可以寫(xiě)出混合樣本的概率密度函數(shù)如下：

y＾=0.4545 N（3.3415，0.81932）+0.3435 N（4.8940，1.20402）+0.2020W（9.5174，0.0685，0）

說(shuō)明：上式中的y＾代表圖2中“混合樣本”頻數(shù)曲線上縱坐標(biāo)的估計(jì)值，僅當(dāng)給定了橫坐標(biāo)上變量的一個(gè)確定的取值，y＾才有一個(gè)具體的數(shù)值與其對(duì)應(yīng)，下同，不再贅述。

圖2 本例資料的頻數(shù)分布直方圖與擬合的頻數(shù)分布曲線圖之一

若利用下面的SAS程序，可以獲得與上面類(lèi)似的結(jié)果，但會(huì)有較明顯的變化：

proc fmm data=cattle gconv=0；

model LogInt=／dist=normal k=2；

model+ ／dist=weibull；

freq count；

run；

ods graphics off；

【SAS程序說(shuō)明】這段 SAS程序與前面那段SAS程序非常相似，其主要區(qū)別在于：前面指定了兩個(gè)正態(tài)分布的“均值”與“方差”，而現(xiàn)在這段SAS程序沒(méi)有指定參數(shù)的具體數(shù)值，完全由實(shí)際的樣本數(shù)據(jù)計(jì)算而得。

【SAS主要輸出結(jié)果如下】

Fit Statistics

以上是擬合統(tǒng)計(jì)量的有關(guān)計(jì)算結(jié)果：前5行都是關(guān)于擬合效果的評(píng)價(jià)指標(biāo)及其取值，與前面相同內(nèi)容作比較，AIC、AICC和BIC的數(shù)值（說(shuō)明：這些數(shù)值越小越好）都變大了，說(shuō)明現(xiàn)在的模型對(duì)資料的擬合效果有所下降。

Parameter Estimates for＇Normal＇Model

第1個(gè)正態(tài)分布為：N（9.2883，0.4158）=N（9.2883，0.64482）；

第2個(gè)正態(tài)分布為：N（4.9106，1.7410）=N（4.9106，1.31952）。

Parameter Estimates for＇Weibull＇Model

第3個(gè)威布爾分布為：W（α，β，δ），其中，α＞0為形狀參數(shù)，β＞0為尺度參數(shù)，δ≥0為位置參數(shù)。上面計(jì)算的結(jié)果為：α=exp（1.2908）=3.6358、β=0.2093、δ=0。

Parameter Estimates for Mixing Probabilities

以上輸出的是各分布在混合分布中出現(xiàn)的概率，第1個(gè)正態(tài)分布出現(xiàn)的概率為0.1902、第2個(gè)正態(tài)分布出現(xiàn)的概率為0.3745，而第3個(gè)威布爾分布出現(xiàn)的概率為1-（0.1902+0.3745）=0.4353。于是，就可以寫(xiě)出混合樣本的概率密度函數(shù)如下：

y＾=0.1902 N（9.2883，0.64482）+0.3745 N（4.9106，1.31952）+0.4353W（3.6358，0.2093，0）

文獻(xiàn)［5］中有一個(gè)關(guān)于“1 000例血清谷丙轉(zhuǎn)氨酶（SGPT）的資料”，請(qǐng)感興趣的讀者運(yùn)用本文提供的SAS程序?qū)Α盎祀s樣本（包括‘非肝病患者’與‘肝病患者’）”進(jìn)行剖分，并與此文獻(xiàn)中基于“G-C級(jí)數(shù)”剖分的結(jié)果進(jìn)行比較。