用mixed模型解決單因素方差分析中方差不齊時的隨機模擬及SAS實現(xiàn)*

復旦大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計與社會醫(yī)學教研室公共安全教育部重點實驗室（200032） 肖林海 趙耐青

用mixed模型解決單因素方差分析中方差不齊時的隨機模擬及SAS實現(xiàn)*

復旦大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計與社會醫(yī)學教研室公共安全教育部重點實驗室（200032） 肖林海 趙耐青△

目的用m ixed模型處理完全隨機設計的三組資料方差不齊時的I類錯誤率和對各組總體均值的估計，并與Kruskal-Wallis秩和檢驗和ANOVA方法比較檢驗的把握度。方法利用SAS9.3產(chǎn)生模擬數(shù)據(jù)，并分別用proc mixed過程、proc npar1way過程和proc glm過程進行m ixed模型參數(shù)估計、Kruskal-Wallis秩和檢驗和one-way ANOVA分析。模擬試驗重復1000次。結(jié)果方差不齊時，對不均衡資料，m ixed法的檢驗把握度要高于Kruskal-Wallis秩和檢驗和one-way ANOVA法。對均衡資料，三種方法的檢驗效能接近，但后兩者的Power受方差不齊的影響較大。結(jié)論研究設計時盡量使各組樣本均衡或接近則三種方法的I類錯誤率最小，Power達到最大。對不均衡資料方差不齊時，m ixed法的檢驗把握度要高于Kruskal-Wallis秩和檢驗和one-way ANOVA法，且其I類錯誤率能得到控制。

方差不齊 mixed模型 Kruskal-Wallis秩和檢驗 把握度

在一般的醫(yī)學統(tǒng)計教科書中，當各組資料滿足獨立、正態(tài)和方差齊性時采用的是方差分析；當資料不服從正態(tài)分布或組間方差不齊時，通常考慮非參數(shù)的檢驗方法，如Kruskal-Wallis秩和檢驗（以下簡稱K-W法）。

本文的研究是基于完全隨機設計的三組資料進行均數(shù)間差異的比較，在滿足了各組資料相互獨立、服從正態(tài)分布的條件下，考慮方差不齊時用mixed方法分析，并綜合考慮樣本量對I類錯誤發(fā)生率的影響，及與Kruskal-Wallis秩和檢驗和ANOVA法比較檢驗的把握度。

研究方法

1.m ixed模型［1］

依據(jù)背景，所采用的統(tǒng)計分析模型如下：

yij＝μ＋τi＋εij，其中i＝1，2，3，表示組別；j＝1，2，…，ni，表示各組第j個體，每組對應有ni個體；yij表示第i組第j個體的觀測值，μ表示總體均數(shù)，τi表示第i組總體均數(shù)與第3組總體均數(shù)的差異（i＝1，2，τ3＝0），εij表示第i組第j個體的觀測誤差，則εij～N表示第i組的總體方差。

在SAS9.3中，利用procm ixed［2］過程即可進行編程分析，所用程序如下：

2.模擬研究

利用SAS9.3產(chǎn)生服從正態(tài)分布的模擬數(shù)據(jù)，分別采用mixed模型、Kruskal-Wallis秩和檢驗和ANOVA法進行分析，并考慮在樣本量均衡和各組樣本量不等的情況下，該mixed法的I類錯誤率和檢驗效能，并與Kruskal-Wallis秩和檢驗和ANOVA法作比較。

為便于模擬，本研究擬將樣本量設定為n1＝n2＝n，以n∶n3分別為1∶1，1∶2，1∶3，1∶4和1∶5進行模擬，各組總體方差亦按照兩種情境進行模擬。重復模擬試驗1000次。

模擬結(jié)果

1.總體方差相同時，mixed方法的Ⅰ類錯誤發(fā)生率

參數(shù)的設定：各組總體均值設為μ1＝μ2＝μ3＝20，對應mixed模型中的τ1＝τ2＝0。各組總體標準差設為σ1＝σ2＝σ3＝1。各組樣本量的設置及相應模擬結(jié)果見表1。

表1 總體方差相同時，樣本量對I類錯誤率的影響（檢驗水準α＝0.05）

由表1中結(jié)果可以看出，樣本量均衡時，三種方法的Ⅰ類錯誤率均很低，而隨著組間樣本量差異變大，則檢驗的Ⅰ類錯誤率均增高，當n∶n3≤1∶4時，Ⅰ類錯誤率無法控制。由此結(jié)果也可以看出，當各組總體方差相同時，ANOVA法相對Kruskal-Wallis秩和檢驗的I類錯誤率稍低，而m ixed法的一類錯誤率較高，因此組間總體方差相同且資料不均衡時，mixed方法不是最優(yōu)的。

2.總體方差相同時，m ixed方法的檢驗效能［3］

參數(shù)的設定：各組總體均值設為μ1＝20，μ2＝21，μ3＝22，對應mixed模型中的τ1＝-2，τ2＝-1。各組總體標準差設為σ1＝σ2＝σ3＝1。各組樣本量的設置及模擬結(jié)果見表2。

表2 總體方差相同時，樣本量對檢驗效能的影響（檢驗水準α＝0.05）

由表2結(jié)果可知，總體方差相同時，三種方法的檢驗效能均較高，且即使各組間樣本量的不均衡性增大，三種方法的把握度并沒有顯著變化。

綜上，由表1和表2的結(jié)果可知，各組總體方差相同時，用mixed法、ANOVA法和Kruskal-Wallis秩和檢驗進行各組總體均數(shù)的比較時，三種方法的檢驗效能都很高，但組間樣本量不均衡時，m ixed法的Ⅰ類錯誤率較高，ANOVA法最低。因此各組總體方差相同時，應優(yōu)先考慮ANOVA法，而不建議采用mixed法。

3.總體方差不同時，mixed方法的Ⅰ類錯誤率

各組總體均值設為μ1＝μ2＝μ3＝20，各組總體標準差分別為σ1＝σ2＜σ3和σ1＜σ2＜σ3的情景進行模擬。各組樣本量和參數(shù)的設置及相應模擬結(jié)果見表3

由表3結(jié)果可以看出，樣本量均衡時，三種方法的Ⅰ類錯誤率均很低，但隨著組間方差差異的增大，I類錯誤率有所增高。組間樣本量不均衡時，mixed法檢驗的I類錯誤率最高，但都在可控制的范圍內(nèi)。而對不均衡資料，固定各組樣本量，各組總體方差差異增大時，三種方法的I類錯誤率略有降低。

表3 總體方差不同時，樣本量和方差不齊對I類錯誤率的影響（檢驗水準α＝0.05）

4.總體方差不同時，m ixed方法的把握度

設定各組總體均值分別為μ1＝20，μ2＝21，μ3＝22。各組樣本量和參數(shù)的設置及相應模擬結(jié)果見表4。

表4 總體方差不同時，樣本量和方差不齊對把握度的影響（檢驗水準α＝0.05）

由表4結(jié)果可以看出，在樣本量均衡時，不管方差的差異性如何，三種方法的Power接近，但方差的差異性越明顯，則三種方法的Power均顯著下降。而資料不均衡時，在方差的差異性較小時，三種方法的把握度均較高，但mixed法的檢驗效能要高于Kruskal-Wallis秩和檢驗和ANOVA法，當方差的差異性增大時，m ixed法的檢驗效能略有下降，而Kruskal-Wallis秩和檢驗和ANOVA法的檢驗效能下降明顯。

結(jié)論和討論

基于以上各種情境的模擬結(jié)果，本研究可以得到如下結(jié)論：

1.不管采用何種統(tǒng)計分析方法，研究設計階段盡量使各組樣本量均衡，才可使統(tǒng)計分析的I類錯誤率降低和把握度達到最大。

2.當各組樣本不均衡時，統(tǒng)計分析的把握度會降低，并且Ⅰ類錯誤率增高。在本研究中，如表1示，當n∶n3＝1∶4時，三種方法的Ⅰ類錯誤率已超過0.05。因此，即使樣本不均衡，考慮到統(tǒng)計檢驗的I類錯誤率和把握度，也應將樣本量的差異控制在一定范圍。

3.對各組樣本不均衡資料，mixed方法的檢驗效能要高于Kruskal-Wallis秩和檢驗和ANOVA法，尤其在方差不齊時，mixed方法的檢驗效能顯著高于后兩者。當組間方差差異增大時，如本例由其中兩組方差相同到三組方差均不同時，m ixed方法檢驗的把握度僅略有下降，而Kruskal-Wallis秩和檢驗和ANOVA法的把握度下降明顯，說明不均衡資料固定各組樣本量后，mixed方法的把握度受方差不齊的影響較小，而秩和檢驗和ANOVA法對方差不齊的變化較為敏感。

綜上，各組總體方差相等時，對不均衡資料的完全隨機設計，mixed方法的Ⅰ類錯誤率高于Kruskal-Wallis秩和檢驗和ANOVA法而無法控制。因此在各組總體方差相等時，不建議采用m ixed方法，而建議采用ANOVA法。

各組總體方差不齊時，mixed法的Ⅰ類錯誤率較低，且mixed法的把握度受各組方差不齊程度的變化不敏感，而Kruskal-Wallis秩和檢驗和ANOVA法受方差不齊的影響較大。另外，m ixed法還能給出各組總體均值和總體方差的估計，相對于Kruskal-Wallis秩和檢驗和ANOVA法有著明顯的優(yōu)勢。但是，當各組樣本量相差較大時，mixed法的Ⅰ類錯誤率增大，因此試驗設計時應盡量使各組樣本均衡。

1.Brown HAPR.Applied M ixed Models in Medicine.New York：John W iley＆Sons，1999.

2.SAS Institute Inc.2011.Base SAS?9.3 Procedures Guide.Cary，NC：SAS Institute Inc.

3.姚嵩坡，劉盛元，王濱有.假設檢驗中檢驗效能的計算及SAS實現(xiàn).中國衛(wèi)生統(tǒng)計，2010，27（4）：434-436.

（責任編輯：劉 壯）

*：國家自然科學基金（81273187）

△通信作者：趙耐青