有限混合模型回歸分析

胡良平

（1.軍事科學院研究生院，北京 100850；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029 *通信作者：胡良平，E-mail：lphu812＠sina.com）

1 概 述

1.1 有限混合模型

有限混合模型［1］見式（1）：

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

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

1.2 有限混合模型回歸分析應用的場合

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

當直方圖顯示只有一個高峰時，就稱其為“單峰分布”。此時，可以根據高峰所處的位置，分為“正偏態分布（高峰位于左側）”“對稱分布（高峰居中，其特例為正態分布）”和“負偏態分布（高峰位于右側）”。此時，若希望擬合頻數分布曲線，可以利用SAS中CAPABILITY過程，見文獻［2］。

當直方圖顯示有兩個或多個峰時，就稱其為“多峰分布”。它們很可能是由多個“單峰分布”疊加而形成，被稱為“有限混合分布樣本”。在實際問題中，這種混合分布樣本很可能來自一個“不同質”的總體，例如正常人樣本、某種疾病不同嚴重程度（輕、中、重度）的患者樣本。若希望擬合頻數分布曲線，可以利用SAS中FMM過程，見文獻［1］。

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

1.3.1 概率分布

對于樣本資料而言，描述單組設計一元計量資料的頻數分布情況，所采用的方法被稱為“擬合頻數分布曲線”；而對于總體資料而言，描述其一元連續性變量的變化規律，所采用的方法被稱為“呈現概率分布密度函數”。在運用統計學解決實際問題時，通常都是通過樣本信息去推論總體的規律，故通常都是以“擬合頻數分布曲線”取代“呈現概率分布密度函數”。換言之，就是以“精確的概率分布密度函數”作為理論依據，描述“頻數分布曲線”的變化規律。

1.3.2 計算原理

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

在式（2）中，有k個樣本混合在一起，k值需要結合專業知識事先給定；求和符號之后的第1項為樣本來自第j個總體的概率，其自變量為z；而其后的那一項為第 j個樣本的“頻數分布（對樣本而言）”或“概率密度函數（對總體而言）”，其自變量為x。式（2）中求和符號之后的第1項還應滿足下式要求，見式（3）：在式（3）中，πj≥0，對于所有的 j（j=1到 k）都成立。

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

2 基于有限混合模型回歸分析解決實際問題

2.1 問題與數據結構

【例1】下面是一個關于牛的喂食時間間隔的數據（計量數據）。按兩種情形對時間進行劃分：第1種情形：分為3類不同的時間段（相當于前面所說的混合樣本個數k=3），即：①各次進食之間的時間間隔很短；②各次進食之間的時間間隔稍長一點，間隙讓牛飲水；③每兩次進食之間的時間間隔很長。第2種情形：在上面的3類時間間隔組中的每一種時間間隔內，不同牛的進食時間是不盡相同的。

測量141 414頭牛中每頭牛每兩次進食之間的“時間間隔數據（int）”，再取其對數，記為“logint”。以其為“計量結果變量（特別說明：選取這樣的指標作為結果變量，只有在特定專業領域內才有意義，在一般的實際問題中，‘時間間隔’是不可能用作結果變量的）”，由于原始數據的精確度很高，保留其精確度為“0.05”，這就導致了相同的數據很多，于是，可用類似“頻數表”簡化地呈現原始數據，其數據格式見表1。

表1 141 414頭牛中每頭牛每兩次進食之間的“時間間隔數據（logint）”

【對數據結構的分析】以上采用“頻數分布表”形式呈現了資料，而資料的原始形式很簡單，即一個計量變量“logint”，它有141 414個取值。通常這樣的數據被稱為“單組設計一元計量資料”；而在本例中，因所有數據分別屬于“3類不同的時間段”。也就是說，若用一個變量k代表“3類不同的時間段”，將此變量“k”及其具體取值也體現在每個“logint”數據之前，相當于多了一個“分組變量”。此時，全部數據就可被視為“單因素三水平設計一元計量資料”了。

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

2.2 創建SAS數據集

創建一個名為“cattle”的臨時SAS數據集的SAS數據步程序：

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過程繪制資料的頻數分布曲線

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

ods graphics on；

proc kde data=cattle；

univar LogInt／bwm=4；

freq count；

run；

【SAS主要輸出結果】

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

圖1 給人的印象是由兩個頻數分布混合而成的，然而，由專業知識可知，它實際上是由三個頻數分布混合而成的。

2.4 利用SAS／STAT中FMM過程擬合“三分量頻數分布曲線”

“三分量頻數分布曲線”就是擬合“由三個不同分布樣本混合而成的一個總樣本”的頻數分布曲線。程序將給出三條各自的頻數分布曲線以及一條混合的頻數分布曲線。

領域專家給出的三個分布分別為：①正態分布，N（3，12）；②正態分布，N（5，12）；③威布爾分布。利用下面的SAS過程步程序，可以擬合并繪制出反映計量變量logint的頻數分布曲線圖。

proc fmm data=cattle gconv=0；

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

model+ ／dist=weibull；

freq count；

run；

【SAS主要輸出結果及解釋】

Fit Statistics

以上是擬合統計量的有關計算結果：前5行都是關于擬合效果的評價指標及其取值，這些數值只有在兩個或多個模型比較時才有參考的價值。

Parameter Estimates for Normal Model

以上是基于“正態分布”假定條件下，計算出兩個正態分布對應的“參數估計”結果。因為對于每一個特定的正態分布而言，只要給定了“均值”與“標準差（或方差）”，該正態分布也就唯一被確定了。

第1個正態分布為：N（3.3415，0.6718）=N（3.3415，0.81932）；

第 2個正態分布為：N（4.8940，1.4497）=N（4.8940，1.20402）。

Parameter Estimates for Weibull Model

以上是基于“威布爾分布”假定條件下，計算出對應的“參數估計”結果。

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

Parameter Estimates for Mixing Probabilities

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

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

說明：上式中的y＾代表圖2中“混合樣本”頻數曲線上縱坐標的估計值，僅當給定了橫坐標上變量的一個確定的取值，y＾才有一個具體的數值與其對應，下同，不再贅述。

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

若利用下面的SAS程序，可以獲得與上面類似的結果，但會有較明顯的變化：

proc fmm data=cattle gconv=0；

model LogInt=／dist=normal k=2；

model+ ／dist=weibull；

freq count；

run；

ods graphics off；

【SAS程序說明】這段 SAS程序與前面那段SAS程序非常相似，其主要區別在于：前面指定了兩個正態分布的“均值”與“方差”，而現在這段SAS程序沒有指定參數的具體數值，完全由實際的樣本數據計算而得。

【SAS主要輸出結果如下】

Fit Statistics

以上是擬合統計量的有關計算結果：前5行都是關于擬合效果的評價指標及其取值，與前面相同內容作比較，AIC、AICC和BIC的數值（說明：這些數值越小越好）都變大了，說明現在的模型對資料的擬合效果有所下降。

Parameter Estimates for＇Normal＇Model

第1個正態分布為：N（9.2883，0.4158）=N（9.2883，0.64482）；

第2個正態分布為：N（4.9106，1.7410）=N（4.9106，1.31952）。

Parameter Estimates for＇Weibull＇Model

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

Parameter Estimates for Mixing Probabilities

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

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

文獻［5］中有一個關于“1 000例血清谷丙轉氨酶（SGPT）的資料”，請感興趣的讀者運用本文提供的SAS程序對“混雜樣本（包括‘非肝病患者’與‘肝病患者’）”進行剖分，并與此文獻中基于“G-C級數”剖分的結果進行比較。