-计算机应用-
*基金项目：国家重点研发计划精准医学研究重点专项(2017 YFC0907200,2017 YFC0907201)#共同第一作者
△通信作者：赵亚玲,E-mail : zhaoyl666 @ xjtu. edu. cn ；颜虹，E-mail : yanhonge@ xjtu. edu. cn
西安交通大学医学部公共卫生学院流行病与卫生统计学系(710061)
吴晨璐#米白冰#陈方尧裴磊磊史青云赵亚玲△ 颜虹^
【提 要】 目的 SAS 软件中目前实现广义线性混合模型的过程步主要包括PROC GLIMMIX 和PROC NLMIXED ，两种方法在实际应用中各有侧重。本文介绍一个可以提高广义线性混合模型运行效率的SAS 宏程序％
HPGLIMMIX 的使用方法及其结果解读。方法通过实例数据，介绍% HPGLIMMIX 分析正态分布和二
项分布数据的过 程，并展示采用%HPGLIMMIX 分析大样本数据的性能优势。结果对于小样本正态分布和二项分布数据，采用％ HPGLIMMIX 和GLIMMIX 、NLMIXED 分析的用法基本一致。对于大样本数据,％ HPGLIMMIX 可进行模型拟合并可有 效节省时间及计算资源。结论 %HPGLIMMIX 可有效提升大样本数据的广义线性混合模型拟合的效率。NLMIXED 过 程可以快速准确地进行参数估计。
【关键词】 广义线性混合模型 % HPGLIMMIX GLIMMIX NLMIXED SAS 宏【中图分类号】R195. 1 【文献标识码】A DOI 10. 3969/j. issn. 1002 -3674.2021.01.036
队列和临床试验等医学研究中经常遇到重复测量 的纵向数据。此类数据不满足观测时点间的独立性假 设，故不宜采用传统的线性模型(linear models , LM ),
而应采用纳入随机效应项的混合模型(mixed model , MM )进行分析［1］。混合模型适用于结局为连续变量、
分类变量或等级变量的数据［2］。随着数据来源及收
集技术的不断发展，大样本队列、多中心临床研究日益
增多。 这些研究常收集十几万例研究对象的长期、多
次观测值，产生高维复杂纵向数据。广义线性混合模
型(generalized linear mixed model , GLMM )广泛应用
于这类数据的分析处理，但其参数估计可能会因为计
算内存不够而无法输出结果［3］，导致其应用受限。针
对上述情况，本文介绍可进行高效广义线性混合模型
拟合估计的SAS 宏程序“ % HPGLIMMIX ”［3］，展示其 在大样本数据参数估计上的优势, 并结合实际数据分
析案例探讨其具体使用方法, 供广大科研工作者进行
大样本数据的 GLMM 拟合及参数估计时参考使用。
广义线性混合模型简介
GLMM 是广义线性模型和线性混合模型的结合, 由于其可应用于指数分布家族的任意一种分布类型，且
引入了随机效应项来解释数据间的相关性、过度离散、 异质性等问题［4］，故可将连续变量(正态)或
离散变量
(二分类或多分类资料)作为反应变量，使其近20年来
在生物医学领域特别是纵向数据分析中广受欢迎［3］。
GLMM 在模型中引入随机效应项“,，一般表达式为［5］:
g (M ij ) = n , = x ；0
+ Z 訥
设计矩阵由固定效应X 与随机效应Z 两部分组
成，随机效应项"，可以解释由于不可测因子引起的类
间异质性和同一类内观测到的相关。
在随机效应存在的情况下，广义线性混合模型默
认通过伪似然法(pseudo-likelihood )来估计参数，也可 以选择最大似然法(maximum likelihood )来进行参数
估计。在SAS 统计软件包中，GLMM 可通过PROC GLMMIX 、PROC NLMIXED 及 PROC GENMOD 等过 程步实现, 以 PROC GLMMIX 最为常用。
SAS 宏程序％HPGLIMMIX 及其参数介绍PROC GLIMMIX 过 程 可 处 理 多 数 情 况 下 的 GLMM 拟合，但当涉及到大样本数据时，特别是随机 效应项水平较多时，该程序拟合常难以收敛。对于大
样本高维纵向数据，若结局变量为满足正态分布的连
续型变量时，可采用自SAS 9.2开始引入的PROC
HPMIXED 过程步构建线性混合模型； 但对结局为分类 和等级变量的大样本纵向数据的GLMM 建模,SAS 统
betal计软件包至今尚未提供相应的分析模块。2014年,Xie
L 和Madden L ［3］编写了“高性能广义线性混合模型宏
( % HPGLIMMIX) ” 来解决大样本数据 GLMM 建模的
问 题。 和 PROC GLIMMIX 过 程 步 类 似, %
HPGLIMMIX 宏也使用限制性残差伪似然法(restricted
pseudo-likelihood,REPL) 拟合模型, 但为适应数据量大 的特征,％ HPGLIMMIX 宏 参考了 HPMIXED 的 过 程［6］，应用了稀疏矩阵技术来解决混合效应，使其适用
于高性能计算（high performance ），改良了计算过程，提 高了计算效率，适用于观测值较多的大样本数据的
GLMM 拟合。
调用％HPGLIMMIX 宏语句的主要框架代码如下丁：
% Hpglimmix （ DATA = ＜ 数据集名称〉
STMTS = %Str （CLASS ＜分类变量〉；
MODEL ＜响应变量〉=〈解释变量1解释
变量 2…〉/SOLUTION ；
RANDOM ＜随机效应〉/SUBJECT = ＜分区
组变量〉；
ESTIMATE/LSMEANS 固定效应/ ＜需输
出的统计量〉•••；）,
ERROR = ＜分布类型〉丄INK = ＜链接函数〉,
TECH = ＜参数估计优化方法〉,OPTIONS = ＜空间分离的选项关键词〉
） RUN ；
在调用宏程序的过程中，各参数意义及其使用如下:
STMTS ：用％ Str 来指定分析步骤，定义方法与
PROC GILMMIX 类似，但残差分布类型和链接函数不 在MODEL 处定义，另由专门参数ERROR 和LINK 定义。
ERROR ：定义残差分布类型，具体选择见表1。当定 义ERROR = User 时，需要同时定义ERRVAR （方差函数）
和ERRDEV （偏差函数）；默认的分布类型为二项分布。
LINK ：定义链接函数，每种分布类型对应默认的
链接函数见表1。
TECH ：参数估计优化方法，默认为Newton-
Raphson Ridge 法（NRRIDG ），其他方法见表 2。
表2 GLMM 中常用的参数估计方法①
参数估计方法全称参数估计方法缩写（TECH ）
Levenberg-Marquardt
LEVMAR Trust Region
TRUREG
Newton-Raphson with Line Search
NEWRAP Newton-Raphson Ridge
NRRIDG
Quasi-Newton QUANEW
Double-Dogleg
DBLDOG
Conjugate Gradient CONGRA Nelder-Mead Simplex
NMSIMP
表1 GLMM 中常用分布及其链接函数L3,7-8-data pr input child gender /*输入妾璧儿童编号
其中性别变壘中1代Smale, 2代表Female ，分布类型
ERROR
Link
二兀分布（Binary ）Binary
Logit 二项分布（Binomial ）
Binomial / Bin / B Logit 多项有序分布（Multinomial ）Multinomial / Multi CLogit 正态分布（Normal ）Normal / Nor / N Identity/ID 泊松分布（Poisson ）
Poisson / Poi / P Log 负二项分布（Negbinomial ）Negbinomial / NB
Log
Gamma 分布(Gamma)Gamma / Gam / G Reciprocal 反高斯分布（Invgaussian ）Invgaussian / IG Power ( -2)几何分布(Geometric )Geometric
Log
用户自定义
User-specified
4个年龄分别的垂体中央到翼突上颔裂的距离*/
array yy yl-y4do tine=l to 4:
OPTIONS ：定义输出选项；PROCOPT 用来指定
HPMIXED 步骤的选项;MAXIT 用来定义最大迭代次
数，默认值为20； OFFSET 用来定义位移变量（offset频率补偿电路
variable ，默认无该变量），仅在采用泊松分布的链接函 数时使用。
SAS 宏程序％HPGLIMMIX 宏应用举例
以下应用实例分别展示:①该宏程序在分析结果 上与GLIMMIX 及NLMIXED 过程的一致性（例1和
例2）；②在大样本数据参数估计上的性能优势
（例 3）。
1. SAS 宏程序％ HPGLIMMIX 基本用法
基于儿童生长曲线研究数据[9]，比较几种方法的
分析结果。该数据系27名儿童（11女，16男）在8、
10、12、14岁的生长发育测量数据，假设反应变量y （儿 童垂体中央到翼突上颌裂的距离，单位为mm ）为正态 分布，自变量child 为儿童序号、gender 为性别（其中1
代表男性,0代表女性）、time 为测量次数、age 为年龄
（岁）。其数据步见图1 -a ，使用% HPGLIMMIX 宏的
程序步见图1-b，使用GLIMMIX 步骤的程序见图1-c, 使用NLMIXED 步骤的程序见图1-d 。
%^£-Iz**zx(data^pr> stntsfstr (
class child time:/*说BE 分类变量*/ model y=gender age / solution: /*构建槓型：将年龄，性别纳入核型*/ random intercept /*说明随机救应I 页：
age = time*2 + 6; y= yyltime] /*分别对应次教，
年龄和垂体中央到翼突上颔裂的距离*/
mean=24. s=2. 928X68:
/*构建複型：将年龄，性别纳入複型*/ random inter type=un sub=child;/*说明随机热 儿童编号*// 认分布类型为正态分布*/
生长曲线数据的GL1MMIX 程序步
»odel y normal (nean,s):/*说朋分布粪型：正态分布” random u " nornal (0.s2u)
图1生长曲线数据SAS
代码及日志
由各程序步运行日志图1-e、图1-f和图1-g可知，分析观测值较少的数据时,％HPGLIMMIX宏的运行时间(1秒)比GLIMMIX(实际时间0.14秒,CPU时间0.10秒)及NLMIXED(实际时间0.51秒,CPU时间0.40秒)长。
图2分别为采用%HPGLIMMIX宏、GLIMMIX和NLMIXED分析27名儿童相关数据的主要参数估计结果。可见，%HPGLIMMIX宏和GLIMMIX程序得出的模型及参数估计基本一致。
Covariance Parameter Estimates
Cov Parm Subject Estimate
UN(1,1)child32668
Residual 2.0495 a.儿童生长曲线研究，数据%HPGLIMM1X宏的主要运行结果:
Covariance Parameter Estimates
Solution for Fixed Effects
Effect Estimate Standard
Error DF t Value Pr>hl Alpha Lower Upper
BLK222Intercept17.70670.833910521.23V.00010.0516.053219.3602 gender-2.32100.7614105-3.05Q.00290.05-3.8308-0.8113 age0.6602□.061611051072V-00010.050.53800.7823 b.儿童生长曲线研沉数宏的主要运行结采:
Solution for Fixed Effects
Covariance Parameter Estimates
Cov Parm Subject Estimate
Standard
UN(1,1)child 3.2668 1.0720
Residual 2.04950.3240
c.儿童牛长曲线研究数据Glimmix过程的上要运行结果:
Covariance Parameter Estimates
Solutions for Fixed Effects
Effect Estimate
Standard
Error DF t Value Pr>|t|
Intercept17.70670.83392521.23<0001
gender-2.32100.76U80-3.050.0031
age0.66020.061S18010.72<0001
dJL童生长曲线研究数iRGlimmix过程的主要运行结果:
Solution for Fixed Effects
Parameter Estimates
Parameter Estimate
Standard
DF t Value Pr>Itl95%Confidence Limits Gradient to19.6350161.1S260.120.9040-311.62350.89-1.06137 tl-0.49650.109726-4.530.0001-0.7219-0.2710-41.7060 120.241726404.401 s2u-11I E-142641.4776 mean22.0947161.21260.140.8920•309.27353.461-06137 s0.3137.26..
e.儿童牛长曲线研究数据NLMIXED过汞
Parameter Estimate的主要运行結果：
-20.8642
图2生长曲线数据的主要结果
2.分类结局变量中%HPGLIMMIX的应用
基于对局部使用乳霜是否有助治愈感染的多中心临床试验数据［8］，介绍％HPGLIMMIX宏在二项分布结局变量分析中的应用。反应变量Y(治愈率)服从二项分布,X表示事件发生数(治愈),n表示试验总数。此数据集共包括8个临床试验中心，研究对象共273例，其中试验组130例，对照组143例。
基础桩%HPGLIMMIX宏的使用见图3-a，基于GLIMMIX的程序见图3-b，基于NLMIXED的程序见图3-c。此例中，％HPGLIMMIX宏因使用的默认起始值算法不同，需要更多迭代次数来达到收敛。三种方法运行的结果基本相同，在此不再赘述。
上述两例说明%HPGLIMMIX宏的基本用法，其分析结果与GLIMMIX步骤一致，但在分析小样本数据时，其节省时间、提高效率的优势难以体现［3］。
%hpglijmix(dat a=inf e ct i o n,
stnrts=Xstr(
class clinic treatment:
/*说明分类变里*/
Model x/n=treatment/solution:
/*构翟模型：将方案分组纳入模型材
random intercept/subject=clinic )/*说明随机效应项：临床试验中心编号刃err or=b i nonial,
定义分布类型为二项分布*/
link=locit proc gliMMix data=infection;
class clinic treatment(ref=‘O')
/*说明分类变里*/
model x/n=treatment/s
dist=binomial link=logit ddfm=residual
/
*构建複型：将方案分组纳入複型；
走义分布类型为二顷分布*/
random intercept/subject=clinic
/*说明随机效应项：临床试验中心编号*/
-proc nlaixed data= infection_2
qpoints= 20tech=newrap;
pains beta0=-l-1464
betal =0.7244s2u=2.0327
/*根据gliiuuix步骤结果设定初始参数*/
eta=betaO+betal*treatment+u;
/*摸型构建：将年龄，性别纳入模型仪
聂铭杰
mu=exp(eta)/(1+exp(eta));
/*说明数：Logitech*/
model pre~binonial(n,hlu):
/*说明分齊类型为二项分布粒
random u~nor>al(0,s2u)subject=clinic
说明随机效应项：临床试验中心编号*/
a.二项分布数据的%HPGLIMMIX步b-二项分布数据的GL1MMIXH序步d二项分布数据的NLMIXED程序步
图3二项分布数据的SAS代码
3.大样本数据中%HPGLIMMIX的应用
为体现%HPGLIMMIX宏在大样本数据分析时的优势，采用“中国健康与营养调查(Chinese health and nutrition surveys,CHNS)”1991年~2011年的营养调查数据［10-13］进行实例验证。该数据属于大样本
高维纵向数据，故引入随机效应项并使用混合模型进行分析［14］。在分析成年人的能量摄入变化时，先使用GLIMMIX构建模型，反应变量为日能量摄入量(d3kcal)，固定效应包括调查年份(wave)、性别(gender)、城乡(urban_cat4)、教育程度(edu_cat4)、地区(area)、收入(income)、年龄(age_group6)、民族(nationality_cat2)及婚姻状况(marry_cat3)，个体ID (indiv)作为随机效应项，链接函数为identity；使用LSMEANS对调整后各调查年份的能量摄入进行估计。SAS程序、日志见图4-a和图4-b。
由于此数据集样本量大，导致GLIMMIX过程因数据溢出无法拟合。NLMIXED与GLIMMIX情况类似，均未能拟合模型。下面采用%HPGLIMMIX宏进行模型构建及参数估计，SAS程序见图4-c。％HPGLIMMIX宏的主要运行结果及能量摄入量的调整值见图4-d和图4-e。
通过上述分析可见,％HPGLIMMIX宏可以进行大样本的模型拟合，解决实际工作中因数据量过大导致内存不够、模型无法拟合的问题。
_proc gliuix dat a=aock.
class indiv wave gende cat4
area incoae edu.cati a;e_{coup6
nationality_cat2narry_cat3: /*说囲分負变塵*/
node!.d3kcal=・ave gender urban_cat4
area incme edu_cat4age_g
nat ionality^cat2narry_cat /*构建播型；入模塑*/
random intercept/subject=indiv
/*说册麻机敘応I页；*体ID*/
%ApgliMix(dat work.chns_nutr,
stnt s=%s t r(
class indiv wave gender urban_cat4
edu_cat4area inc upS
nat i o nality_cit2marxy.oa't3;
/
*说明分类殘■*/
model d3kcal=vave gender uzban_cat4
edu_cat4area income age_group6
nationality.eat2marry_cat3/so lut ion:
/*构建fli型；
random intercept/sub ject=indiv
/*说阴pa机敘应项；个体“*/
lsneans wave/cl;
/*对调盤后SiBSS份的入及!HfllS间SB行fSi+“
error=noraal,lirit=identity
/^HPGLIMMIX宏獸认为二顶分布，St说明正态分布的魅按幽数"
刮膜棒
d.采用％HPGLIMMIX宏分术丿FCHNS数据的主要结果；
Covariance Parameter Estimates
c.CHNS研究数据的％HPGL1MM1X程序步
e.米用%HPGLIMMIX宏分析CHNS数据的主耍结采：
Wave Least Squares Means
Covariance Parameter Estimates
Cov Parm Subject Estimate
Intercept indiv51603
Residual322546
图4CHNS数据的SAS代码、日志及结果
讨论
在大样本高维复杂纵向数据的分析中,％HPGLIMMIX宏解决了SAS自带的HPMIXED语句无法拟合离散型数据GLMM模型的问题。其参数估计默认采用的双重迭代线性化法(doubly iterativelinearization method,DILM)，可在保证足够参数估计精度的情况下大量节约内存和运行时间，降低时间成本。既往研究证明在高维复杂数据的分析实践中,%HPGLIMMIX宏可以节约高达90%的运行时间，特别是在考虑随机效应时，节省时间的效果更明显［3］。但需注意的是，尽管在模型参数估计时使用了相同的伪似然算法DILM,%HPGLIMMIX宏和GLIMMIX过程在使用中仍存在一些差异［3］。此外，尽管DILM和MLE(maximum likelihood estimation)在方法学上差别不大，但MLE无法处理大样本数据，且二者在参数估计上的差别仅在样本量特别小等极端情况下才会出现［15］。
综上所述,%HPGLIMMIX宏是一个新的、高效、可靠的广义线性混合模型建模方法，适合进行大样本纵向数据的分析，研究者可根据数据的特点合理选用。
参考文献
［1］汤宁，宋秋月，易东，等.医学纵向数据建模方法及其统计分析策略.中国卫生统计,2019,36(3)：441-444.
［2］刘星星.广义线性混合模型在二分类纵向数据中的探索研究.复旦大学,2014.[3]Xie L,Madden LV.%HPGLIMMIX：A High-Performance SAS
Macro for GLMM Estimation.Journal of Statistical Software,2014,8
(58)：1-25.
[4]王玲•广义线性混合模型在多中心中药临床试验中的应用•北京
中医药大学,2013.
[5]李丽霞，郜艳晖，张丕德，等•广义线性混合效应模型及其应用•现
代预防医学,2007,34(11)：2103-2104.
[6]Kiernan K,Tao J,Gibbs P.Tips and Strategies for Mixed Modeling
with SAS/STAT?Procedures.https：//support,sas/resources/ papers/proceedings12/332-2012.pdf.
[7]杨宇弘•广义线性混合模型在死亡数据分析中的应用•上海交通
大学,2015.
[8]罗天娥•非正态及非线性重复测量资料分析模型及其医学应用.
山西医科大学,2007.
[9]Pothoff RF,Roy SN.A Generalized Multivariate Analysis of
Variance Model Useful Especially for Growth Curve Problems.
Biometrika,1964.
[10]Popkin BM,Du S,Zhai F,et al.Cohort Profile：The China Health and
Nutrition Survey-monitoring and Understanding Socio-economic and Health Change in China,1989-2011.International Journal of Epidemiology,2010,39(6)：1435-1440.
[11]Zhang B,Zhai FY,Du SF,et al.The China Health and Nutrition
Survey,1989-2011.Obesity Reviews,2014,15：2-7.
[12]朴建华，张坚，赵文华，等.中国居民营养与健康状况调查的质量
控制•中华流行病学杂志,2005,(7)：474477.
[13]杨晓光，孔灵芝，翟凤英，等.中国居民营养与健康状况调查的总
体方案.中华流行病学杂志,2005,(7)：471-474.
[14]Su C,Zhao J,Wu Y,et al.Temporal Trends in Dietary Macronutrient
Intakes among Adults in Rural China from1991to2011：Findings from the CHNS.Nutrients,2017,9(3)：227.
[15]Stroup.Generalized Linear Mixed Models：Modern Concepts,
Methods and Applications//Statistical Science.Chapman&Hall/ CRC,2012.
(责任编辑：刘壮
)
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议