西安交通大學(xué)醫(yī)學(xué)部公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(710061)

吳晨璐# 米白冰# 陳方堯 裴磊磊 史青云 趙亞玲△ 顏 虹△

【提 要】 目的 SAS軟件中目前實(shí)現(xiàn)廣義線性混合模型的過(guò)程步主要包括PROC GLIMMIX和PROC NLMIXED，兩種方法在實(shí)際應(yīng)用中各有側(cè)重。本文介紹一個(gè)可以提高廣義線性混合模型運(yùn)行效率的SAS宏程序%HPGLIMMIX的使用方法及其結(jié)果解讀。方法 通過(guò)實(shí)例數(shù)據(jù)，介紹%HPGLIMMIX分析正態(tài)分布和二項(xiàng)分布數(shù)據(jù)的過(guò)程，并展示采用%HPGLIMMIX分析大樣本數(shù)據(jù)的性能優(yōu)勢(shì)。結(jié)果 對(duì)于小樣本正態(tài)分布和二項(xiàng)分布數(shù)據(jù)，采用%HPGLIMMIX和GLIMMIX、NLMIXED分析的用法基本一致。對(duì)于大樣本數(shù)據(jù)，%HPGLIMMIX可進(jìn)行模型擬合并可有效節(jié)省時(shí)間及計(jì)算資源。結(jié)論 %HPGLIMMIX可有效提升大樣本數(shù)據(jù)的廣義線性混合模型擬合的效率。NLMIXED過(guò)程可以快速準(zhǔn)確地進(jìn)行參數(shù)估計(jì)。

隊(duì)列和臨床試驗(yàn)等醫(yī)學(xué)研究中經(jīng)常遇到重復(fù)測(cè)量的縱向數(shù)據(jù)。此類數(shù)據(jù)不滿足觀測(cè)時(shí)點(diǎn)間的獨(dú)立性假設(shè)，故不宜采用傳統(tǒng)的線性模型(linear models，LM)，而應(yīng)采用納入隨機(jī)效應(yīng)項(xiàng)的混合模型(mixed model，MM)進(jìn)行分析[1]?；旌夏Ｐ瓦m用于結(jié)局為連續(xù)變量、分類變量或等級(jí)變量的數(shù)據(jù)[2]。隨著數(shù)據(jù)來(lái)源及收集技術(shù)的不斷發(fā)展，大樣本隊(duì)列、多中心臨床研究日益增多。這些研究常收集十幾萬(wàn)例研究對(duì)象的長(zhǎng)期、多次觀測(cè)值，產(chǎn)生高維復(fù)雜縱向數(shù)據(jù)。廣義線性混合模型(generalized linear mixed model，GLMM)廣泛應(yīng)用于這類數(shù)據(jù)的分析處理，但其參數(shù)估計(jì)可能會(huì)因?yàn)橛?jì)算內(nèi)存不夠而無(wú)法輸出結(jié)果[3]，導(dǎo)致其應(yīng)用受限。針對(duì)上述情況，本文介紹可進(jìn)行高效廣義線性混合模型擬合估計(jì)的SAS宏程序“%HPGLIMMIX”[3]，展示其在大樣本數(shù)據(jù)參數(shù)估計(jì)上的優(yōu)勢(shì)，并結(jié)合實(shí)際數(shù)據(jù)分析案例探討其具體使用方法，供廣大科研工作者進(jìn)行大樣本數(shù)據(jù)的GLMM擬合及參數(shù)估計(jì)時(shí)參考使用。

廣義線性混合模型簡(jiǎn)介

GLMM是廣義線性模型和線性混合模型的結(jié)合，由于其可應(yīng)用于指數(shù)分布家族的任意一種分布類型，且引入了隨機(jī)效應(yīng)項(xiàng)來(lái)解釋數(shù)據(jù)間的相關(guān)性、過(guò)度離散、異質(zhì)性等問(wèn)題[4]，故可將連續(xù)變量(正態(tài))或離散變量(二分類或多分類資料)作為反應(yīng)變量，使其近20年來(lái)在生物醫(yī)學(xué)領(lǐng)域特別是縱向數(shù)據(jù)分析中廣受歡迎[3]。

GLMM在模型中引入隨機(jī)效應(yīng)項(xiàng)ui，一般表達(dá)式為[5]：

設(shè)計(jì)矩陣由固定效應(yīng)X與隨機(jī)效應(yīng)Z兩部分組成，隨機(jī)效應(yīng)項(xiàng)ui可以解釋由于不可測(cè)因子引起的類間異質(zhì)性和同一類內(nèi)觀測(cè)到的相關(guān)。

在隨機(jī)效應(yīng)存在的情況下，廣義線性混合模型默認(rèn)通過(guò)偽似然法(pseudo-likelihood)來(lái)估計(jì)參數(shù)，也可以選擇最大似然法(maximum likelihood)來(lái)進(jìn)行參數(shù)估計(jì)。在SAS統(tǒng)計(jì)軟件包中，GLMM可通過(guò)PROC GLMMIX、PROC NLMIXED及PROC GENMOD等過(guò)程步實(shí)現(xiàn)，以PROC GLMMIX最為常用。

SAS宏程序%HPGLIMMIX及其參數(shù)介紹

PROC GLIMMIX過(guò)程可處理多數(shù)情況下的GLMM擬合，但當(dāng)涉及到大樣本數(shù)據(jù)時(shí)，特別是隨機(jī)效應(yīng)項(xiàng)水平較多時(shí)，該程序擬合常難以收斂。對(duì)于大樣本高維縱向數(shù)據(jù)，若結(jié)局變量為滿足正態(tài)分布的連續(xù)型變量時(shí)，可采用自SAS 9.2開(kāi)始引入的PROC HPMIXED過(guò)程步構(gòu)建線性混合模型；但對(duì)結(jié)局為分類和等級(jí)變量的大樣本縱向數(shù)據(jù)的GLMM建模，SAS統(tǒng)計(jì)軟件包至今尚未提供相應(yīng)的分析模塊。2014年，Xie L和Madden L[3]編寫(xiě)了“高性能廣義線性混合模型宏(%HPGLIMMIX)”來(lái)解決大樣本數(shù)據(jù)GLMM建模的問(wèn)題。和PROC GLIMMIX過(guò)程步類似，%HPGLIMMIX宏也使用限制性殘差偽似然法(restricted pseudo-likelihood，REPL)擬合模型，但為適應(yīng)數(shù)據(jù)量大的特征，%HPGLIMMIX宏參考了HPMIXED的過(guò)程[6]，應(yīng)用了稀疏矩陣技術(shù)來(lái)解決混合效應(yīng)，使其適用于高性能計(jì)算(high performance)，改良了計(jì)算過(guò)程，提高了計(jì)算效率，適用于觀測(cè)值較多的大樣本數(shù)據(jù)的GLMM擬合。

調(diào)用%HPGLIMMIX宏語(yǔ)句的主要框架代碼如下[3]：

%Hpglimmix(DATA=〈數(shù)據(jù)集名稱〉

STMTS=%Str(

CLASS〈分類變量〉；

MODEL〈響應(yīng)變量〉=〈解釋變量1 解釋變量2…〉/SOLUTION;

RANDOM〈隨機(jī)效應(yīng)〉/SUBJECT=〈分區(qū)組變量〉；

ESTIMATE/LSMEANS固定效應(yīng)/〈需輸出的統(tǒng)計(jì)量〉…；)，

ERROR=〈分布類型〉，LINK=〈鏈接函數(shù)〉，

TECH=〈參數(shù)估計(jì)優(yōu)化方法〉，

OPTIONS=〈空間分離的選項(xiàng)關(guān)鍵詞〉

)RUN；

在調(diào)用宏程序的過(guò)程中，各參數(shù)意義及其使用如下：

STMTS：用%Str來(lái)指定分析步驟，定義方法與PROC GILMMIX類似，但殘差分布類型和鏈接函數(shù)不在MODEL處定義，另由專門參數(shù)ERROR和LINK定義。

ERROR：定義殘差分布類型，具體選擇見(jiàn)表1。當(dāng)定義ERROR=User時(shí)，需要同時(shí)定義ERRVAR(方差函數(shù))和ERRDEV(偏差函數(shù))；默認(rèn)的分布類型為二項(xiàng)分布。

LINK：定義鏈接函數(shù)，每種分布類型對(duì)應(yīng)默認(rèn)的鏈接函數(shù)見(jiàn)表1。

表1 GLMM中常用分布及其鏈接函數(shù)[3，7-8]

TECH：參數(shù)估計(jì)優(yōu)化方法，默認(rèn)為Newton-Raphson Ridge法(NRRIDG)，其他方法見(jiàn)表2。

表2 GLMM中常用的參數(shù)估計(jì)方法[3]

OPTIONS：定義輸出選項(xiàng)；PROCOPT用來(lái)指定HPMIXED步驟的選項(xiàng)；MAXIT用來(lái)定義最大迭代次數(shù)，默認(rèn)值為20；OFFSET用來(lái)定義位移變量(offset variable，默認(rèn)無(wú)該變量)，僅在采用泊松分布的鏈接函數(shù)時(shí)使用。

SAS宏程序%HPGLIMMIX宏應(yīng)用舉例

以下應(yīng)用實(shí)例分別展示：①該宏程序在分析結(jié)果上與GLIMMIX及NLMIXED過(guò)程的一致性(例1和例2)；②在大樣本數(shù)據(jù)參數(shù)估計(jì)上的性能優(yōu)勢(shì)(例3)。

1.SAS宏程序%HPGLIMMIX基本用法

基于兒童生長(zhǎng)曲線研究數(shù)據(jù)[9]，比較幾種方法的分析結(jié)果。該數(shù)據(jù)系27名兒童(11女，16男)在8、10、12、14歲的生長(zhǎng)發(fā)育測(cè)量數(shù)據(jù)，假設(shè)反應(yīng)變量y(兒童垂體中央到翼突上頜裂的距離，單位為mm)為正態(tài)分布，自變量child為兒童序號(hào)、gender為性別(其中1代表男性，0代表女性)、time為測(cè)量次數(shù)、age為年齡(歲)。其數(shù)據(jù)步見(jiàn)圖1-a，使用%HPGLIMMIX宏的程序步見(jiàn)圖1-b，使用GLIMMIX步驟的程序見(jiàn)圖1-c，使用NLMIXED步驟的程序見(jiàn)圖1-d。

圖1 生長(zhǎng)曲線數(shù)據(jù)SAS代碼及日志

由各程序步運(yùn)行日志圖1-e、圖1-f和圖1-g可知，分析觀測(cè)值較少的數(shù)據(jù)時(shí)，%HPGLIMMIX宏的運(yùn)行時(shí)間(1秒)比GLIMMIX(實(shí)際時(shí)間0.14秒，CPU時(shí)間0.10秒)及NLMIXED(實(shí)際時(shí)間0.51秒，CPU時(shí)間0.40秒)長(zhǎng)。

圖2分別為采用%HPGLIMMIX宏、GLIMMIX和NLMIXED分析27名兒童相關(guān)數(shù)據(jù)的主要參數(shù)估計(jì)結(jié)果。可見(jiàn)，%HPGLIMMIX宏和GLIMMIX程序得出的模型及參數(shù)估計(jì)基本一致。

圖2 生長(zhǎng)曲線數(shù)據(jù)的主要結(jié)果

2.分類結(jié)局變量中%HPGLIMMIX的應(yīng)用

基于對(duì)局部使用乳霜是否有助治愈感染的多中心臨床試驗(yàn)數(shù)據(jù)[8]，介紹%HPGLIMMIX宏在二項(xiàng)分布結(jié)局變量分析中的應(yīng)用。反應(yīng)變量Y(治愈率)服從二項(xiàng)分布，X表示事件發(fā)生數(shù)(治愈)，n表示試驗(yàn)總數(shù)。此數(shù)據(jù)集共包括8個(gè)臨床試驗(yàn)中心，研究對(duì)象共273例，其中試驗(yàn)組130例，對(duì)照組143例。

%HPGLIMMIX宏的使用見(jiàn)圖3-a，基于GLIMMIX的程序見(jiàn)圖3-b，基于NLMIXED的程序見(jiàn)圖3-c。此例中，%HPGLIMMIX宏因使用的默認(rèn)起始值算法不同，需要更多迭代次數(shù)來(lái)達(dá)到收斂。三種方法運(yùn)行的結(jié)果基本相同，在此不再贅述。

圖3 二項(xiàng)分布數(shù)據(jù)的SAS代碼

上述兩例說(shuō)明%HPGLIMMIX宏的基本用法，其分析結(jié)果與GLIMMIX步驟一致，但在分析小樣本數(shù)據(jù)時(shí)，其節(jié)省時(shí)間、提高效率的優(yōu)勢(shì)難以體現(xiàn)[3]。

3.大樣本數(shù)據(jù)中%HPGLIMMIX的應(yīng)用

為體現(xiàn)%HPGLIMMIX宏在大樣本數(shù)據(jù)分析時(shí)的優(yōu)勢(shì)，采用“中國(guó)健康與營(yíng)養(yǎng)調(diào)查(Chinese health and nutrition surveys，CHNS)”1991年～2011年的營(yíng)養(yǎng)調(diào)查數(shù)據(jù)[10-13]進(jìn)行實(shí)例驗(yàn)證。該數(shù)據(jù)屬于大樣本高維縱向數(shù)據(jù)，故引入隨機(jī)效應(yīng)項(xiàng)并使用混合模型進(jìn)行分析[14]。在分析成年人的能量攝入變化時(shí)，先使用GLIMMIX構(gòu)建模型，反應(yīng)變量為日能量攝入量(d3kcal)，固定效應(yīng)包括調(diào)查年份(wave)、性別(gender)、城鄉(xiāng)(urban_cat4)、教育程度(edu_cat4)、地區(qū)(area)、收入(income)、年齡(age_group6)、民族(nationality_cat2)及婚姻狀況(marry_cat3)，個(gè)體ID(indiv)作為隨機(jī)效應(yīng)項(xiàng)，鏈接函數(shù)為identity；使用LSMEANS對(duì)調(diào)整后各調(diào)查年份的能量攝入進(jìn)行估計(jì)。SAS程序、日志見(jiàn)圖4-a和圖4-b。

由于此數(shù)據(jù)集樣本量大，導(dǎo)致GLIMMIX過(guò)程因數(shù)據(jù)溢出無(wú)法擬合。NLMIXED與GLIMMIX情況類似，均未能擬合模型。下面采用%HPGLIMMIX宏進(jìn)行模型構(gòu)建及參數(shù)估計(jì)，SAS程序見(jiàn)圖4-c。%HPGLIMMIX宏的主要運(yùn)行結(jié)果及能量攝入量的調(diào)整值見(jiàn)圖4-d和圖4-e。

圖4 CHNS數(shù)據(jù)的SAS代碼、日志及結(jié)果

通過(guò)上述分析可見(jiàn)，%HPGLIMMIX宏可以進(jìn)行大樣本的模型擬合，解決實(shí)際工作中因數(shù)據(jù)量過(guò)大導(dǎo)致內(nèi)存不夠、模型無(wú)法擬合的問(wèn)題。

討 論

在大樣本高維復(fù)雜縱向數(shù)據(jù)的分析中，%HPGLIMMIX宏解決了SAS自帶的HPMIXED語(yǔ)句無(wú)法擬合離散型數(shù)據(jù)GLMM模型的問(wèn)題。其參數(shù)估計(jì)默認(rèn)采用的雙重迭代線性化法(doubly iterativelinearization method，DILM)，可在保證足夠參數(shù)估計(jì)精度的情況下大量節(jié)約內(nèi)存和運(yùn)行時(shí)間，降低時(shí)間成本。既往研究證明在高維復(fù)雜數(shù)據(jù)的分析實(shí)踐中，%HPGLIMMIX宏可以節(jié)約高達(dá)90%的運(yùn)行時(shí)間，特別是在考慮隨機(jī)效應(yīng)時(shí)，節(jié)省時(shí)間的效果更明顯[3]。但需注意的是，盡管在模型參數(shù)估計(jì)時(shí)使用了相同的偽似然算法DILM，%HPGLIMMIX宏和GLIMMIX過(guò)程在使用中仍存在一些差異[3]。此外，盡管DILM和MLE(maximum likelihood estimation)在方法學(xué)上差別不大，但MLE無(wú)法處理大樣本數(shù)據(jù)，且二者在參數(shù)估計(jì)上的差別僅在樣本量特別小等極端情況下才會(huì)出現(xiàn)[15]。

綜上所述，%HPGLIMMIX宏是一個(gè)新的、高效、可靠的廣義線性混合模型建模方法，適合進(jìn)行大樣本縱向數(shù)據(jù)的分析，研究者可根據(jù)數(shù)據(jù)的特點(diǎn)合理選用。