唐明生黃水平金英良張耀東趙 楊陳 峰曾 平△

重抽樣方差成分檢驗(yàn)的多位點(diǎn)關(guān)聯(lián)分析*

唐明生1黃水平1金英良1張耀東1趙 楊2陳 峰2曾 平1△

目的在線性混合模型框架下探索基于重抽樣方差成分的似然比檢驗(yàn)多位點(diǎn)關(guān)聯(lián)分析方法。方法假設(shè)SNP位點(diǎn)效應(yīng)服從正態(tài)分布，將多位點(diǎn)關(guān)聯(lián)分析轉(zhuǎn)化為對(duì)隨機(jī)效應(yīng)方差的似然比檢驗(yàn)，采用重抽樣方法獲得統(tǒng)計(jì)量的零分布，通過混合分布提高計(jì)算速度。結(jié)果模擬研究顯示重抽樣檢驗(yàn)以及混合分布表現(xiàn)良好，能夠有效控制I型錯(cuò)誤，統(tǒng)計(jì)效能優(yōu)于現(xiàn)有方法。結(jié)論似然比檢驗(yàn)是一種高效能的多位點(diǎn)關(guān)聯(lián)檢驗(yàn)方法，通過對(duì)置換和bootstrap分布的有效近似，基于重抽樣的似然比方差成分檢驗(yàn)，在保持良好統(tǒng)計(jì)效能的同時(shí)，能明顯降低計(jì)算負(fù)荷。

關(guān)聯(lián)研究 方差成分 線性混合模型 似然比檢驗(yàn) 置換法 非參數(shù)bootstrap法

在全基因組關(guān)聯(lián)性研究（genome w ide association study，GWAS）中，相對(duì)于單位點(diǎn)分析，多位點(diǎn)分析已經(jīng)成為一種有效的關(guān)聯(lián)性檢驗(yàn)方法［1－8］。單位點(diǎn)分析對(duì)成千上萬個(gè)單核苷酸多態(tài)性（single nucleotide polymorphisms，SNP）依次檢驗(yàn)，要求十分嚴(yán)格的多重檢驗(yàn)水準(zhǔn)［9］。相反，多位點(diǎn)分析同時(shí)檢驗(yàn)一組SNP是否與表型相關(guān)，有諸多優(yōu)點(diǎn)：能顯著降低多重比較負(fù)擔(dān)，能利用SNP間連鎖不平衡信息來提高效能。考慮到位點(diǎn)間聯(lián)合效應(yīng)，SNP集合常根據(jù)基因選定，也即位于相同基因內(nèi)的變異位點(diǎn)組成一個(gè)集合［10］。

多位點(diǎn)分析可采用固定效應(yīng)模型和固定效應(yīng)檢驗(yàn)執(zhí)行，如Wald檢驗(yàn)；但固定效應(yīng)檢驗(yàn)在SNP位點(diǎn)個(gè)數(shù)較多時(shí)效能低下。另一方面，在線性混合模型下假設(shè)SNP具有隨機(jī)效應(yīng)，能夠?qū)⒍辔稽c(diǎn)分析轉(zhuǎn)化為對(duì)隨機(jī)效應(yīng)的方差成分檢驗(yàn)；現(xiàn)階段主要通過得分檢驗(yàn)來實(shí)現(xiàn)這一目的［1，8，11－16］。得分檢驗(yàn)只需擬合零模型，計(jì)算效率高；但結(jié)果相對(duì)保守，特別是在小樣本的情況下［16］。本文嘗試采用似然比檢驗(yàn)（likelihood ratio test，LRT）及限制性似然比檢驗(yàn)（restricted likelihood ratio test，ReLRT）對(duì)SNP集合進(jìn)行關(guān)聯(lián)性檢驗(yàn)［17－19］，通過置換及bootstrap［20－23］等重抽樣方法獲得統(tǒng)計(jì)量的零分布。事實(shí)上，在關(guān)聯(lián)性研究中為了避免導(dǎo)出解析解，重抽樣技術(shù)已經(jīng)廣泛應(yīng)用于復(fù)雜假設(shè)檢驗(yàn)［6，24－29］。然而，重抽樣屬于計(jì)算密集型統(tǒng)計(jì)方法，為了減少計(jì)算負(fù)荷，本文進(jìn)一步通過混合分布來近似重抽樣分布，以提高似然比檢驗(yàn)的計(jì)算速度。最后通過模擬和實(shí)際數(shù)據(jù)來展示本文的方法。

資料與方法

1．線性混合模型和似然的推斷

用線性混合模型［30－32］來描述多個(gè)SNP位點(diǎn)和表型之間的關(guān)系：

Y是n維連續(xù)性表型，X為n×m的協(xié)變量矩陣（m為協(xié)變量個(gè)數(shù)），具有固定效應(yīng)α，G表示n×p的SNP基因型矩陣（編碼0，1，2，p為SNP個(gè)數(shù)），具有隨機(jī)效應(yīng)b，服從均數(shù)為0、方差為τ的正態(tài)分布。ε為殘差項(xiàng)，服從均數(shù)為0、方差為σ2的正態(tài)分布。

檢驗(yàn)G與表型Y之間的關(guān)系等價(jià)于檢驗(yàn)隨機(jī)效應(yīng)H0：b＝0，也等價(jià)于檢驗(yàn)方差成分H0：τ＝0。這樣，通過混合效應(yīng)模型SNP集合的關(guān)聯(lián)分析就轉(zhuǎn)化為對(duì)隨機(jī)效應(yīng)方差成分的假設(shè)檢驗(yàn)。

省略常數(shù)項(xiàng)，模型（1）的對(duì)數(shù)似然函數(shù)和限制性對(duì)數(shù)似然函數(shù)分別為：

與得分檢驗(yàn)不同，似然比檢驗(yàn)和限制性似然比要求同時(shí)擬合零模型和備擇模型，本文通過R軟件的lm和lme軟件包［33－34］來估計(jì)線性模型和混合效應(yīng)模型。

在H0的條件下，τ位于參數(shù)空間的邊緣；由于這種限制，TLRT和TReLRT的零分布不再漸近服從標(biāo)準(zhǔn)的卡方分布［35］，先前的研究顯示在一定條件下，TLRT和TReLRT漸近服從的混合分布［18］。但是，該混合分布的假定條件在實(shí)際中很難得到滿足。Pinheiro和Bates提出了0．65：0．35的混合分布［36］。本文模擬研究表明，在多位點(diǎn)似然比關(guān)聯(lián)性分析中兩種混合比例都不正確。事實(shí)上，混合分布比例參數(shù)依具體情況而定，依賴于特殊矩陣的特征值［37－39］。本文通過重抽樣方法來獲得似然比統(tǒng)計(jì)量和限制性似然比統(tǒng)計(jì)量的零分布。

2．重抽樣算法

當(dāng)烘絲機(jī)入料電子秤不再檢測(cè)到煙絲，烘絲機(jī)即由生產(chǎn)狀態(tài)切換至收尾狀態(tài)，筒壁Ⅰ區(qū)、Ⅱ區(qū)溫度控制器停止工作，Ⅰ區(qū)、Ⅱ區(qū)蒸汽壓力控制閥保持其最后適用的修正變量60 s后；Ⅰ區(qū)、Ⅱ區(qū)蒸汽壓力控制閥開始逐漸關(guān)閉Ⅰ區(qū)、Ⅱ區(qū)蒸汽閥門開度，在20 s內(nèi)完全關(guān)閉Ⅰ區(qū)、Ⅱ區(qū)，蒸汽停止加熱筒壁。在此過程中，因?yàn)楹娼z機(jī)蒸汽壓力控制閥對(duì)筒壁Ⅰ區(qū)、Ⅱ區(qū)蒸汽壓力保持以及對(duì)Ⅰ區(qū)、Ⅱ區(qū)蒸汽的關(guān)閉均為同時(shí)進(jìn)行，且烘絲機(jī)由收尾狀態(tài)至蒸汽壓力保持的延時(shí)時(shí)間以及對(duì)筒壁Ⅰ區(qū)和Ⅱ區(qū)蒸汽壓力的保持時(shí)間較長(zhǎng)，使烘絲機(jī)熱慣性過大，導(dǎo)致了收尾狀態(tài)產(chǎn)生大量的干尾煙絲。

置換和bootstrap法都是通過重復(fù)抽樣獲得零分布，不同之處在于兩者的抽樣方式。置換法通過擾亂原始數(shù)據(jù)的標(biāo)簽產(chǎn)生，需要將Y和X視為一個(gè)整體，即置換是對(duì)Y和X同時(shí)進(jìn)行，而保持G固定不變；這樣可以保持Y與X的相關(guān)結(jié)構(gòu)。這里潛在假設(shè)G和X是獨(dú)立的。bootstrap法采用非參數(shù)有放回抽樣產(chǎn)生，即是從中有放回抽得的樣本。以上兩種產(chǎn)生數(shù)據(jù)D*的方法是等價(jià)的［21－22］。在bootstrap算法中，最自然的抽樣是從原始?xì)埐頧中進(jìn)行抽樣；然而Davison和Hinkley［21］指出，當(dāng)是異方差性時(shí)，從修正殘差中 抽樣更好，h為杠桿值。置換和bootstrap檢驗(yàn)的P值采用蒙特卡洛方法獲得，即通過有限次數(shù)（設(shè)為S）的重抽樣方式代替。很多研究表明，在檢驗(yàn)水準(zhǔn)較大時(shí)（如0．05），S＝1000或2000時(shí)P值就趨于相對(duì)穩(wěn)定［20－22］。

3．近似分布

假設(shè)TLRT和TReLRT服從混合分布

κ是未知的尺度參數(shù)分別為自由度為0和1的卡方分布。假設(shè)近似零分布為卡方混合的原因在于：（I）在適當(dāng)?shù)臈l件下，TLRT和TReLRT取0的概率不為零，因此包括作為退化到0的分布；（II）在標(biāo)準(zhǔn)條件下，即參數(shù)在其參數(shù)空間的內(nèi)部時(shí)，TLRT和TReLRT漸近服從自由度為1的卡方分布［40］，因此包括表示可能的偏離；類似的混合分布也被應(yīng)用在其他地方［17，36］。κ的估計(jì)值為：

本文采用局部概率法來估計(jì)π［38－39］：

μ是G′P0G的特征值，P0＝In－X（X′X）－1X′，ξ是G′G的特征值，u是獨(dú)立標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量。局部概率法相對(duì)于矩估計(jì)法的優(yōu)勢(shì)在于：其π估計(jì)值與L的選擇無關(guān)，因此π估計(jì)值是穩(wěn)定的，并會(huì)導(dǎo)致更加精確的估計(jì)。對(duì)（6）式可采用Davies法［42］獲得π的估計(jì)值然后將代入（5）中可算得式（6）還進(jìn)一步表明：TLRT和TReLRT的零分布依具體情況而定，固定比例的卡方混合分布是不恰當(dāng)?shù)摹?/p>

結(jié) 果

用FEV1下降率的數(shù)據(jù)［43］來展示本文提出的LRT和ReLRT多位點(diǎn)關(guān)聯(lián)性分析方法，樣本量為301，用FEV1下降的斜率ES作為表型［43］。結(jié)果表明，位點(diǎn)rs9469089與FEV1相關(guān)聯(lián)；該位點(diǎn)位于染色體6p21．32，在基因RNF5的第一個(gè)內(nèi)含子區(qū)域內(nèi)，而RNF5編碼膜結(jié)合型泛素連接酶。在該數(shù)據(jù)中RNF5一共包含14個(gè)SNP。

首先評(píng)價(jià)I類錯(cuò)誤控制。假設(shè)表型來自于

X1和X2分別為標(biāo)準(zhǔn)正態(tài)變量和二分類變量，模擬次數(shù)為105次。重抽樣中S＝2000，對(duì)于近似混合分布，取L＝2000、1500、1000、800、500、300和100。用m lp1－m lp7對(duì)應(yīng)以上各L取值。在評(píng)價(jià)統(tǒng)計(jì)效能時(shí)，假設(shè)表型來自于

這里，j取1到14，即每個(gè)SNP位點(diǎn)依次被設(shè)為關(guān)聯(lián)位點(diǎn)，效應(yīng)值為0．20，0．30和0．40，重復(fù)103次，則一共運(yùn)行1．4×104次。P值用小于檢驗(yàn)水準(zhǔn)α＝0．05的比例估計(jì)。

2．解釋和說明

模擬結(jié)果顯示，在LRT和ReLRT中κ與π的平均值分別為1．118和0．909與0．865和0．604；這些數(shù)值和本文特定的基因型矩陣和協(xié)變量矩陣相關(guān)，其他不同的選擇會(huì)導(dǎo)致不同的κ和π估計(jì)值。這表明Self和Liang［17］及Pinheiro和Bates［36］提出的零分布是不合理的，同時(shí)也表明LRT和ReLRT服從不同的零分布。圖1表明參數(shù)π在有限的范圍內(nèi)變化（這是因?yàn)樵谒械哪M中使用相同的G和相似的X），將其固定取某一常數(shù)是不適當(dāng)?shù)?。根?jù)公式（6），π隨數(shù)據(jù)集而變化，因此直接根據(jù)數(shù)據(jù)估計(jì)π更合理。此外，模擬還表明，相對(duì)于在ReLRT中，κ估計(jì)值在LRT中更大（圖1）。κ估計(jì)值的精度隨著L減小而降低；然而，不同L值通常產(chǎn)生非常相似的結(jié)果。

圖2表明對(duì)LRT而言，模擬算法［39，45］估計(jì)的I類錯(cuò)誤率結(jié)果偏于保守，而置換和bootstrap檢驗(yàn)可以調(diào)整LRT的這種保守性；混合分布對(duì)I類錯(cuò)誤的控制與L取值無關(guān)。ReLRT在所有情形下對(duì)I類錯(cuò)誤控制都表現(xiàn)良好；這可能是因?yàn)長(zhǎng)RT對(duì)方差成分的估計(jì)是有偏估計(jì)，而ReLRT對(duì)方差成分的估計(jì)是無偏估計(jì)［46－47］。LRT模擬算法的統(tǒng)計(jì)效能低于相應(yīng)置換和bootstrap檢驗(yàn)的統(tǒng)計(jì)效能；在ReLRT中，模擬算法、置換和bootstrap檢驗(yàn)的統(tǒng)計(jì)效差異微小（圖3）。此外，通常次要等位基因頻率高，統(tǒng)計(jì)效能也較高（圖3）。

圖1 在不同的L值下，LRT和ReLRT近似混合分布的尺度參數(shù)κ和混合比例參數(shù)π的估計(jì)值

圖2 LRT和ReLRT I類錯(cuò)誤率的置換或bootstrap算法估計(jì)

LRT和ReLRT的統(tǒng)計(jì)效能通常高于得分檢驗(yàn)的統(tǒng)計(jì)效能。對(duì)于得分檢驗(yàn)、LRT和ReLRT，在模擬效應(yīng)為0．40時(shí)三者的平均統(tǒng)計(jì)效能分別為0．474、0．544和0．545；模擬效應(yīng)為0．30時(shí)為0．315，0．453和0．510；模擬效應(yīng)為0．20時(shí)為0．179，0．232和0．250。圖3表明：對(duì)于LRT（或ReLRT），置換和bootstrap兩種檢驗(yàn)方法的統(tǒng)計(jì)效能相似。

在FEV1數(shù)據(jù)中將斜率（即ES）作為表型［43］，內(nèi)毒素暴露和BMI作為協(xié)變量，采用LRT和ReLRT以及得分檢驗(yàn)來分析基因RNF5的關(guān)聯(lián)性（表1）。從表1可見，在LRT中由模擬算法獲得的P值比其他情形下的P值大，置換和bootstrap檢驗(yàn)的P值更小；ReLRT的P值變化較小。這些結(jié)果和模擬的結(jié)論一致。得分檢驗(yàn)的P值為0．027。

圖3 LRT和ReLRT統(tǒng)計(jì)效能的置換或bootstrap算法估計(jì)

表1 FEV1數(shù)據(jù)中基因RNF5的P值

討 論

本文在線性混合模型框架下研究了基于重抽樣的似然比多位點(diǎn)關(guān)聯(lián)分析，在該框架下，一組SNP位點(diǎn)與表型的關(guān)聯(lián)性分析被轉(zhuǎn)化為對(duì)隨機(jī)效應(yīng)的方差成分檢驗(yàn)［1，11－16］。本文采用重抽樣技術(shù)（置換和bootstrap）獲得似然比統(tǒng)計(jì)量的零分布，避免了復(fù)雜的數(shù)學(xué)推導(dǎo)。此外，還通過混合分布近似置換或bootstrap分布。本文采用非參數(shù)而不是參數(shù)的bootstrap法，其原因在于非參數(shù)bootstrap法更加穩(wěn)健［20－21］。模擬表明，對(duì)于限制似然比檢驗(yàn)，置換和bootstrap法能有效控制I類錯(cuò)誤且統(tǒng)計(jì)效能高于現(xiàn)有的得分檢驗(yàn)方法；然而，對(duì)于似然比檢驗(yàn)，其I類錯(cuò)誤未能得到有效控制。研究還發(fā)現(xiàn)，對(duì)于小樣本基于模擬算法的似然比檢驗(yàn)［39，41］會(huì)導(dǎo)致保守的結(jié)果，而隨著樣本量的增加，似然比檢驗(yàn)表現(xiàn)趨于正常［45］。

重抽樣方法的主要缺點(diǎn)是耗時(shí)，本文利用混合分布近似來減少計(jì)算負(fù)荷。模擬表明，近似分布能顯著減少計(jì)算時(shí)間。例如，L＝100時(shí)的混合分布比S＝2000時(shí)的重抽樣檢驗(yàn)計(jì)算時(shí)間大約減少20倍。更重要的是，在大規(guī)模多位點(diǎn)關(guān)聯(lián)分析時(shí)，要求更加精確的P值，這需要成千上萬次的置換或bootstrap抽樣，導(dǎo)致重抽樣方法不可行；而近似混合分布能夠適合這種情形。如前所述，在混合分布的參數(shù)估計(jì)中，和其他方法相比，局部概率法具有數(shù)值穩(wěn)定的優(yōu)點(diǎn)。

類似的重抽樣方法也被應(yīng)用在其他情形，例如，F(xiàn)araway［48］和Samuh等［49］提出了線性混合模型似然比檢驗(yàn)的參數(shù)bootstrap法，F(xiàn)itzmaurice等［50］，Lee和Braun［51］，Samuh［49］提出了似然比檢驗(yàn)置換法，Sinha［52］在廣義線性混合模型下通過參數(shù)bootstrap法建立了方差成分的得分檢驗(yàn)。本文的研究與上述文獻(xiàn)的區(qū)別在于：（I）上述研究主要針對(duì)的是縱向數(shù)據(jù)，且這些數(shù)據(jù)在各水平內(nèi)是相關(guān)的，而模型（1）事實(shí)上應(yīng)用在看上去獨(dú)立的基于總體人群設(shè)計(jì)的遺傳數(shù)據(jù)；（II）對(duì)于本文的似然比檢驗(yàn)，同時(shí)采用置換和bootstrap法，且采用相對(duì)穩(wěn)健的非參數(shù)bootstrap法而不是參數(shù)bootstrap法；（III）將基于重抽樣的方法與其他方法（基于模擬的算法）進(jìn)行比較［39］，結(jié)果表明重抽樣方法更有效，雖然計(jì)算負(fù)荷更重；（IV）采用了有效的近似方法來提高計(jì)算效率。最后，有待進(jìn)一步研究將其他方法（如，基于參數(shù)bootstrap的得分檢驗(yàn)［52］）運(yùn)用到本文的研究中并進(jìn)行比較。

［1］Wu MC，Kraft P，Epstein MP，et al．Powerful SNP-Set Analysis for Case-Control Genome-wide Association Studies．American Journal of Human Genetics，2010，86（6）：929-942．

［2］Cai T，Lin X，Carroll RJ．Identifying genetic marker sets associated w ith phenotypes via an efficient adaptive score test．Biostatistics，2012，13（4）：776-790．

［3］Maity A，Sullivan PF，Tzeng Ji．Multivariate Phenotype Association Analysis by Marker-Set Kernel Machine Regression．Genetic Epidem iology，2012，36（7）：686-695．

［4］Schifano ED，Epstein MP，Bielak LF，et al．SNP Set Association A-nalysis for Fam ilial Data．Genetic Epidem iology，2012，36（8）：797-810．

［5］Dai H，Zhao Y，Qian C，et al．Weighted SNP Set Analysis in Genome-W ide Association Study．PLoSONE，2013，8（9）：e75897．

［6］Huang YT，Lin X．Gene set analysis using variance component tests．BMC Bioinformatics，2013，14（1）：210．

［7］Jiao S，Hsu L，Bézieau S，et al．SBERIA：Set-Based Gene-Environment Interaction Test for Rare and Common Variants in Complex Diseases．Genetic Epidem iology，2013，37（5）：452-464．

［8］Wu MC，Maity A，Lee S，etal．KernelMachine SNP-Set Testing Under Multiple Candidate Kernels．Genetic Epidemiology，2013，37（3）：267-275．

［9］Balding D．A tutorial on statisticalmethods for population association studies．Nature Reviews．Genetics，2006，7（10）：781-791．

［10］Prescott NJ，Lehne B，Stone K，et al．Pooled Sequencing of 531 Genes in Inflammatory Bowel Disease Identifies an Associated Rare Variant in BTNL2 and Implicates Other Immune Related Genes．PLoSGenetics，2015，11（2）：e1004955．

［11］Liu D，Lin X，Ghosh D．Semiparametric Regression of Multidimensional Genetic Pathway Data：Least-Squares Kernel Machines and Linear M ixed Models．Biometrics，2007，63（4）：1079-1088．

［12］Tzeng J，Zhang D．Haplotype-based association analysis via variancecomponents score test．American Journal of Human Genetics，2007，81（5）：927-938．

［13］Kwee L，Liu D，Lin X，etal．A Powerful and Flexible Multilocus Association Test for Quantitative Traits．American Journal of Human Genetics，2008，82（2）：386-397．

［14］Liu D，Ghosh D，Lin X．Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernelmachine regression via logistic m ixed models．BMC Bioinformatics，2008，9（1）：292．

［15］Tzeng J，Zhang D，Chang SM，et al．Gene-trait sim ilarity regression formultimarker-based association analysis．Biometrics，2009，65（3）：822-832．

［16］Wu MC，Lee S，Cai T，etal．Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test．American Journal of Human Genetics，2011，89（1）：82-93．

［17］Self SG，Liang KY．Asymptotic Properties of Maximum Likelihood Estimators and Likelihood Ratio Tests under Nonstandard Conditions．JAM STAT ASSOC，1987，82（398）：605-610．

［18］Stram DO，Lee JW．Variance Components Testing in the Longitudinal M ixed Effects Model．Biometrics，1994，50（4）：1171-1177．

［19］Liang KY，Self SG．On the Asymptotic Behaviour of the Pseudolikelihood Ratio Test Statistic．JRoy Stat Soc，1996，58（4）：785-796．

［20］Efron B，Tibshirani R．An Introduction to the bootstrap．New York：Chapman and Hall，1993．

［21］Davison AC，Hinkley DV．bootstrap Methods and Their Application．Cambridge：Cambridge University Press，1997．

［22］Good P．Permutation，Parametric，and bootstrap Tests of Hypotheses．3 edition．New York：Springer，2005．

［23］Efron B．The Jackknife，the bootstrap and Other Resampling Plans．Philadelphia：Society for industrial and applied mathematics，1982．

［24］Han F，Pan W．A data-adaptive sum test for disease association w ith multiple common or rare variants．Human Heredity，2010，70：42-54．

［25］Lin D．An efficientMonte Carlo approach to assesssing statistical significance in genom ic studies．Bioinformatics，2005，21（6）：781-787．

［26］Madsen BE，Browning SR．A Groupw ise Association Test for Rare Mutations Using a Weighted Sum Statistic．PLoS Genetics，2009，5（2）：e1000384．

［27］Lee S，Emond MJ，Bamshad MJ，et al．Optimal Unified Approach for Rare-Variant Association Testing w ith Application to Small-Sample Case-ControlWhole-Exome Sequencing Studies．American Journal of Human Genetics，2012，91（2）：224-237．

［28］Ferkingstad E，Holden L，Sandve GK．Monte Carlo nullmodels for genom ic data．Statistical Science，2015，30（1）：59-71．

［29］Lin D，Tang Z．A General Framework for Detecting Disease Associations w ith Rare Variants in Sequencing Studies．American Journal of Human Genetics，2011，89（3）：354-367．

［30］Laird NM，Ware JH．Random-effects models for longitudinal data．Biometrics，1982，38（4）：963-974．

［31］Littel R，M illiken GA，Stroup WW，et al．SAS for mixed models．SAS Institute Inc．Cary，NC，2006．

［32］Verbeke G，MolenberghsG．Linearmixedmodels for longitudinal data．New York：Springer，2009．

［33］Pinheiro J，Bates D，DebRoy S，et al．nlme：Linear and Nonlinear M ixed Effects Models．R package version 3．1-113．2013．

［34］R Core Team．R：A language and environment for statistical computing．R Foundation for Statistical Computing，Vienna，Austria．URL http：//www．R-project．org/，2014．

［35］Molenberghs G，Verbeke G．Likelihood Ratio，Score，and Wald Tests in a Constrained Parameter Space．Am Stat，2007，61（1）：22-27．

［36］Pinheiro JC，BatesD．M ixed-EffectsModels in S and S-PLUS，2nd edition．New York：Springer，2009．

［37］Kuo BS．Asymptotics of ML estimator for regression models with a stochastic trend component．Economet Theor，1999，15（1）：24-49．

［38］Claeskens G．Restricted likelihood ratio lack-of-fit tests using mixed splinemodels．JR Stat Soc B，2004，66（4）：909-926．

［39］Crainiceanu CM，Ruppert D．Likelihood ratio tests in linear mixed modelswith one variance component．J Roy Stat Soc B，2004，66（1）：165-185．

［40］Lehmann EL，Romano JP．Testing statistical hypotheses，Third edition．New York：Springer，2006．

［41］Greven S，Crainiceanu CM，Küchenhoff H，et al．Restricted Likelihood Ratio Testing for Zero Variance Components in Linear M ixed Models．JComput Graph Stat，2008，17（4）：870-891．

［42］Davies RB．Algorithm AS 155：The Distribution of a Linear Combination of chi-2 Random Variables．JR Stat Soc C，1980，29（3）：323-333．

［43］Zhang R，Zhao Y，Chu M，et al．A Large Scale Gene-Centric Association Study of Lung Function in New ly-Hired Female Cotton Textile Workers w ith Endotoxin Exposure．PLoSONE，2013，8（3）：e59035．

［44］Barrett JC，F(xiàn)ry B，Maller J，et al．Haploview：analysis and visualization of LD and haplotype maps．Bioinformatics，2005，21（2）：263-265．

［45］Zeng P，Zhao Y，Liu J，et al．Likelihood Ratio Tests in Rare Variant Detection for Continuous Phenotypes．Annals of Human Genetics，2014，78（5）：320-332．

［46］Patterson HD，Thompson R．Recovery of interblock information when block sizes are unqual．Biometrika，1971，58（3）：545-555．

［47］Harville DA．Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems．JR Stat Soc B，1977，72（358）：320-338．

［48］Faraway JJ．Extending the linear model w ith R：generalized linear，m ixed effects and nonparametric regressionmodels．New York：Chapman＆Hall/CRC，2005．

［49］Samuh MH，Grilli L，Rampichini C，et al．The Use of Permutation Tests for Variance Components in Linear M ixed Models．Commun Stat-Theor M，2012，41（16-17）：3020-3029．

［50］Fitzmaurice GM，Lipsitz SR，Ibrahim JG．A Note on Permutation Tests for Variance Components in Multilevel Generalized Linear M ixed Models．Biometrics，2007，63（3）：942-946．

［51］Lee OE，Braun TM．Permutation Tests for Random Effects in Linear M ixed Models．Biometrics，2012，68（2）：486-493．

［52］Sinha SK．Bootstrap tests for variance components in generalized linearm ixed models．Can JStat，2009，37（2）：219-234．

（責(zé)任編輯：郭海強(qiáng)）

國(guó)家自然科學(xué)基金項(xiàng)目（81402765）；國(guó)家統(tǒng)計(jì)局統(tǒng)計(jì)科學(xué)研究項(xiàng)目（2014LY112）

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

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

△通信作者：曾平，E-mail：zpstat＠xzhmu．edu．cn