混合效應(yīng)模型中的方差成分檢驗*

曾 平 趙 楊 陳 峰

混合效應(yīng)模型中的方差成分檢驗*

曾 平1趙 楊2陳 峰2

在很多科學(xué)問題中需要在混合效應(yīng)模型框架下對隨機(jī)效應(yīng)方差成分(暫記為τ2)進(jìn)行假設(shè)檢驗[1-6]，也即檢驗H0:τ2=0。除直接科學(xué)興趣外，許多間接醫(yī)學(xué)問題也能轉(zhuǎn)化為對方差成分的檢驗。例如，為判斷在懲罰樣條回歸中是參數(shù)模型還是非參數(shù)模型更合適，Claeskens[7]首先建立混合效應(yīng)模型，將模型選擇問題轉(zhuǎn)化為對隨機(jī)效應(yīng)方差成分是否等于零的假設(shè)檢驗問題，最后通過限制性似然比檢驗H0:τ2=0。其他研究者也用同樣的方法處理過類似問題[8-12]。然而，方差成分為非負(fù)參數(shù)，對方差成分的假設(shè)檢驗是非標(biāo)準(zhǔn)的：在H0下τ2位于參數(shù)空間邊緣。由于這種限制，常用的漸近χ2無效分布對似然比統(tǒng)計量不再成立[1-3,8]?；旌闲?yīng)模型中的方差成分檢驗吸引了廣泛研究興趣[7-8，13-19]。本文主要介紹方差成分檢驗的似然比方法，方差成分得分檢驗或Wald檢驗可參考其他相關(guān)文獻(xiàn)[13,20,21]。

線性混合效應(yīng)模型的方差成分檢驗

設(shè)Y為連續(xù)變量，X=(X1,…,Xp)和Z=(Z1,…,Zq)分別為n×p和n×q的矩陣，n為樣本量。在廣泛意義上Y、X和Z既適于非獨立數(shù)據(jù)結(jié)構(gòu)[22-25]，如重復(fù)測量數(shù)據(jù)或聚群數(shù)據(jù)[26-29]，也適于獨立數(shù)據(jù)情形[7-11,30,31]，如懲罰樣條回歸。更加詳細(xì)的例子可參考上述文獻(xiàn)。Y、X和Z之間的關(guān)系由線性混合效應(yīng)模型[25]刻畫

(1)

用矩陣的形式模型(1)可表達(dá)為

Y=Xβ+Zb+ε,b～N(0,τ2K),ε～N(0,σ2In)

(2)

式中β為固定效應(yīng)，b為隨機(jī)效應(yīng)，τ2為未知方差成分，K為已知q×q矩陣描述隨機(jī)效應(yīng)間相互關(guān)系，In為n維單位陣，σ2為殘差方差，與b不相關(guān)。模型中Y的期望和方差分別為

E(Y)=Xβ,Var(Y)=τ2ZKZ′+σ2In=σ2Vλ

(3)

其中Vλ=λZKZ′+In,λ=τ2/σ2。設(shè)參數(shù)為θ=(β,λ,σ2)，模型(2)的對數(shù)似然函數(shù)為

(4)

似然比統(tǒng)計量定義為

LRTn=2sup[H1(θ)-H0(θ)]

(5)

1.漸近混合卡方分布

2.基于譜分解的模擬分布

為了獲得似然比統(tǒng)計量的無效分布，Crainiceanu

圖1 模擬計算的似然比統(tǒng)計量和的分位數(shù)。

和Ruppert在譜分解基礎(chǔ)上提出了一種基于模擬抽樣算法[8]。假設(shè)模型(2)中λ已知，可獲得β和σ2的極大似然估計值

(6)

帶入對數(shù)似然函數(shù)，獲得剖面似然函數(shù)

(7)

定義

(8)

ξ是矩陣K1/2Z′ZK1/2的特征根，其中，

(9)

3.重抽樣方法

統(tǒng)計分析中當(dāng)面對復(fù)雜假設(shè)檢驗時訴諸于重抽樣方法是一種常用策略[34-39]。重抽樣避免了數(shù)學(xué)推導(dǎo)，概念上簡單并且容易操作，屬于計算密集型統(tǒng)計技術(shù)。Fitzmaurice等[14]以及Lee和Braun[40]提出了似然比方差成分的permutation檢驗，F(xiàn)araway[23]建立了似然比方差成分的參數(shù)bootstrap檢驗。其他研究者對似然比方差成分提出了類似的重抽樣檢驗方法[15-16]。

廣義線性混合效應(yīng)模型的方差成分檢驗

1.logistic混合效應(yīng)模型和PQL算法

目前幾乎所有似然比方差成分檢驗都是在線性混合效應(yīng)模型框架下完成，即都是針對連續(xù)應(yīng)變量的；而對離散應(yīng)變量似然比方差成分檢驗的研究還很少。以logistic回歸為例，接下來討論廣義混合效應(yīng)模型的似然比方差成分檢驗。設(shè)Y為二分類應(yīng)變量，X和Z同前。Y、X和Z之間的關(guān)系通過logistic混合效應(yīng)模型[22,41-42]描述

(10)

2[Lτ2≥0(β,τ2)-L(β,τ2=0)]

(11)

然而，由于對數(shù)似然函數(shù)涉及積分運算，實際中除極有限和簡單情況外，大多數(shù)時候都無法精確計算Lτ2≥0(β,τ2)和獲得精確的對數(shù)似然值。因此不能直接通過式(11)進(jìn)行似然比檢驗。

logistic混合模型的參數(shù)估計算法包括數(shù)值積分[43]，懲罰擬似然(penalized quasi-likelihood，PQL)和邊際擬似然(marginal quasi-likelihood，MQL)[42]，偽似然估計[41]，Monte Carlo估計[44-45]和Gibbs抽樣[46]。數(shù)值積分只適用于低維度的積分運算，一般很難處理超過五維的積分[22]；Monte Carlo估計和Gibbs抽樣的不足在于計算繁重。PQL和MQL能夠通過近似方法有效處理高維積分問題，計算相對簡單，且能夠通過重復(fù)調(diào)用已有的線性混合模型程序執(zhí)行。對logistic混合模型進(jìn)行方差成分的似然比檢驗至少包含以下的三個困難[47]：與線性混合模型不同，通常不能獲得關(guān)于logistic混合模型對數(shù)函數(shù)的閉型表達(dá)式，因此精確計算logistic混合模型的對數(shù)似然值和似然比統(tǒng)計量變得很困難；由于涉及高維計算運算，很難獲得logistic混合模型參數(shù)的精確估計值；即使能夠有效估計logistic混合模型和計算其對數(shù)似然比統(tǒng)計量，如何獲得該統(tǒng)計量的無效分布也不清楚。基于logistic混合模型已有的理論發(fā)展，本文僅僅專注其方差成分的似然比檢驗。

2.基于重抽樣方法的偽似然比方差成分檢驗

為了處理上面的問題，可以通過基于重抽樣的偽似然比方法進(jìn)行方差成分檢驗[47]。具體步驟如下：①采用PQL算法估計logistic混合模型參數(shù)；PQL算法可通過R軟件中的函數(shù)glmmPQL進(jìn)行，該函數(shù)包含在MASS軟件包中[48]；②當(dāng)PQL算法收斂時，有工作應(yīng)變量

(12)

和偽似然函數(shù)

(13)

其他問題和可能的研究方向

目前，由于以下幾方面的原因，似然比方差成分檢驗主要集中在線性混合效應(yīng)模型并且模型中只包含一個方差成分[8,12,17-18]：①能夠較為容易估計線性混合效應(yīng)模型的參數(shù)和計算其對數(shù)似然值，進(jìn)而計算似然比統(tǒng)計量；②多個方差成分共存時，不存在簡單的對數(shù)似然函數(shù)譜分解形式，因此很難獲得似然比統(tǒng)計量的無效分布；③由于可能的高維積分運算，一般難以精確估計廣義混合效應(yīng)模型的參數(shù)和計算其對數(shù)似然值，以及獲得似然比統(tǒng)計量的無效分布；④重抽樣方法屬于計算密集型統(tǒng)計方法，難以大規(guī)模運用。

在標(biāo)準(zhǔn)的參數(shù)假設(shè)檢驗中似然比檢驗被證明是在小樣本時最優(yōu)的[49]。事實上，由于似然比檢驗充分利用了H0條件下的模型和H1條件下的模型，相對得分檢驗和Wald檢驗，似然比檢驗被證明在非標(biāo)準(zhǔn)參數(shù)假設(shè)檢驗也具有更高統(tǒng)計效能[7,30,50]。因此，推廣和發(fā)展似然比檢驗以適合更加廣泛的情形具有理論和實際意義。具體而言，以下幾方面問題值得進(jìn)一步探索：①在線性混合效應(yīng)模型中當(dāng)存在多個方差成分而僅對其中某個成分感興趣時，Greven[12]等建議可以首先從應(yīng)變量中消除冗余隨機(jī)效應(yīng)；此時，尚包括應(yīng)變量的殘差和待檢驗方差成分，然后按照前文只包含單個方差成分的似然比檢驗進(jìn)行統(tǒng)計分析。然而，這樣做的缺點在于，沒有考慮到隨機(jī)效應(yīng)之間可能存在的相關(guān)，而且在模型估計多個方差成分時還可能存在數(shù)值問題，如模型不收斂或計算不穩(wěn)定。②在線性混合效應(yīng)模型中假設(shè)殘差是獨立的；當(dāng)存在相關(guān)協(xié)方差結(jié)構(gòu)，即ε～N(0,R(α))，R為一般形式協(xié)方差矩陣、包含未知參數(shù)向量α；在這種情況下可先估計殘差協(xié)方差，在此基礎(chǔ)上建立偽數(shù)據(jù)，然后將偽數(shù)據(jù)當(dāng)做原始數(shù)據(jù)按照線性混合效應(yīng)模型框架下相似的方法進(jìn)行似然比方差成分檢驗[18]。但這種方法在協(xié)方差矩陣R的選擇以及其適用性等方面需要更多的研究。③前文在廣義線性混合效應(yīng)模型框架下建立的偽似然比方差成分檢驗其有效性取決于PQL算法；PQL算法本身屬于有偏估計；雖然提出了校正PQL偏倚的算法[51-55]，有偏和無偏的參數(shù)估計算法對偽似然比方差成分檢驗的影響如何還需進(jìn)一步研究[47]。④當(dāng)前討論的方差成分檢驗都假設(shè)被檢驗的方差成分和其他冗余成分之間獨立，如何在考慮方差成分相關(guān)的前提下建立有效的似然比檢驗？⑤如何進(jìn)一步在廣義線性混合效應(yīng)模型框架下發(fā)展在計算上更加有效的似然比檢驗、而不僅僅依賴于重抽樣方法？⑥雖然似然比檢驗具有更高的統(tǒng)計效能，但是由于其需要同時估計H0條件下的模型和H1條件下的模型，所以相對得分檢驗，計算速度明顯較慢；因此，需要發(fā)展有效的計算方法提高似然比檢驗的計算速度，以適合大規(guī)模的應(yīng)用[8,12]。

[1]Self SG,Liang KY.Asymptotic Properties of Maximum Likelihood Estimators and Likelihood Ratio Tests under Nonstandard Conditions.J R Stat Soc B,1987,82(398):605-610.

[2]Stram DO,Lee JW.Variance Components Testing in the Longitudinal Mixed Effects Model.Biometrics,1994,50(4):1171-1177.

[3]Liang KY,Self SG.On the Asymptotic Behaviour of the Pseudolikelihood Ratio Test Statistic.J R Stat Soc B,1996,58(4):785-796.

[4]Lindquist MA,Spicer J,Asllani I,et al.Estimating and testing variance components in a multi-level GLM.Neuroimage,2012,59(1):490-501.

[5]Drikvandi R,Verbeke G,Khodadadi A,et al.Testing multiple variance components in linear mixed-effects models.Biostatistics,2013,14(1):144-159.

[6]Nobre J,Singer J,Sen P.U-tests for variance components in linear mixed models.Test,2013,22(4):580-605.

[7]Claeskens G.Restricted likelihood ratio lack-of-fit tests using mixed spline models.J R Stat Soc B,2004,66(4):909-926.

[8]Crainiceanu CM,Ruppert D.Likelihood ratio tests in linear mixed models with one variance component.J R Stat Soc B,2004,66(1):165-185.

[9]Crainiceanu CM,Ruppert D.Restricted likelihood ratio tests in nonparametric longitudinal models.Stat Sinica,2004,14(3):713-730.

[10]Crainiceanu CM,Ruppert D.Likelihood ratio tests for goodness-of-fit of a nonlinear regression model.J Multivariate Anal,2004,91(1):35-52.

[11]Crainiceanu C,Ruppert D,Claeskens G,et al.Exact likelihood ratio tests for penalised splines.Biometrika,2005,92(1):91-103.

[12]Greven S,Crainiceanu CM,Küchenhoff H,et al.Restricted Likelihood Ratio Testing for Zero Variance Components in Linear Mixed Models.J Comput Graph Stat,2008,17(4):870-891.

[13]Lin X.Variance component testing in generalised linear models with random effects.Biometrika,1997,84(2):309-326.

[14]Fitzmaurice GM,Lipsitz SR,Ibrahim JG.A Note on Permutation Tests for Variance Components in Multilevel Generalized Linear Mixed Models.Biometrics,2007,63(3):942-946.

[15]Sinha SK.Bootstrap tests for variance components in generalized linear mixed models.Scand J Stat,2009,37(2):219-234.

[16]Samuh MH,Grilli L,Rampichini C,et al.The Use of Permutation Tests for Variance Components in Linear Mixed Models.Commun Stat Theor M,2012,41(16-17):3020-3029.

[17]Staicu A-M,Li Y,Crainiceanu CM,et al.Likelihood ratio tests for dependent data with applications to longitudinal and functional data analysis.Scand J Stat,2014,41(4):932-949.

[18]Wiencierz A,Greven S,Küchenhoff H.Restricted likelihood ratio testing in linear mixed models with general error covariance structure.Electron J Stat,2011,5:1718-1734.

[19]Chen Y,Liang KY.On the asymptotic behaviour of the pseudolikelihood ratio test statistic with boundary problems.Biometrika,2010,97(3):603-620.

[20]Lin X,Zhang D.Inference in generalized additive mixed models by using smoothing splines.J Roy Stat Soc,B,1999,61(2):381-400.

[21]Zhang D,Lin X.Hypothesis testing in semiparametric additive mixed models.Biostatistics,2003,4(1):57-74.

[22]Diggle P,Heagerty P,Liang KY,et al.Analysis of longitudinal data,2nd Edition.New York:Oxford University Press,2002.

[23]Faraway JJ.Extending the linear model with R:generalized linear,mixed effects and nonparametric regression models.New York:Chapman & Hall,2005.

[24]Pinheiro JC,Bates D.Mixed-Effects Models in S and S-PLUS,2nd Edition.New York:Springer,2009.

[25]Laird NM,Ware JH.Random-effects models for longitudinal data.Biometrics,1982,38(4):963-974.

[26]張巖波,何大衛(wèi),劉桂芬,等.重復(fù)測量數(shù)據(jù)的混合模型及其MIXED過程實現(xiàn)-混合線性模型及其SAS軟件實現(xiàn).中國衛(wèi)生統(tǒng)計,2001,18(5):272-275.

[27]楊珉.多元分析的發(fā)展-多水平模型簡介.中國衛(wèi)生統(tǒng)計,1994,11(5):32-35.

[28]曾平,曹紅艷,劉桂芬.重復(fù)測量計數(shù)資料的隨機(jī)效應(yīng)ZIP模型.中國衛(wèi)生統(tǒng)計,2009,26(4):344-346,349.

[29]任仕泉,陳峰.多水平統(tǒng)計模型在多中心臨床試驗評價中的應(yīng)用.中國衛(wèi)生統(tǒng)計,2000,17(2):108-110.

[30]Zeng P,Zhao Y,Liu J,et al.Likelihood Ratio Tests in Rare Variant Detection for Continuous Phenotypes.Ann Hum Genet,2014,78(5):320-332.

[31]Wu MC,Lee S,Cai T,et al.Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test.Am J Hum Genet,2011,89(1):82-93.

[32]Lehmann EL,Romano JP.Testing statistical hypotheses,3rd Edition.New York:Springer,2006.

[33]Scheipl F,Greven S,Küchenhoff H.Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed models.Comput Stat Data Anal,2008,52(7):3283-3299.

[34]Davison AC,Hinkley DV.Bootstrap Methods and Their Application.Cambridge:Cambridge University Press,1997.

[35]Efron B,Tibshirani R.An Introduction to the Bootstrap.New York:Chapman and Hall,1993.

[36]Good P.Permutation,Parametric,and Bootstrap Tests of Hypotheses,3rd Edition:New York:Springer,2005.

[37]荀鵬程,趙楊,易洪剛,等.Permutation Test 在假設(shè)檢驗中的應(yīng)用.數(shù)理統(tǒng)計與管理,2006,25(5):616-621.

[38]陳峰,陸守曾,楊珉.Bootstrap 估計及其應(yīng)用.中國衛(wèi)生統(tǒng)計,1997,14(5):5-7.

[39]劉沛,劉元寶,王燦楠.膳食暴露評估模型及其構(gòu)建方法.中華預(yù)防醫(yī)學(xué)雜志,2007,41(6):502-504.

[40]Lee OE,Braun TM.Permutation Tests for Random Effects in Linear Mixed Models.Biometrics,2012,68(2):486-493.

[41]Wolfinger R,O′.Connell M.Generalized linear mixed models a pseudo-likelihood approach.J Stat Comput Sim,1993,48(3-4):233-243.

[42]Breslow NE,Clayton DG.Approximate Inference in Generalized Linear Mixed Models.J Am Stat Assoc,1993,88(421):9-25.

[43]曹紅艷,劉桂芬,曾平,等.預(yù)測性偽似然法和貝葉斯法廣義線性混合模型估計.中國藥物與臨床,2008,8(11):861-863.

[44]McCulloch CE.Maximum Likelihood Variance Components Estimation for Binary Data.J Am Stat Assoc,1994,89(425):330-335.

[45]McCulloch CE.Maximum likelihood algorithms for generalized linear mixed models.J Am Stat Assoc,1997,92(437):162-170.

[46]Zeger SL,Karim MR.Generalized linear models with random effects:a Gibbs sampling approach.J Am Stat Assoc,1991,86(413):79-86.

[47]Zeng P,Zhao Y,Chen F.Permutation-based variance component test in generalized linear mixed model with application to multilocus genetic association study.BMC Med Res Methodol,2015,15:37.

[48]Venables WN,Ripley BD.Modern Applied Statistics with S,4th Edition.New York:Springer,2002.

[49]Donoho D,Jin J.Higher criticism for detecting sparse heterogeneous mixtures.Ann Stat,2004,32(3):962-994.

[50]Zeng P,Zhao Y,Zhang L,et al.Rare Variants Detection with Kernel Machine Learning Based on Likelihood Ratio Test.PLoS ONE,2014,9(3):e93355.

[51]Lin X,Breslow NE.Bias Correction in Generalized Linear Mixed Models with Multiple Components of Dispersion.J Am Stat Assoc,1996,91(435):1007-1016.

[52]Wang N,Lin X,Gutierrez RG,et al.Bias Analysis and SIMEX Approach in Generalized Linear Mixed Measurement Error Models.J Am Stat Assoc,1998,93(441):249-261.

[53]Pawitan Y.Two-staged estimation of variance components in generalized linear mixed models.J Stat Comput Sim,2001,69(1):1-17.

[54]Breslow NE,Lin X.Bias correction in generalised linear mixed models with a single component of dispersion.Biometrika,1995,82(1):81-91.

[55]Jiang J.Consistent Estimators in Generalized Linear Mixed Models.J Am Stat Assoc,1998,93(442):720-729.

(責(zé)任編輯：鄧 妍)

*：國家自然科學(xué)基金項目(81402765，81473070，81373102)；國家統(tǒng)計局全國統(tǒng)計科學(xué)研究項目(2014LY112)；江蘇省教育廳高校哲學(xué)社會科學(xué)研究基金項目(2013SJD790032，2013SJB790059)；南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院優(yōu)秀博士論文培育項目

1. 徐州醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)教研室(221004)

2. 南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系