SAS宏程序%HPGLIMMIX在大样本数据广义线性混合模型参数估计中的应用*

西安交通大学医学部公共卫生学院流行病与卫生统计学系(710061)

吴晨璐# 米白冰# 陈方尧 裴磊磊 史青云 赵亚玲△ 颜 虹△

【提 要】 目的 SAS软件中目前实现广义线性混合模型的过程步主要包括PROC GLIMMIX和PROC NLMIXED，两种方法在实际应用中各有侧重。本文介绍一个可以提高广义线性混合模型运行效率的SAS宏程序%HPGLIMMIX的使用方法及其结果解读。方法 通过实例数据，介绍%HPGLIMMIX分析正态分布和二项分布数据的过程，并展示采用%HPGLIMMIX分析大样本数据的性能优势。结果 对于小样本正态分布和二项分布数据，采用%HPGLIMMIX和GLIMMIX、NLMIXED分析的用法基本一致。对于大样本数据，%HPGLIMMIX可进行模型拟合并可有效节省时间及计算资源。结论 %HPGLIMMIX可有效提升大样本数据的广义线性混合模型拟合的效率。NLMIXED过程可以快速准确地进行参数估计。

队列和临床试验等医学研究中经常遇到重复测量的纵向数据。此类数据不满足观测时点间的独立性假设，故不宜采用传统的线性模型(linear models，LM)，而应采用纳入随机效应项的混合模型(mixed model，MM)进行分析[1]。混合模型适用于结局为连续变量、分类变量或等级变量的数据[2]。随着数据来源及收集技术的不断发展，大样本队列、多中心临床研究日益增多。这些研究常收集十几万例研究对象的长期、多次观测值，产生高维复杂纵向数据。广义线性混合模型(generalized linear mixed model，GLMM)广泛应用于这类数据的分析处理，但其参数估计可能会因为计算内存不够而无法输出结果[3]，导致其应用受限。针对上述情况，本文介绍可进行高效广义线性混合模型拟合估计的SAS宏程序“%HPGLIMMIX”[3]，展示其在大样本数据参数估计上的优势，并结合实际数据分析案例探讨其具体使用方法，供广大科研工作者进行大样本数据的GLMM拟合及参数估计时参考使用。

广义线性混合模型简介

GLMM是广义线性模型和线性混合模型的结合，由于其可应用于指数分布家族的任意一种分布类型，且引入了随机效应项来解释数据间的相关性、过度离散、异质性等问题[4]，故可将连续变量(正态)或离散变量(二分类或多分类资料)作为反应变量，使其近20年来在生物医学领域特别是纵向数据分析中广受欢迎[3]。

GLMM在模型中引入随机效应项ui，一般表达式为[5]：

设计矩阵由固定效应X与随机效应Z两部分组成，随机效应项ui可以解释由于不可测因子引起的类间异质性和同一类内观测到的相关。

在随机效应存在的情况下，广义线性混合模型默认通过伪似然法(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统计软件包至今尚未提供相应的分析模块。2014年，Xie L和Madden L[3]编写了“高性能广义线性混合模型宏(%HPGLIMMIX)”来解决大样本数据GLMM建模的问题。和PROC GLIMMIX过程步类似，%HPGLIMMIX宏也使用限制性残差伪似然法(restricted pseudo-likelihood，REPL)拟合模型，但为适应数据量大的特征，%HPGLIMMIX宏参考了HPMIXED的过程[6]，应用了稀疏矩阵技术来解决混合效应，使其适用于高性能计算(high performance)，改良了计算过程，提高了计算效率，适用于观测值较多的大样本数据的GLMM拟合。

调用%HPGLIMMIX宏语句的主要框架代码如下[3]：

%Hpglimmix(DATA=〈数据集名称〉

STMTS=%Str(

CLASS〈分类变量〉；

MODEL〈响应变量〉=〈解释变量1 解释变量2…〉/SOLUTION;

RANDOM〈随机效应〉/SUBJECT=〈分区组变量〉；

ESTIMATE/LSMEANS固定效应/〈需输出的统计量〉…；)，

ERROR=〈分布类型〉，LINK=〈链接函数〉，

TECH=〈参数估计优化方法〉，

OPTIONS=〈空间分离的选项关键词〉

)RUN；

在调用宏程序的过程中，各参数意义及其使用如下：

STMTS：用%Str来指定分析步骤，定义方法与PROC GILMMIX类似，但残差分布类型和链接函数不在MODEL处定义，另由专门参数ERROR和LINK定义。

ERROR：定义残差分布类型，具体选择见表1。当定义ERROR=User时，需要同时定义ERRVAR(方差函数)和ERRDEV(偏差函数)；默认的分布类型为二项分布。

LINK：定义链接函数，每种分布类型对应默认的链接函数见表1。

表1 GLMM中常用分布及其链接函数[3，7-8]

TECH：参数估计优化方法，默认为Newton-Raphson Ridge法(NRRIDG)，其他方法见表2。

表2 GLMM中常用的参数估计方法[3]

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。

图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程序得出的模型及参数估计基本一致。

图2 生长曲线数据的主要结果

2.分类结局变量中%HPGLIMMIX的应用

基于对局部使用乳霜是否有助治愈感染的多中心临床试验数据[8]，介绍%HPGLIMMIX宏在二项分布结局变量分析中的应用。反应变量Y(治愈率)服从二项分布，X表示事件发生数(治愈)，n表示试验总数。此数据集共包括8个临床试验中心，研究对象共273例，其中试验组130例，对照组143例。

%HPGLIMMIX宏的使用见图3-a，基于GLIMMIX的程序见图3-b，基于NLMIXED的程序见图3-c。此例中，%HPGLIMMIX宏因使用的默认起始值算法不同，需要更多迭代次数来达到收敛。三种方法运行的结果基本相同，在此不再赘述。

图3 二项分布数据的SAS代码

上述两例说明%HPGLIMMIX宏的基本用法，其分析结果与GLIMMIX步骤一致，但在分析小样本数据时，其节省时间、提高效率的优势难以体现[3]。

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。

图4 CHNS数据的SAS代码、日志及结果

通过上述分析可见，%HPGLIMMIX宏可以进行大样本的模型拟合，解决实际工作中因数据量过大导致内存不够、模型无法拟合的问题。

讨 论

在大样本高维复杂纵向数据的分析中，%HPGLIMMIX宏解决了SAS自带的HPMIXED语句无法拟合离散型数据GLMM模型的问题。其参数估计默认采用的双重迭代线性化法(doubly iterativelinearization method，DILM)，可在保证足够参数估计精度的情况下大量节约内存和运行时间，降低时间成本。既往研究证明在高维复杂数据的分析实践中，%HPGLIMMIX宏可以节约高达90%的运行时间，特别是在考虑随机效应时，节省时间的效果更明显[3]。但需注意的是，尽管在模型参数估计时使用了相同的伪似然算法DILM，%HPGLIMMIX宏和GLIMMIX过程在使用中仍存在一些差异[3]。此外，尽管DILM和MLE(maximum likelihood estimation)在方法学上差别不大，但MLE无法处理大样本数据，且二者在参数估计上的差别仅在样本量特别小等极端情况下才会出现[15]。

综上所述，%HPGLIMMIX宏是一个新的、高效、可靠的广义线性混合模型建模方法，适合进行大样本纵向数据的分析，研究者可根据数据的特点合理选用。