混合效应模型中的方差成分检验

・１６６・　中国卫生统计２０１７年２月第３４卷第１期　混合效应模型中的方差成分检验　曾　平　赵杨　陈峰。　

在很多科学问题中需要在混合效应模型框架下对　随机效应方差成分（暂记为　）进行假设检验　ｊ，也　即检验　：ｒ　＝０。除直接科学兴趣外，许多间接医学　问题也能转化为对方差成分的检验。例如，为判断在　惩罚样条回归中是参数模型还是非参数模型更合适，　Ｃｌａｅｓｋｅｎｓ　首先建立混合效应模型，将模型选择问题　转化为对随机效应方差成分是否等于零的假设检验问　题，最后通过限制性似然比检验　：丁　＝０。其他研究　者也用同样的方法处理过类似问题Ⅲ８　。然而，方差　成分为非负参数，对方差成分的假设检验是非标准的：　在　下　位于参数空间边缘。由于这种限制，常用　的渐近　无效分布对似然比统计量不再成立¨　’　。　混合效应模型中的方差成分检验吸引了广泛研究兴　趣ｌ卜８．　。本文主要介绍方差成分检验的似然比方　法，方差成分得分检验或Ｗａｌｄ检验可参考其他相关　文献　１３，２０，２１］。　线性混合效应模型的方差成分检验　设ｙ为连续变量，　＝（Ｘ　－－，　）和Ｚ＝（ｚ　，…，　Ｚｑ）分别为，ｚ×Ｐ和ｎ×ｑ的矩阵，　为样本量。在广　泛意义上】，、　和Ｚ既适于非独立数据结构　Ｉ２　，如　重复测量数据或聚群数据　，也适于独立数据情　形　¨，如惩罚样条回归。更加详细的例子可参　考上述文献。ｙ、　和ｚ之间的关系由线性混合效应　模型　刻画　ｐ　ｑ　Ｙ＝∑Ｘｊｆｌｊ＋∑Ｚｋｂ　＋　，　（１）　用矩阵的形式模型（１）可表达为　Ｙ＝　＋ｚ易＋ｓ，ｂ～Ⅳ（０，ｒ　Ｋ），占～Ⅳ（０，　Ｉｎ）　（２）　式中　为固定效应，ｂ为随机效应，　为未知方差成　分，　为已知ｑ×ｑ矩阵描述随机效应间相互关系，Ｌ　：国家自然科学基金项目（８１４０２７６５，８１４７３０７０，８１３７３１０２）；国家统计　局全国统计科学研究项目（２０１４ＬＹ１　１２）；江苏省教育厅高校哲学社会　科学研究基金项目（２０１３ＳＪＤ７９００３２，２０１３ＳＪＢ７９００５９）；南京医科大学公　共卫生学院优秀博士论文培育项目　１．徐州医科大学公共卫生学院流行病与卫生统计学教研室（２２１００４）　２．南京医科大学公共卫生学院生物统计学系　为ｎ维单位阵，　为残差方差，与ｂ不相关。模型中　ｙ的期望和方差分别为　Ｅ（ｙ）＝　，Ｖａｒ（１，）＝．１－２ＺＫＺ　＋　Ｉ　＝ｏｒ　＾　（３）　其中　＾＝ＡＺＫＺ　＋Ｊ　，Ａ＝，／２／ｏ－。。设参数为０＝（　，Ａ，　ｏｒ　），模型（２）的对数似然函数为　’　ｒ　１　（０）：一　１｛，ｚ　ｌｏｇ（ｏ＂　）＋　（ｙ一　）　（ｙ—　

Ｌ　）＋ｌｏｇｌ　１）　（４）　似然比统计量定义为　

＝２ｓｕｐ［　（０）一　（０）］　（５）　１．渐近混合卡方分布　显然检验／４０：ｂ＝０等价于检验／４０：　＝０。如前　文所述虽然定义似然比统计量并不难，但是　：　＝０　是非标准形式检验，渐近　无效分布不再成立。这　里，与非标准形式检验对应是标准形式检验，后者指参　数位于参数空间内部　。Ｓｅｌｆ和Ｌｉａｎｇ…、Ｌｉａｎｇ和　Ｓｅｌｆｌ３　以及Ｓｔｒａｍ和Ｌｅｅ＿２　假设：如ｙ能够划分为许多　独立同分布的亚向量、亚向量维度趋于无穷且研究为　均匀设计，那么职　渐近服从０．５Ｑ　＋０．５ｏｘ　，　为　退化到零的点分布，　是自由度为１的卡方分布。然　而，Ｃｒａｉｎｉｃｅａｎｕ和Ｒｕｐｐｅｒｔ　认为，Ｓｔｒａｍ和Ｌｅｅ的假设　条件在实际中很难满足，并显示０．５ｏｘ　＋０．５ｏｘ　会导　致保守结果。为验证这一点，图１给出了模拟例子。　对特殊纵向数据，Ｐｉｎｈｅｉｒｏ和Ｂａｔｅｓ　发现非等比例卡　方混合分布更加合适，如０．６５ｘ　＋０．３５ｘ　。但对更一　般的似然比方差成分检验，Ｃｒａｉｎｉｃｅａｎｕ和Ｒｕｐｐｅｒｔ＿８　证　明　和　的混合比例取决于具体数据，没有满足所　有情形的混合比例参数Ｈ　。　数据根据模型（１）模拟　，设丁　＝０，然后对比　Ｌ尺　和０．５ｏｘ　＋０．５ｏｘ　的分位数。０．５ｏｘ　＋０．５ｏｘ　的分位数总是大于对应　的分位数（散点和虚线　重叠表示两分布一致），意味着在０．５０Ｘｏ＋０．５ｏｘ　提　供的临界值过大。因此如果采用０．５　＋０．５　作为　似然比统计量的理论无效分布将产生保守结果，降低　统计效能和导致假阴性结果。　２．基于谱分解的模拟分布　为了获得似然比统计量的无效分布，Ｃｒａｉｎｉｃｅａｎｕ　Ｃｈｉｎｅｓｅ　Ｊｏｕｍａ１　ｏｆＨｅａｌｔｈ　Ｓｔａｔｉｓｔｉｃｓ．Ｆｅｂ．２０１７．Ｖｏ１．３４．Ｎｏ．１　模拟计算的似然比统计量　模拟计算的似然比统计量　模拟计算的似然比统计量模拟计算的似然比统计量　图１模拟计算的似然比统计量和０．５　＋０．５０ｘ　的分位数。　和Ｒｕｐｐｅｒｔ在谱分解基础上提出了一种基于模拟抽样　残差，由　条件下的模型计算获得；ｈ。　为矩阵Ｈ＝　算法　。假设模型（２）中Ａ已知，可获得卢和　的极ｘ（ｘ　Ｘ）　Ｘ　对角线上元素。在　：　＝０条件下，模型　大似　估计值　（１）退化为简单线性模型Ｙ：∑７　，ｑ－Ｅｏ假设台和　卢＝（　ｌＸ）　ｙ，　＝（ｙ　Ｐ　１ｐ￣ｙ）／　分别是ＪＢ和　ｚ在　模型下的　值。基于重抽样　

，　获得似然比统计量无效分布可执行为：①根据原始数　带人对数似然函数，获得剖面似然函数　据ｊ［）获得似然比统计量的观察值，记做职　如；②获得　

（Ａ）＝÷｛，ｚ　１ｏｇ（　ｌ，）＋ｌｏｇｆ　ｆ｝　（７）　重抽样数据　，采用相同的模型和估计算法获得关　于Ｄ　似然比统计量，记做职　；③重复第二步　次，　一ｆ　ＲＴ　ｒ　一　得到Ｂ个　；则似然比统计量统计量的无效分布由　

～　职　组成，其对应的Ｐ值由ＬＲｎ　７　等于或大于职　，　、　…，：ＬＬ，　，７　／’Ｊ　Ｈ　Ｌｔｉ　Ｌ　Ｊ　ｎ　、　Ｊ　’～／　Ｊ　…ｎ　

｛，ｚｌｏｇ［¨　ｇ　¨　））　比例　ｐ：　。具体而言在　是矩阵　ｚ　ＺＫ　的特征根，其中，　第二步中如果通过分布Ｎ（　，　）获得重抽样】，　、同　（Ａ）　２，Ｄｎ（Ａ）　丁＿　２＋　时保持　和Ｚ不变，则为参数ｂ。。ｔｓｔｒａｐ检验；如果Ｙ　

ｎ叩　＝Ｘ　＋台　、且台　通过有放回方式从台。随机获得，那　Ｔ，　，０、　厶＋１“　么为非参数ｂｏｏｔｓｔｒａｐ检验；如果Ｙ　＝　ｏ＋　、且　

是　ｚ　ＰｏＺＫ　特征根，Ｐ。＝，　一Ｘ（Ｘ　Ｘ）　Ｘ　。则　ＬＲＴ　以分布收敛于．厂ＬＫ　。　基于上述谱分解可建立模拟算法获得职　的无　效分布。需要注意的是，该分布是有限样本下的精确　分布而非大样本近似分布。谱分解模拟算法只与随机　效应维度有关、与样本量无关，因此在实际中具有较高　计算速度　。。。　。该算法获得的无效分布能够保证似　然比检验有效控制Ｉ型错误率，与其他方差成分检验　方法相比具有更高统计效能　，　。’　。当Ｓｔｒａｍ和　Ｌｅｅ假设条件满足时，基于模拟得到的精确分布弱收　敛于０．５０￣　＋０．５０￣　。　３．重抽样方法　统计分析中当面对复杂假设检验时诉诸于重抽样　方法是一种常用策略　。重抽样避免了数学推导，　概念上简单并且容易操作，属于计算密集型统计技术。　Ｆｉｔｚｍａｕｒｉｃｅ等　以及Ｌｅｅ和Ｂｒａｕｎ＿４　提出了似然比　方差成分的ｐｅｒｍｕｔａｔｉｏｎ检验，Ｆａｒａｗａｙ　建立了似然　比方差成分的参数ｂｏｏｔｓｔｒａｐ检验。其他研究者对似　然比方差成分提出了类似的重抽样检验方法¨　Ｊ。　

设　。　＝　￣／１一ｈ。　称为调整残差，称台。为原始　

通过无放回方式从台。随机获得，那么为ｐｅｒｍｕｔａｔｉｏｎ检　验。由此可见，参数ｂｏｏｔｓｔｒａｐ检验、非参数ｂｏｏｔｓｔｒａｐ　检验和ｐｅｒｍｕｔａｔｉｏｎ检验的主要差别在于如何产生新　的数据。这三种重抽样方法往往获得相似的结果，但　相对参数和非参数ｂｏｏｔｓｔｒａｐ检验，ｐｅｒｍｕｔａｔｉｏｎ检验假　设条件更少　，　。另外进行重抽样的对象是调整　残差　而非原始残差　，这是因为　。是同方差，而　是异方差　。一般而言，实际数据分析　取值１０００　左右就足够　’　．３　５１　。　

广义线性混合效应模型的方差成分检验　１．１ｏｇｉｓｔｉｃ混合效应模型和ＰＱＬ算法　目前几乎所有似然比方差成分检验都是在线性混　合效应模型框架下完成，即都是针对连续应变量的；而　对离散应变量似然比方差成分检验的研究还很少。以　ｌｏｇｉｓｔｉｃ回归为例，接下来讨论广义混合效应模型的似　然比方差成分检验。设ｙ为二分类应变量，　和ｚ同　前。ｙ、　和ｚ之间的关系通过ｌｏｇｉｓｔｉｃ混合效应模　型　卜　描述　

（】，），ｌｏｇ【　）＝　＋　＋Ｚｂ，　・１６８・　ｂ—Ｎ（Ｏ，Ｔ２Ｋ）　（１０）　为ｙ期望。感兴趣的检验是１４ｏ：西＝０，同样等价于　Ｈｎ：丁　＝０。模型（１Ｏ）的对数似然函数可表示为：Ｌ　ｎ　（／３，丁　）＝∑Ｊ．Ｅｆ（Ｙｉ　ｌ／３，ｂ）ｆ（ｂ　ｌ丁　）］　。似然比检验　

ｉ＝ｌ　统计量　

２［Ｌ　。（Ｊ８，　）一Ｌ（／３，　＝０）］　（１１）　然而，由于对数似然函数涉及积分运算，实际中除极有　限和简单情况外，大多数时候都无法精确计算　Ｌ　。（卢，　）和获得精确的对数似然值。因此不能直　接通过式（１１）进行似然比检验。　ｌｏｇｉｓｔｉｃ混合模型的参数估计算法包括数值积　分　，惩罚拟似然（ｐｅｎａｌｉｚｅｄ　ｑｕａｓｉ—ｌｉｋｅｌｉｈｏｏｄ，ＰＱＬ）和　边际拟似然（ｍａｒｇｉｎａｌ　ｑｕａｓｉ—ｌｉｋｅｌｉｈｏｏｄ，ＭＱＬ）　，伪似　然估计　川，Ｍｏｎｔｅ　Ｃａｒｌｏ估计　和Ｇｉｂｂｓ抽样　。　数值积分只适用于低维度的积分运算，一般很难处理　超过五维的积分　；Ｍｏｎｔｅ　Ｃａｒｌｏ估计和Ｇｉｂｂｓ抽样的　不足在于计算繁重。ＰＱＬ和ＭＱＬ能够通过近似方法　有效处理高维积分问题，计算相对简单，且能够通过重　复调用已有的线性混合模型程序执行。对ｌｏｇｉｓｔｉｃ混　合模型进行方差成分的似然比检验至少包含以下的三　个困难　］：与线性混合模型不同，通常不能获得关于　ｌｏｇｉｓｔｉｃ混合模型对数函数的闭型表达式，因此精确计　算ｌｏｇｉｓｔｉｃ混合模型的对数似然值和似然比统计量变　得很困难；由于涉及高维计算运算，很难获得ｌｏｇｉｓｔｉｃ　混合模型参数的精确估计值；即使能够有效估计ｌｏｇｉｓ—　ｃｉｃ混合模型和计算其对数似然比统计量，如何获得该　统计量的无效分布也不清楚。基于ｌｏｇｉｓｔｉｃ混合模型　已有的理论发展，本文仅仅专注其方差成分的似然比　检验。　２．基于重抽样方法的伪似然比方差成分检验　为了处理上面的问题，可以通过基于重抽样的伪　似然比方法进行方差成分检验　。具体步骤如下：①　采用ＰＱＬ算法估计ｌｏｇｉｓｔｉｃ混合模型参数；ＰＱＬ算法　可通过Ｒ软件中的函数ｇｌｍｍＰＱＬ进行，该函数包含　在ＭＡＳＳ软件包中　；②当ＰＱＬ算法收敛时，有工作　应变量　Ｙ　＝　＋Ｚ西＋８，ｂ～Ｎ（０，Ｔ２Ｋ），　—Ｎ（０，Ｒ）　（１２）　和伪似然函数　１　．　．　（　）＝一　一［１ｏｇ　ｌ　Ｄ　Ｉ＋（ｙ　一Ｘ卢）　Ｄ一　（ｙ　一Ｘ　）］　厶　＋Ｃ　（１３）　ｃ为一个常数，Ｄ＝ｒ　ＺＫＺ　＋Ｒ，　＝（　Ｄ　）　Ｄ　ｙ　。③将　（．ｒ　）看做ｙ　的对数似然函数，计算伪　似然比统计量ＬＲＴ　＝２ｓｕｐ　ｏ［　ｖＨ！（７－　）一　，Ｈ０　（　）］；④通过重抽样的方法获得伪似然比统计量　ＬＲＴ　的无效分布。　中国卫生统计２０１７年２月第３４卷第１期　其他问题和可能的研究方向　目前，由于以下几方面的原因，似然比方差成分检　验主要集中在线性混合效应模型并且模型中只包含一　个方差成分　’他’＂　：①能够较为容易估计线性混合　效应模型的参数和计算其对数似然值，进而计算似然　比统计量；②多个方差成分共存时，不存在简单的对数　似然函数谱分解形式，因此很难获得似然比统计量的　无效分布；③由于可能的高维积分运算，一般难以精确　估计广义混合效应模型的参数和计算其对数似然值，　以及获得似然比统计量的无效分布；④重抽样方法属　于计算密集型统计方法，难以大规模运用。　在标准的参数假设检验中似然比检验被证明是在　小样本时最优的　。事实上，由于似然比检验充分利　用了Ｈ。条件下的模型和Ｈ　条件下的模型，相对得分　检验和Ｗａｌｄ检验，似然比检验被证明在非标准参数　假设检验也具有更高统计效能　．３　Ｕ＇如Ｊ。因此，推广和　发展似然比检验以适合更加广泛的情形具有理论和实　际意义。具体而言，以下几方面问题值得进一步探索：　①在线性混合效应模型中当存在多个方差成分而仅对　其中某个成分感兴趣时，Ｇｒｅｖｅｎ　等建议可以首先从　应变量中消除冗余随机效应；此时，尚包括应变量的残　差和待检验方差成分，然后按照前文只包含单个方差　成分的似然比检验进行统计分析。然而，这样做的缺　点在于，没有考虑到随机效应之间可能存在的相关，而　且在模型估计多个方差成分时还可能存在数值问题，　如模型不收敛或计算不稳定。②在线性混合效应模型　中假设残差是独立的；当存在相关协方差结构，即　～　Ｎ（０，　（　）），　为一般形式协方差矩阵、包含未知参　数向量　；在这种情况下可先估计残差协方差，在此基　础上建立伪数据，然后将伪数据当做原始数据按照线　性混合效应模型框架下相似的方法进行似然比方差成　分检验　。但这种方法在协方差矩阵　的选择以及　其适用性等方面需要更多的研究。③前文在广义线性　混合效应模型框架下建立的伪似然比方差成分检验其　有效性取决于ＰＱＬ算法；ＰＱＬ算法本身属于有偏估　计；虽然提出了校正ＰＱＬ偏倚的算法　卜　Ｊ，有偏和无