混合效應模型中的方差成分檢驗*

曾 平 趙 楊 陳 峰

混合效應模型中的方差成分檢驗*

曾 平1趙 楊2陳 峰2

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

線性混合效應模型的方差成分檢驗

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

(1)

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

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

(2)

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

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

(3)

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

(4)

似然比統計量定義為

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

(5)

1.漸近混合卡方分布

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

為了獲得似然比統計量的無效分布，Crainiceanu

圖1 模擬計算的似然比統計量和的分位數。

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

(6)

帶入對數似然函數，獲得剖面似然函數

(7)

定義

(8)

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

(9)

3.重抽樣方法

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

廣義線性混合效應模型的方差成分檢驗

1.logistic混合效應模型和PQL算法

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

(10)

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

(11)

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

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

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

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

(12)

和偽似然函數

(13)

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

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

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

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

[28]曾平,曹紅艷,劉桂芬.重復測量計數資料的隨機效應ZIP模型.中國衛生統計,2009,26(4):344-346,349.

[29]任仕泉,陳峰.多水平統計模型在多中心臨床試驗評價中的應用.中國衛生統計,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 在假設檢驗中的應用.數理統計與管理,2006,25(5):616-621.

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

[39]劉沛,劉元寶,王燦楠.膳食暴露評估模型及其構建方法.中華預防醫學雜志,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]曹紅艷,劉桂芬,曾平,等.預測性偽似然法和貝葉斯法廣義線性混合模型估計.中國藥物與臨床,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.

(責任編輯：鄧 妍)

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

1. 徐州醫科大學公共衛生學院流行病與衛生統計學教研室(221004)

2. 南京醫科大學公共衛生學院生物統計學系