胡良平

(1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029

1 概 述

1.1 有限混合模型

有限混合模型[1]見式(1)：

(1)

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

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

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

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

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

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

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

1.3.1 概率分布

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

1.3.2 計算原理

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

(2)

在式(2)中，有k個樣本混合在一起，k值需要結(jié)合專業(yè)知識事先給定；求和符號之后的第1項(xiàng)為樣本來自第j個總體的概率，其自變量為z；而其后的那一項(xiàng)為第j個樣本的“頻數(shù)分布(對樣本而言)”或“概率密度函數(shù)(對總體而言)”，其自變量為x。式(2)中求和符號之后的第1項(xiàng)還應(yīng)滿足下式要求，見式(3)：

(3)

在式(3)中，πj≥0，對于所有的j(j=1到k)都成立。

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

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

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

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

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

表1 141 414頭牛中每頭牛每兩次進(jìn)食之間的“時間間隔數(shù)據(jù)(logint)”

注：表中僅列出了極少部分?jǐn)?shù)據(jù)，詳細(xì)數(shù)據(jù)見后面SAS程序

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

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

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

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

data cattle;

input LogInt Count @@;

datalines;

0.701951.102331.403551.60563 1.808221.959262.1010182.201712 2.3031902.4022122.5016922.551558 2.6516222.7016372.7515682.851599 2.9015752.9515263.0015373.051561 3.1015553.1514273.2028523.251396 3.3013433.3524733.4013103.452453 3.5011683.5523003.6021743.652050 3.7019263.7518493.8016873.852416 3.9014493.9520954.0012784.051864 4.1016724.1521044.2014434.251341 4.3016854.3514454.4013694.451284 4.5015234.5513674.6010274.651491 4.7010574.7511554.8010954.851019 4.9011584.9510885.0010755.05912 5.1010735.158035.209245.25916 5.307845.357515.407665.45833 5.507485.557255.606745.65690 5.706595.756955.805295.85639 5.905805.955576.005246.05473 6.105386.154446.204566.25453 6.303746.354066.404096.45371 6.503206.553346.603536.65305 6.703026.753016.802636.85218 6.902556.952407.002197.05202 7.101927.151807.201627.25126 7.301487.351737.401427.45163 7.501527.551497.601397.65161 7.701747.751797.801887.85239 7.902257.952138.002358.05256 8.102728.152908.203208.25355 8.303078.353118.403178.45335 8.503698.553658.603658.65396 8.704198.754678.804688.85515 8.905588.956239.007129.05716 9.108299.158039.208349.25856 9.308389.358429.408269.45834 9.507989.558019.607809.65849 9.707799.757379.806839.85686 9.906269.9558210.0052210.0545010.1044310.1537510.2034210.2528510.3025410.3523110.4019510.4518610.5014310.5510010.607310.654910.702810.753610.801610.85910.90510.95611.00411.05111.15111.25411.30211.35511.40411.45311.501

;

run;

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

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

ods graphics on;

proc kde data=cattle;

univar LogInt / bwm=4;

freq count;

run;

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

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

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

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

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

領(lǐng)域?qū)＜医o出的三個分布分別為：①正態(tài)分布，N(3,12)；②正態(tài)分布，N(5,12)；③威布爾分布。利用下面的SAS過程步程序，可以擬合并繪制出反映計量變量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-2 Log Likelihood563153AIC (smaller is better)563169AICC (smaller is better)563169BIC (smaller is better)563248Pearson Statistic141458Effective Parameters8Effective Components3

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

Parameter Estimates for Normal Model成分參數(shù)估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|1Intercept3.34150.01260265.16<.00012Intercept4.89400.0544789.84<.00011Variance0.67180.012872Variance1.44970.05247

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

第1個正態(tài)分布為：N(3.3415,0.6718)=N(3.3415,0.81932)；

第2個正態(tài)分布為：N(4.8940,1.4497)=N(4.8940,1.20402)。

Parameter Estimates for Weibull Model成分參數(shù)估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|逆關(guān)聯(lián)估計3Intercept2.25310.0005064452.11<.00019.51743Scale0.068480.000427

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

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

Parameter Estimates for Mixing Probabilities成分參數(shù)鏈接尺度估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|概率1Probability0.81060.0340923.78<.00010.45452Probability0.53050.0464011.43<.00010.3435

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

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

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

proc fmm data=cattle gconv=0;

model LogInt = / dist=normal k=2;

model + / dist=weibull;

freq count;

run;

ods graphics off;

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

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

Fit Statistics-2 Log Likelihood564431AIC (smaller is better)564447AICC (smaller is better)564447BIC (smaller is better)564526Pearson Statistic141228Effective Parameters8Effective Components3

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

Parameter Estimates for 'Normal' Model成分參數(shù)估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|1Intercept9.28830.0050311846.28<.00012Intercept4.91060.02604188.56<.00011Variance0.41580.0050862Variance1.74100.02753

第1個正態(tài)分布為：N(9.2883,0.4158)=N(9.2883,0.64482)；

第2個正態(tài)分布為：N(4.9106,1.7410)=N(4.9106,1.31952)。

Parameter Estimates for ' Weibull' Model成分參數(shù)估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|逆關(guān)聯(lián)估計3Intercept1.29080.002790462.71<.00013.63583Scale0.20930.001311

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

Parameter Estimates for Mixing Probabilities成分參數(shù)鏈接尺度估計值標(biāo)準(zhǔn)誤差z 值Pr > |z|概率1Probability-0.82800.01922-43.08<.00010.19022Probability-0.15050.03678-4.09<.00010.3745

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

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

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