懲罰廣義估計方程在縱向數據基因關聯分析中的應用*

曹紅艷 曾 平 李 治 崔躍華，4張巖波△

懲罰廣義估計方程在縱向數據基因關聯分析中的應用*

曹紅艷1曾 平2李 治3崔躍華1，4張巖波1△

目的 探討懲罰廣義估計方程（pGEE）在縱向數據基因關聯分析的應用，為縱向數據基因關聯分析提供方法學參考。方法 以小鼠糖尿病發病相關的數量性狀位點識別為例，分別采用廣義估計方程（GEE）和pGEE進行分析。結果 pGEE篩選出糖尿病發病關聯位點，為分子生物學研究提供了重要的候選位點。結論 pGEE能有效的實現高維縱向數據的變量選擇，識別出有意義的關聯位點。

懲罰廣義估計方程 縱向數據 SCAD 基因關聯分析 數量性狀位點

隨著分子生物學測序技術的不斷發展，基因關聯分析已經成為了復雜疾病研究的最重要手段之一，成功的識別了多種人類復雜疾?。ㄈ纰蛐吞悄虿?，高血壓等）的遺傳變異位點［1-2］。為了更好地了解復雜疾病和遺傳變異之間的關系，越來越多的基因關聯分析通過縱向的方式進行，對同一觀察對象的某觀察指標在不同時間上進行重復測量，數據形式表現為非獨立結構。縱向數據基因關聯分析相比于橫斷面研究，可以研究復雜性狀隨時間變化的關系，從而增強統計效能，提高遺傳變異對復雜疾病的解釋程度［3-5］。

縱向基因數據和一般的縱向數據相比，更為復雜，高維度，不同SNP之間存在復雜的LD結構，表現為強相關性和多維共線性；同時，大部分的SNP為冗余信息，只有少量的SNP為關聯位點且信號強度弱，將導致參數估計和統計推斷的準確性和有效性大大降低，高維縱向基因關聯分析面臨著巨大的挑戰。因此，發展新的變量選擇方法尤為重要。近年來，基于懲罰的變量選擇方法備受關注，如Lasso［6］，SCAD［7］，自適應Lasso［8］等。懲罰方法能有效地應用于高維數據分析，通過收縮將弱效應估計為0，參數估計和變量選擇同時進行，廣泛的應用于不同模型的變量選擇。經典的縱向數據分析中，廣義估計方程（generalized estimating equations，GEE）［9-11］只需定義偽得分方程（quasi -score equations），當作業相關矩陣設置不正確時，仍然能得到一致性的參數估計，在應用中獨占優勢。因此，Wang，Zhou和Qu發展了基于SCAD的懲罰廣義估計方程（penalized generalized estimating equations，pGEE）［12］，不僅保持了GEE的重要特性，同時將GEE推廣到高維數據分析，適用于協變量個數p隨樣本例數n同階變化的情況，即P＝O（n）。本文以小鼠的糖尿病發病相關的數量性狀位點（quantitative trait locus，QTL）識別為例，進行pGEE分析，識別出小鼠的糖尿病發病相關的QTL，為QTL分析提供方法支持。

資料與方法

1．資料來源

數據來源于2型糖尿病發病相關的QTL遺傳研究［13］，由于肥胖是2型糖尿病主要危險因素，Reifsnyder等將肥胖且有糖尿病傾向的NZO（new zealand obese）／HILt小鼠與瘦且無糖尿病傾向的NON（nonobese nondiabetic）／Lt小鼠進行遠交（outcross），產生的F1代再與親本NON／Lt小鼠回交（backcross）生成多基因實驗小鼠模型。對203例雄性實驗小鼠，測量其在4、8、12、16、20、24周的體重，同時，在已知QTL位點周圍20 cM范圍內的高頻標記中，挑選了83個微衛星位點，研究微衛星位點對體重的影響，數據形式如表1所示。

表1 體重和微衛星標記數據集形式

2．方法

（1）廣義估計方程

在縱向數據基因關聯分析中，設Yit為個體i在時間t的觀察值，Xit為p維協變量向量，令Yi＝（Yi1，…Yit）T為個體i的觀察向量，Xi＝（Xi1，…，Xit）T為個體i的協變量矩陣，則縱向數據分析的模型框架可表達如下：

E（Yit）＝μi，g（μi）＝X′iβ

基于全似然函數的參數估計方法，只有在一些特殊情況下，如當Yit服從多元正態分布時，才能直接算出，對于其他的分布，全似然函數的計算非常復雜。GEE采用類似于廣義線性模型的得分方程的方法，定義偽得分方程為：

（2）懲罰廣義估計方程

基于懲罰的變量選擇方法有很多，Fan和Li指出一個好的懲罰函數估計值應具備無偏性、稀疏性和連續性，即Oracle性質［7］。SCAD懲罰能保留較大的系數，同時將較小的系數收縮為0，具有Oracle性質，其懲罰函數pλ（θ）的導數如下：

其中，I為指示函數，a為預先選擇的常數，Fan和Li推薦a＝3．7。

對GEE的得分函數進行SCAD懲罰，得到pGEE的懲罰表達式為：

采用minorization-maximization算法得到（2）式的漸進表達式，

其中Sj（β＾）為S（β＾）的第j個得分函數，ε的取值較小，一般取ε＝10-6。進一步采用New ton-Raphson算法對（3）進行迭代，計算得到pGEE的參數估計值。pGEE估計值依賴于懲罰參數λ，在此，采用交叉驗證（cross validation，CV）選擇最優λ。

（3）統計分析采用R軟件中的“PGEE 1．4”軟件包，可實現GEE和pGEE分析，GEE分析也可采用“geepack”包?！癙GEE 1．4”包含3個函數：MGEE，CVfit和PGEE?！癕GEE”函數用于GEE分析，“CVfit”用于計算交叉驗證的最優懲罰參數λ，得到最優λ后，采用“PGEE”函數計算pGEE的參數估計值。適用于測量結果呈正態分布、二項分布、Poisson分布以及Gamma分布的情況，為縱向數據基因關聯分析提供了強有力的分析工具。

本例中，以體重為應變量，以83個微衛星位點、測量時間及時間的平方為協變量，共85個協變量，X＝［1，time，time2，x1，x2，…，x83］，其中xj編碼為0和1，這是典型的縱向數據基因關聯分析。對所有變量均進行標準化，采用正態鏈接函數的GEE和pGEE進行分析。由于每只實驗小鼠在等距的時間間隔內重復測量了6次，兩次測量相隔時間越長，相應的相關關系越小，因此，假定作業相關矩陣為1階自回歸結構。截距項不進行懲罰。通過“CVfit”函數選擇出最優的懲罰參數為λ＝0．0212，進一步利用“PGEE”函數進行擬合，得到pGEE參數估計結果。GEE分析在“geepack”軟件包中實現，GEE和pGEE參數設置均采用軟件中的默認設置。

結 果

1．實驗小鼠體重一般情況

將每只小鼠的體重測量值按時間進行連接，得體重隨時間變化趨勢圖，如圖1?？梢娦∈笤?4周內的個體生長趨勢總體上一致，呈先快后慢的趨勢，但也存在個體間的差異，個體增長幅度可能不一致，其差別隨時間的變化而增加。進一步分析每個微衛星位點的單獨效應，對每個位點分別擬合GEE，得到單個位點的效應系數，如圖2，單個位點系數分布在［-0．176，0．099］之間，大部分位點的效應非常弱。

圖1 小鼠不同周次體重變化趨勢圖

圖2 單個QTL系數圖

2．GEE和pGEE分析

在考慮了時間效應的情況下，觀察GEE和pGEE對83個微衛星位點的參數估計散點圖（圖3），pGEE將大部分具有微弱效應的位點收縮為0，有效的進行了降維，而未懲罰的GEE卻不滿足稀疏性。GEE采用Wald卡方檢驗進行統計推斷，選出了14個有意義的QTL（檢驗水準），pGEE篩選出11個QTL，其參數估計結果見表2。比較兩方法共同篩選出的QTL參數估計值，可看出pGEE估計值明顯小于GEE，從QTL效應的微弱性而言，pGEE估計出的QTL效應更符合分子生物學特點。

圖3 83個微衛星位點的GEE和pGEE參數估計散點圖

在pGEE篩選出的11個QTL中，根據每個QTL系數的正負，認為D1M it211，D5M it158，D6M it275及D15M it29的突變為危險因素，將導致小鼠肥胖，其余7個QTL的變異對小鼠肥胖具有保護效應。

肥胖與2型糖尿病的發病顯著相關，將誘發糖尿病的發生。肥胖的發生、發展受多個微效基因及其復雜的交互作用控制，同時也受環境因子的調控，是一種復雜的多因素、多基因疾?。?4］。本例中通過pGEE分析篩選出的QTL，和Reifsnyder采用單因素方差分析和置換檢驗識別的體重相關的7個第一染色體位點相比［13］，共同識別的位點為D1M it211，D1M it411，D1M it76。陳峰采用線性混合效應模型（linear mixed model，LMM）的Lasso和SCAD懲罰對小鼠體重QTL進行了識別［15］，分別識別出14和3個QTL，和本文結果相比較，基于Lasso的LMM、基于SCAD的LMM以及本文的pGEE，三種不同懲罰模型共同識別出了D1M it411和D5M it158位點，提示這兩個位點的重要性，同時pGEE識別的其他9個QTL，也將為糖尿病發病相關QTL分子生物學研究提供重要的候選位點。

表2 GEE和pGEE參數估計值

總之，pGEE成功的對高維縱向基因數據進行了降維，其參數估計具有oracle性質，有效的實現了高維縱向基因數據的變量選擇，篩選出了有意義的遺傳變異。

討 論

pGEE分析借助作業相關矩陣考慮了不同時間點測量值之間的內部相關關系，通過基于SCAD的懲罰方法，實現了高維縱向數據的變量選擇。本文以糖尿病縱向數據的QTL分析為例，識別出糖尿病發病相關的QTL，說明了pGEE在縱向數據基因關聯性分析中應用的科學性和可行性。

基于Wald卡方檢驗的GEE統計推斷，其前提是在保證其他變量不變的情況下，估計變量的偏效應，當維度較高且樣本量較小時，基于偏效應的統計推斷無法揭示變量的真實效應?；趹土P的方法參數估計和變量選擇同時進行，pGEE中采用的SCAD懲罰具有一致性和oracle性質，保證了非零系數變量的正確推斷和選擇。在高維縱向基因關聯分析中，pGEE的變量選擇明顯優于GEE的參數估計和統計推斷。

pGEE變量選擇的一致性和稀疏性依賴于懲罰參數的選擇，基于交叉驗證的懲罰參數選擇方法容易出現過擬合，從而將無效變量選為有意義的變量，出現錯誤選擇［16］。Hyunkeun等通過數據模擬和真實數據分析，指出基于貝葉斯信息準則（bayesian information criterion，BIC）的懲罰參數選擇方法得到的參數估計值具有一致性［17］。因此，今后可進一步研究pGEE框架下的BIC懲罰參數選擇方法，選擇出最優的懲罰參數，以得到更為確切的變量選擇結果。需要注意的是，經pGEE篩選出的有意義的位點，還需要結合分子生物學實驗進一步驗證其生物學意義。

總之，pGEE作為縱向數據經典方法GEE的懲罰模型，一方面，延續了GEE在縱向數據分析中的重要優勢：只需設定一階和二階矩以及作業相關矩陣，避免了高維縱向數據分析中更為復雜的全似然函數計算；當作業相關矩陣指定不當時，仍能保持參數估計的一致性。另一方面，pGEE將GEE推廣到高維數據分析，既能考慮縱向數據之間的相關關系，又能有效的進行降維，識別出高維縱向基因數據的關聯位點，為高維縱向數據分析提供了強有力的分析方法。隨著分子生物學技術成本的降低，縱向基因數據必將日益增多，pGEE將在縱向數據基因關聯性分析中發揮及其重要的作用。

［1］M ccarthy M I，Zeggini E．Genome-w ide association studies in type 2 diabetes．Curr Diab Rep，2009，9（2）：164-171．

［2］Ehret GB．Genome-w ide association studies：contribution of genom ics to understanding blood pressure and essential hypertension．Curr Hypertens Rep，2010，12（1）：17-25．

［3］Sitlani CM，Rice K M，Lum ley T，et al．Generalized estimating equations for genome-w ide association studies using longitudinal phenotype data．Stat Med，2015，34（1）：118-130．

［4］Furlotte NA，Eskin E，Eyheramendy S．Genome-w ide association mapping w ith longitudinal data．Genet Epidem iol，2012，36（5）：463-471．

［5］ShiG，Rao D．Ignoring temporal trends in genetic effects substantially reduces power of quantitative trait linkage analysis．Genetic epidem iology，2008，32（1）：61-72．

［6］Tibshirani R．Regression shrinkage and selection via the lasso．Journal of the Royal Statistical Society．Series B（Methodological），1996：267-288．

［7］Fan J，Li R．Variable selection via nonconcave penalized likelihood and its oracle properties．Journal of the American statistical Association，2001，96（456）：1348-1360．

［8］Zou H．The adaptive lasso and its oracle properties．Journal of the A-merican statistical association，2006，101（476）：1418-1429．

［9］Liang KY，Zeger SL．Longitudinal data analysis using generalized linearmodels．Biometrika，1986，73（1）：13-22．

［10］陳啟光．縱向研究中重復測量資料的廣義估計方程分析．中國衛生統計，1995，12（1）：22-25．

［11］李洪艷，譚珊，高曉，等．基于廣義估計方程的嬰兒超重的影響因素分析．中國衛生統計，2016，33（2）：222-225．

［12］Wang L，Zhou J，Qu A．Penalized generalized estimating equations for high-dimensional longitudinal data analysis．Biometrics，2012，68（2）：353-360．

［13］Reifsnyder PC，Churchill G，Leiter EH．Maternal environment and genotype interact to establish diabesity in m ice．Genome Res，2000，10（10）：1568-1578．

［14］Reifsnyder PC，Leiter EH．Deconstructing and reconstructing obesityinduced diabetes（diabesity）in m ice．Diabetes，2002，51（3）：825-832．

［15］陳峰．線性混合效應模型的懲罰變量選擇．中國衛生信息管理雜志，2014，11（3）：278-284．

［16］Wang H，Li R，Tsai CL．Tuning parameter selectors for the smoothly clipped absolute deviation method．Biometrika，2007，94（3）：553-568．

［17］Cho H，Qu A．Model selection for correlated data with diverging number of parameters．Statistica Sinica，2013，23：901-927．

（責任編輯：張 悅）

The Application of Penalized Generalized Estimating Equations in Genetic Association w ith Longitudinal Data

Cao Hongyan，Zeng Ping，Li Zhi，et al（Department of Health Statistics，ShanxiMedical University（030001），Taiyuan）

Objective To explore the application of penalized generalized estimating equations（pGEE）in genetic association w ith longitudinal data，and provide new statistical solutions for longitudinal genetic data analysis．M ethods We applied the generalized estimating equations（GEE）and pGEEmethods to a type II diabetes dataset for quantitative trait locus identification．Results Several loci associated w ith the development of type IIdiabeteswere identified using the pGEEmethod，providing important candidatemakers for further biological validation．Conclusion The pGEEmethod provides a powerful tool for high dimensional longitudinal genetic association studies．

Penalized generalized estimating equations；Longitudinal data；SCAD；Genetic association studies；Quantitative trait locus

*：國家自然科學基金資助項目（30972553，31371336）

1．山西醫科大學衛生統計教研室（030001）

2．徐州醫科大學流行病與衛生統計學教研室

3．中北大學體育學院

4．美國密西根州立大學統計與概率系

△通信作者：張巖波，E-mail：sxmuzyb＠126．com