成长曲线模型原理及其应用

傅 强

美国塔夫茨大学文理学院社区健康系，美国马萨诸塞州梅德福市，02155

在社会科学及医学领域中，数据往往具有分级嵌套结构。例如在医学研究中，研究者先随机选择一定数量的医院，再从入选的每家医院中随机选取一定数量的医生，随后从入选的每个医生接诊的病人中随机选取一定数量的病人，最后对每位病人每个月随访一次，共随访三次。如果将每一次随访的数据确定为第一水平分析单位，则病人可确定为第二水平分析单位，医生为第三水平分析单位，医院为第四水平分析单位。每一水平的数据代表相应水平的独特特征，下一水平的数据嵌套(nested)或聚集(clustered)在上一水平数据结构中。又如在社会科学研究中，研究者先随机选择一定数量的城市，再从入选的每个城市中随机选取一定数量的社区，随后从入选的每个社区中随机选取一定数量的居民，最后对每位居民每周随访一次，多次随访。由此城市、社区、居民、每次随访构成了不同的水平分析单位，每一水平的数据都代表相应水平的独特特征。如果研究者仅仅随机选取病人或居民并多次随访，收集相关数据，这种数据结构则形成最简单的两级嵌套结构，也可以被视为是最简单的聚集数据(clustered data)。上述这些数据被称为纵向数据(longitudinal data)，也称为队列数据(cohort data)或面板数据(panel data)。

纵向数据包含对随机抽样个体的特征多次重复的观测值。这种研究设计的优点是可以发现一些特征随时间变化的趋势，例如病人临床转归的结局或居民生命质量和幸福感随时间变化的趋势以及影响其变化的因素，从而回答许多社会科学和医学领域中的问题。例如居民的生命质量在开始追踪调查时是否有差异？随着时间推移生命质量变化的轨迹(trajectory)怎样？不同居民生命质量变化的轨迹有区别吗？什么因素可能解释居民生命质量平均变化趋势的差异？但是，这种研究设计对运用一般线性模型(general linear models, GLM)的分析方法提出了新的挑战。一般线性模型假设任何两个观测值之间都是独立的，它可以用来分析个体水平的数据。然而，在两级嵌套结构的纵向设计中，来自不同研究个体的观测值之间保持独立，但来自同一研究个体重复观测所得到的观测值之间，通常由于共享相同个体水平的特征而彼此相互关联，因此同一研究个体内的观测值之间不再独立。纵向研究设计得到的重复测量数据不适合使用一般线性模型进行分析，它要求对一般线性模型加以改进以适应纵向数据的特殊性。为此，在生物统计学领域，混合效应模型(mixed-effects model)又称随机效应模型(random-effects model)应运而生[1-3]。自20世纪80年代起，这种模型才开始流行开来，不同学科对该模型的命名也有所不同，例如社会科学领域的多水平线性模型(multilevel linear models)[4-5]、 教育学领域的分级线性模型(hierarchical linear models)[6]、经济学领域的随机系数回归模型(random-coefficient regression models)[7-8]、统计学领域的协方差成份模型(covariance components models)[9-10]。该模型应用到纵向数据上又称为成长曲线模型(growth curve models)[11]。同一个统计方法在不同领域有不同的名称，显示同一种统计方法可以从不同的角度去理解和使用。

1 因变量变化的特征和要求

在现实世界中变化是永恒的，而静止是暂时的，不变可以看作是变化的一种特殊形态。对因变量随时间推移产生的变化进行观测时要注意3点[3]：①被衡量特征的数值应随时间变化而变化，而且该特征需满足两个条件。一是该特征的变化在理论上或实际中有意义；二是该特征的含义在不同时间衡量时都是相同的，即特征的含义不随时间变化而变化。例如观测体重指数(body mass index, BMI)随时间推移发生的变化，无论何时，体重指数的含义始终都是一样的，都等于体重除以平方身高。又如观察生活满意度随时间推移发生的变化，如果受试者对生活满意的定义和条件依时间推移而变化，诸如20世纪80年代拥有老三大件即冰箱、彩电和洗衣机的话，大多数中国人就觉得满意了；20世纪90年代变为拥有新三大件即空调、电脑和录像机才满意。在这种情况下，就无法确定所观测到的生活满意度变化是由于时间变化带来的还是受试者对生活达到满意的要求变化了而导致的。②重复观测或随访的间隔时间需要考虑其适宜性和敏感性。对变化频繁的特征，随访的间隔时间太长就有失去观测到变化的可能，反之对变化不频繁的特征，随访的间隔时间太短就可能观测不到实际变化。例如血压和血糖变化较频繁，以小时为单位比用周为单位更合适。观测抑郁症的变化以抑郁症发作的间歇期为单位比用天或周为单位更合适。由于时间变量是成长曲线模型中不可或缺的基本变量，研究者应在仔细考量之后确定能可靠地、准确地、敏感地、有意义地观测出变化特征的间隔时间。③重复观测的次数应在3或4次以上。对于连续型变量重复观测的次数至少3次，对于二分类变量重复观测的次数至少4次，以提供足够信息来估计模型中的参数。

2 纵向数据的储存格式

纵向数据分析的第一步是将收集的数据以一个合适的格式储存在数据库里。储存纵向数据的数据库一般有两种格式：一种是宽格式(wide format)，即每一个个体只有一行记录其观测的数据，每次重复或随访的观测值均按不同的变量记录在数据库的不同列中。这种格式也可称为个体级别格式(person-level format)；另一种是长格式(long format)，即每一行记录每一个个体每一次随访的观测值，每一个个体因有多次随访数据故有多行记录，每一个个体都有一个特定的个体编号加以识别。对于不随时间变化的个体特征变量(time-invariant variables)，例如性别和民族等，则在同一个个体的每一行重复同样的观测值。对于随时间变化的个体特征变量(time-varying variables)，例如体重和血压等，则在同一个个体的每一行记录不同的观测值。这种格式也可称为个体时间格式(person-period format)。表1显示用长格式记录纵向数据的例子，可以看出每一个个体都有一个编号，重复观测的间隔时间按年计算，年龄和耐受性是随时间变化的个体特征变量，故每人每次的值可能不同。性别和暴露因素是不随时间变化的个体特征变量，故每人每次的值不变。大多数分析纵向数据的软件要求使用长格式的数据储存格式。

表1 长格式纵向数据示例

3 成长曲线模型的构成

成长曲线模型有两种等同的数学表达形式，每种形式各有其优点，没有对错之分。第一种形式是按不同水平的顺序以等级模型(hierarchical model)形式表达。这种形式通俗易懂，但所用方程的数量较多。最低水平是每次重复的时间，上一水平是被研究的个体。下一等级嵌套在上一等级的结构之中，故此模型为最简单的多水平模型(multilevel model)。统计学文献习惯用希腊字母代表需估计的参数，用大写的英文字母代表随机变量，小写的英文字母代表变量取值，本文亦遵循这一传统。

第一水平模型：每个个体重复观测值随时间变化而变化，又称个体内变化(within-person change)

Yij=π0i+π1iTIMEij+π2iXij+εij

(1)

第二水平模型：个体间变化(between-person change)

π0i=γ00+γ01Zi+ζ0i

(2)

π1i=γ10+γ11Zi+ζ1i

(3)

方程(2)中π0i代表第i个个体的截距。γ00代表样本所属群体的截距的平均值。γ01是第二水平自变量Z的斜率，代表自变量Z每增加一个单位时所对应的截距平均变化值。自变量Z反映的是不随时间变化的个体特征，例如社会人口特征等性质较稳定的一类自变量。ζ0i代表截距的残差，即自变量Z无法解释的截距部分。方程(3)中π1i代表时间变量TIME的斜率。γ10代表样本所属群体的时间变量TIME的斜率的平均值。γ11是自变量Z的斜率，代表自变量Z每增加一个单位时所对应的斜率变化，即成长率的平均变化值。ζ1i代表时间变量TIME的斜率的残差，即自变量Z无法解释的成长率部分。总之，每个个体有各自的截距和时间变量斜率，故共有i个截距和i个时间变量斜率。但从统计学角度看，这些不同的截距和时间变量斜率可以被看作是来自一个共同群体的抽样样本，这些截距和时间变量斜率的残差符合多元正态分布(multivariate normal distribution)，

其中τ代表方差/协方差。注意尽管方程(1)中Xij变量的观测值依时间不同而不同，如果Xij变量与Yij之间的关系在每个个体都一样，则π2i是一个常数，即代表固定效应回归系数，因此π2i没有对应的第二水平方程来估计其方差。如果Xij变量与Yij之间的关系在每个个体都不一样，则π2i是一个随机数，代表随机效应回归系数，因此π2i也可像π1i一样有对应的第二水平方程来估计其方差。由于本文重点介绍成长系数，主要关注作为随机效应的时间变量，为了简洁易懂，暂且将π2i作为固定效应回归系数。将方程(2)和(3)代入方程(1)可得到一个复合模型(composite model)，即方程(4)。

Yij=(γ00+γ01Zi+ζ0i)+(γ10+γ11Zi+ζ1i)TIMEij+π2iXij+εij

(4)

方程(4)是成长曲线模型第二种表达形式，它的优点是可以直观地区别出随机效应和固定效应参数，缺点是方程较为复杂。方程(4)的第一个括号中是复合截距，复合截距里包含截距的残差，表明复合截距不是一个常数，而是服从一个分布的复合变量，该变量描述了不同个体的Yij值在起始时间点的差异。复合截距里还包含了自变量Z，表示复合截距受自变量Z的影响，即第二水平的自变量影响第一水平的因变量。第二个括号中是时间变量TIME的复合斜率，复合斜率里包含斜率的残差，表明不同个体在成长速度上不是一个常数，而是服从一个分布的复合变量。复合斜率里也包含了第二水平的Z自变量，表示复合斜率也受Z变量的影响。将上述复合方程重新组合排列以后得到方程(5)，其参数可分为结构部分(structural component)和随机部分(stochastic component)两类。

Yij=γ00+γ10TIMEij+γ01Zi+γ11TIMEijZi+π2iXij+(ζ0i+ζ1iTIMEij+εij)

(5)

方程(5)的括号外部分为结构部分，括号内部分为随机部分。这也是混合效应模型名字的由来。在这个方程中可以直观地看出γ11系数代表了跨水平变量之间的交互效应，即时间变量TIME与Y之间的关系受Z变量的影响。t检验用于对结构部分的参数(即γ00、γ10、γ01、γ11、π2i)是否分别为零进行假设检验。方程的复合残差包含两种独立的残差来源且随时间变化而变化，表明各组内观测值之间不独立，同时也包含异方差(heteroscedasticity)，这是成长曲线模型与经典的一般线性回归方程的根本区别所在。Z检验用于对随机部分的参数(即ζ0i的方差τ00、ζ1i的方差τ11、ζ0i和ζ1i之间的协方差τ01)是否分别为零进行假设检验。成长曲线模型在实际运用过程中会根据具体的研究问题或假设检验的需要加以增减变化。目前主要统计软件包给出的模型结果对应方程(5)。本文之所以对上述两种模型的数学表达形式均作介绍是方便读者理解成长曲线模型。

4 成长曲线模型参数估计方法

方程中参数的常用估计方法有以下两种。

第一种是最大似然估计(maximum likelihood estimation)方法，估计混合效应模型的未知参数[2,10,12]。1970年代后随着期望-最大(expectation-maximization, EM)算法的诞生[13]，通过优化似然函数找到参数的无偏估计被广泛应用。该方法的基本思路是先确定包含参数的似然函数，再给方程中截距、斜率、第一和第二水平的方差各自赋予一任意初始值，得到一个可以观测到实际因变量值的似然值，然后不断尝试其他参数值的组合，并与之前得到的似然值对比，看是否优化，如此循环往复，直到得到可以观测到实际因变量值的最大似然值。该方法在样本足够大时一致性和有效性均佳。但也发现有收敛慢的缺点[14]。随后又有基于Newton-Raphson和Fisher scoring的估计方法诞生[15-16]。这些估计方法已经被各种软件采用。最大似然估计法有一大优点是如果缺失数据的概率与因变量无关，即符合随机缺失(missing at random, MAR)机制，即使有部分资料缺失，最大似然估计法仍然会给出无偏参数估计值，因此成长曲线模型允许部分样本存在部分缺失的随访资料。

最大似然估计法在实际应用中又分为两种情况。其一是对全模型即对结构部分和随机部分一起进行运算，称为全信息最大似然估计法(full information maximum likelihood estimation, FIML)；其二是仅对随机部分即复合残差部分进行运算，称为限制性最大似然估计法(restricted maximum likelihood estimation, REML)。当需要比较一系列模型并用最大似然值或其演化出来的统计指标评估哪一个模型拟合更好时，所比较的模型只有在使用同一最大似然估计方法的条件下才具有可比性，即要么比较用FIML得出的似然值，要么比较用REML得出的似然值。

第二种是广义最小二乘法(generalized least squares, GLS)[17]。这种方法是对一般线性回归模型中使用的最小二乘法的扩展。第一步忽略样本的等级结构，用普通最小二乘法(ordinary least squares)估计其残差的方差/协方差矩阵，假定从样本得到的方差/协方差矩阵就是估计的总体的方差/协方差矩阵。第二步根据估计的方差/协方差矩阵再估计方程中的固定效应参数及相应的标准误。再将第二步得到的新参数放回到第一步重新估计残差的方差/协方差矩阵，然后再回到第二步，根据新估计的方差/协方差矩阵再重新估计方程中的固定效应参数及相应的标准误，如此重复，直到基于预先确定的标准，使得估计的参数无法得到改进为止，专业术语上称之为收敛(convergence)。

除上述两种方法外，估计混合效应模型参数的方法还有很多，并且不断有新方法出现。这些方法已经编入统计软件中以便使用，故对于广大应用者，只需对其估计方法有一个基本的了解即可。

5 成长曲线模型的应用

当研究者收集完数据并以长格式的方式录入到数据库后，就可以开始分析。由于成长曲线模型结构较复杂，一般主张由简到繁，分步骤检验不同的成长曲线模型。本文以一项纵向研究为例，展示如何应用成长曲线模型[11]。

美国青少年至成年人健康纵向研究(the national longitudinal study of adolescent to adult health，简称Add Health)，自1994-1995学年起对一组超过两万名具有全美代表性的七至十二年级的青少年进行跟踪随访，观测这些青少年的发育、成长、成熟、老龄化的全生命过程。该研究收集的数据包括人口资料、社会因素、家庭因素、社会经济状况、行为、心理、认知、生理和遗传、用药。除了调查在校青少年外，它还调查了青少年的家长、学校、社区以及地理环境数据。该研究在1994-1995学年、1996-1997学年、2001-2002学年、2008-2009年、2016-2018年共进行了5次调查。除了可能暴露研究对象身份的敏感数据(包括受访者的生物遗传信息、居住社区及地理信息)需提出申请并经过审查后才能获取使用外，该研究的绝大部分数据已公布在互联网上供研究和教学之用(https://addhealth.cpc.unc.edu/data)。至今已有数以千计使用该数据库的论文发表。

本文选取该研究前四次调研的一部分数据(n=6291)作为学习探讨使用成长曲线模型的实例。20世纪90年代至今的30年期间，美国面临的一个公共卫生危机就是青少年肥胖率迅猛增长。由于肥胖症已被证明是导致心脑血管疾病、肿瘤、精神疾病等健康问题的重要成因，并且很多成年人疾病起因于青少年时期，因此从青少年开始研究肥胖症的发生和发展规律，对于慢性病和精神疾病预防有着重大意义。本文拟从以下两个方面来研究肥胖症：①从1994至2008年，美国青少年体重总体上变化的趋势是怎样的？②是否所有青少年体重变化的趋势都是一致的？哪些因素可以解释青少年体重变化的趋势？体重指数常用来衡量肥胖程度。从表2可以看出在本文使用的样本中，从少年到青年，体重指数从初始时的正常到后来的超重呈现增长的趋势，标准差随时间推移而变大，男女变化的趋势一致。下面本文分步分析数据以回答上述问题，所有结果均使用SAS软件包获得[18]。感兴趣的读者可通过微信公众号(医学与社会)、中国知网、万方等渠道下载本文的附件文件，获取相关数据和SAS程序编码。

表2 四次调查样本的体重指数描述性统计结果

5.1 体重指数变化趋势的图形展示

图1为运用惩罚B样条曲线(penalized B-spline curves)拟合方法显示随机抽取的12位青少年每人在4次随访中体重指数的变化趋势。横轴表示随访时间，纵轴表示体重指数的范围。这12条曲线的形状各异，每次随访时体重指数数值均不同。由于存在失访，有的受访者只有2次或3次体重指数记录。因此，如果仅看每个个体，很难把握体重指数的变化趋势。

图1 12位青少年在4次随访中体重指数的变化趋势(惩罚B样条曲线)

为了更好地看清每个个体的体重指数变化趋势，改用线性回归方程分析每位受访者体重指数随时间变化的线性趋势。图2显示这12位青少年中大部分人体重指数呈直线增长趋势，也有下降和平缓趋势，这12条直线的斜率及每人的初始体重指数均不同。

图2 12位青少年在4次随访中体重指数的变化趋势(线性回归方程)

当随机抽取的样本量增加到100，就有100条体重指数变化直线，如果叠加在一张图中，就很难用肉眼看清这100人体重指数变化的趋势了，如果样本更多难度则更大。研究者可以借助线性回归方程，对所有个体样本的斜率加以平均，归纳出所有受试者的平均变化趋势。图3分别展示了各100位男性和女性青少年体重指数随时间的变化趋势。图3中的左图显示男性结果，右图显示女性结果，加粗的直线代表男性和女性青少年各自的平均线性回归方程，其斜率分别为男性和女性青少年学生体重指数的平均成长率。

图3 100位男性(GENDER=1)和女性(GENDER=2)青少年在4次随访中体重指数的变化趋势(线性回归方程)

从图3中我们可以看到，左图男性青少年的平均体重指数在初始时比右图的女性青少年稍高，14年后到第4次观测时，女性的平均体重指数已经接近男性的平均体重指数，提示男性和女性青少年体重指数状况有别。具体差别可以用成长曲线模型对所有样本体重指数在这14年中的变化规律进行研究。

5.2 无条件均值模型

首先需要检验在面板数据中是否有证据支持引入混合效应模型。如果把方程(1)和(2)简化为没有任何自变量的方程，即：

BMIij=π0i+εij

(6)

π0i=γ00+ζ0i

(7)

BMIij=γ00+ζ0i+εij

(8)

(9)

ICC的值介于0～1之间，其值越接近1，表示使用混合效应模型的必要性越大；其值等于零，则提示没有必要使用混合效应模型。ICC的大小及显著性检验的结果提供了一个是否需要使用混合效应模型的统计判断。方程(8)的主要参数估计结果见表3。

表3中截距显示所有样本在所有时间的平均体重指数为24.73，略小于25.00即在正常范围。第一水平方程残差的方差为14.74，Z检验结果显著(14.74/0.17=86.71，P<0.0001)，表示拒绝该方差等于零的原假设，说明每个学生在不同时间的体重指数观测值围绕个体体重指数平均值波动，个体内差异有显著性。第二水平方程残差的方差为21.32，Z检验结果显著(21.32/0.47=45.36，P<0.0001)，表示学生个体间的截距的差异有显著性，说明个体存在组间异质性(heterogeneity)，组内个体存在聚集现象，即观测值不独立。将表3中方差值代入方程(9)算出ICC值等于21.32/(14.74+21.32)=0.59，它说明学生之间体重指数的差异，即组间的变异可以解释59%的总体体重指数变异，还有41%的体重指数变异是由于时间不同而导致。此结果说明在研究体重指数与时间的关系时，需要把个体内(within-person)因时间不同而导致的体重指数变异和个体之间(between-person)体重指数的变异结合起来考虑，而混合效应模型是目前处理这种两水平数据最佳的统计方法之一，可以将这两种变异均纳入分析。

表3 无条件均值模型结果

表3中还列出了三项评估模型与样本资料匹配程度的指标，即模型拟合度。其中Deviance为偏差，AIC是Akaike's information criterion的缩写，BIC是Bayesian information criterion的缩写。这些指标可用于模型之间的比较，即将无条件均值模型与后续其它演化模型进行比较(model comparison)，模型拟合度指标值较小的模型更好[20-21]。尽管无条件均值模型提供了重要的基本信息，但该模型仅提示体重指数会随着个体的不同而不同，同一个体体重指数也会随着时间的不同而不同，但体重指数与时间的具体关系并未给出。我们需要进一步拟合其它模型，并与无条件均值模型进行比较。

5.3 无条件成长曲线模型

BMIij=π0i+π1iTIMEij+εij

(10)

方程(10)在无条件均值模型的基础上加入了时间自变量TIME。方程(10)为第一水平模型，描述了不同个体在不同时间的体重指数与时间的关系。在本研究中，每次随访的间隔时间不等，第一次与第二次之间相隔2年，第二次与第三次之间相隔5年，第三次与第四次之间相隔7年。间隔时间不相等不会影响混合效应模型的应用，混合效应模型可以灵活处理间隔时间相等或不相等的数据。在分析中对时间变量的编码有两种方法。一是按调查的时间顺序依次编码，即第一次记为0，第二次记为1，第三次记为2，第四次记为3。这种方法可以用来估计体重指数随调查时间的推移变化的轨迹。该编码方法既可用于间隔时间相等，也可用于间隔时间不等的纵向数据。另一种编码方法是按调查时间的实际间隔编码，以体现间隔时间不等。本例可以1994-1995学年为初始点，第一次记为0，第二次记为2，第三次年为7，第四次记为14。该方法估计出每一年体重指数的变化轨迹。两种编码方法估计的变化轨迹的具体数值不一样，但总体方向在多数情况下是一致的。本文采用第二种时间编码方法，因其更接近实际情况。

π0i=γ00+ζ0i

(7)

π1i=γ10+ζ1i

(11)

方程(7)与(11)均为第二水平模型。方程(11)的右端包含了时间变量的总体斜率和个体斜率与总体斜率之间的残差。

将方程(7)和(11)代入方程(10)后得到方程(12)

BMIij=γ00+γ10TIMEij+(ζ0i+ζ1i*TIMEij+εij)

(12)

方程(12)即为无条件成长曲线模型，显示体重指数与时间的关系，括号中是包含有时间变量的复合残差，它显示复合残差会随时间的变化而变化，即存在异方差。方程(12)的主要参数估计结果见表4。

表4中截距估计值表示在1994-1995学年美国青少年的平均体重指数为22.52，小于25.0即正常。时间斜率估计值为正，表明美国青少年的体重指数以每年平均0.42的速度增长。截距(22.52/0.06=375.33，P<0.0001)和时间斜率(0.42/0.01=42.00，P<0.0001)的t检验结果均显著，表示拒绝截距和时间斜率等于零的原假设。表4中第一水平方程残差的方差为4.23，与表3中相对应的第一水平方程残差的方差14.74相比下降了71%，表示71%的个体内体重指数的差异可以被时间变量所解释。由于该方差Z检验结果显著(4.23/0.06=70.50，P<0.0001)，表示不接受该方差等于零的原假设，提示还有其它第一水平的自变量可以引入方程来解释剩余29%的个体内体重指数的变异。

表4 无条件成长曲线模型结果

对于随时间变化的自变量，其变化可分为频繁变化和相对稳定两种。例如许多生理指标，包括血压、脉搏、血脂、血糖、激素水平等，变化频繁，其变化可以分钟或小时为单位。这类因素一般都作为随时间变化的变量处理，即每一次随访都记录下当时的数值，在数据库里存储的每一条记录可能相同也可能不同。又如生活方式和某些行为特征，随时间变化而变化，但这些特征变化不是很频繁，其变化多以年甚至5年为单位计算。如果在一段相对较短的时间内反复观测，这些特征可能没有变化。在这种情况下，这些特征可视为不随时间变化的变量，属于第二水平的个体变量。但是如果随访时间较长，则这类变量可视为随时间变化的变量，属于第一水平的变量。总之，随时间变化的自变量和不随时间变化的自变量的区分不是绝对的，需依具体情况做出最合适的判断。

第二水平初始值的方差Z检验结果显著(4.16/0.04=104.00，P<0.0001)，表示不接受该方差等于零的原假设(P<0.0001)，说明不同的学生个体在研究开始时体重指数有明显的差异。第二水平成长率即时间斜率的方差Z检验结果显著(0.29/0.01=29.00，P<0.0001)，表示不接受该方差等于零的原假设，说明不同的学生个体每年体重指数增长速度明显不同。协方差表示截距和斜率残差间的相关度，表4中协方差估计值为0.08且Z检验有显著性，表示初始体重指数水平与体重指数增长速度之间呈正相关，即初始体重指数大的学生，后续体重指数的成长率越快。假如协方差为负值且Z检验有显著性，则结论恰好相反。值得指出的是根据实际情况，初始值和成长率之间的协方差可能为零也可能不为零，可以通过检验加以验证。无条件成长曲线模型的模型拟合度指标与无条件均值模型相比具有较小的指标值，表明无条件成长曲线模型对数据的拟合更好。

图4显示由方程(12)得出的预测的体重指数与时间的平均线性增长关系，成长率方差0.29描绘了其变异度，即不同的人有不同的增长率。横轴代表时间，单位为年；纵轴代表模型预测的体重指数。

图4 方程(12)预测的体重指数与时间的平均线性增长关系

尽管无条件成长曲线模型显示美国青少年初始的体重指数水平及后来体重指数增长速度因人而异，但是没有揭示哪些因素影响学生个体之间体重指数的差异，需要进一步用其他模型加以探讨。

5.4 加入第二水平性别变量的条件成长曲线模型

众所周知男性与女性在生理、心理、行为、生活方式上都有差别，因此首先检验性别(GENDER)是否影响体重指数。性别不随时间变化，属于第二水平的自变量。在本研究数据录入的长格式中，同一个青少年有多至四条观测记录，故每个青少年的性别有多至四条重复的记录，表明性别不随时间变化而变化。其他类似的第二水平的自变量有籍贯、民族、家庭出身、出生地、血型、出生体重、基线调查的背景特征等。

BMIij=π0i+π1iTIMEij+εij

(10)

与无条件成长曲线模型一样，第一水平模型仍然是方程(10)，显示体重指数与时间之间的关系。

π0i=γ00+γ01GENDERi+ζ0i

(13)

π1i=γ10+γ11GENDERi+ζ1i

(14)

方程(13)表示截距与性别之间的关系，方程(14)表示时间斜率与性别之间的关系，方程(13)和(14)为第二水平模型。将方程(13)和(14)代入方程(10)后得到方程(15)。

BMIij=γ00+γ01GENDERi+γ10TIMEij+γ11TIMEij*GENDERi+(ζ0i+ζ1i*TIMEij+εij)

(15)

方程(15)在一个方程中清楚地展示了分属两个水平的自变量如何统一对第一水平因变量体重指数产生影响。时间和性别的主效应分别显示第一水平和第二水平的自变量对体重指数的影响，时间和性别的交互效应为跨水平变量之间的交互效应，表明时间对体重指数的效应可能依性别不同而不同。方程(15)的主要参数估计结果见表5。

表5 加入第二水平性别变量的条件成长曲线模型结果

表5中截距估计值表示在1994-1995学年，美国女性青少年的平均体重指数为22.39。性别的主效应斜率表明男性比女性在1994-1995学年时的平均体重指数高0.27，t检验结果显著(0.27/0.11=2.45，P=0.017)，表示不接受该斜率等于零的原假设。成长率的主效应估计值表明美国女性青少年的体重指数每年平均增长0.44。t检验结果显示女性初始值(22.39/0.08=279.88，P<0.0001)及其成长率(0.44/0.01=44.00，P<0.0001)的检验均显著，表示不接受女性青少年的初始值及其成长率等于零的原假设。因为时间与性别交互效应的系数为-0.03，所以男性青少年的体重指数每年平均增长率为0.44-0.03=0.41，较女生慢，其t检验结果显著(-0.03/0.01=-3.01，P=0.003)，表示拒绝时间与性别无交互效应的原假设。图5显示由方程(15)得出的预测的男性和女性青少年体重指数与时间的不同的平均线性增长关系。横轴代表时间，单位为年；纵轴代表模型预测的体重指数。由图5可见女性青少年在初始1994-1995学年时平均体重指数比男性青少年低，但体重指数的增长速度较男性青少年快，到2008-2009学年时女性青少年平均体重指数比男性青少年高。

图5 方程(15)按男性(GENDER=1)和女性(GENDER=2)分别预测的体重指数与时间的平均线性增长关系

表5中第一水平的方差Z检验结果显著(4.09/0.06=68.17，P<0.0001)，表示不接受该方差等于零的原假设，说明每个青少年在不同时间的体重指数观测值围绕个体体重指数平均值波动，个体内差异有显著性，再次提示还有其它第一水平随时间变化的自变量可以引入方程来解释剩余的个体内体重指数的变异。第二水平截距即初始值的方差Z检验结果显著(18.41/0.37=49.76，P<0.0001)，表示不接受该方差等于零的原假设，说明不同的青少年个体在研究开始时体重指数在考虑性别的影响后仍有明显的差异，与表4中对应的第二水平截距即初始值的方差相比，数值仅减少了0.02，提示需进一步探索除性别以外的影响因素。第二水平成长率的方差Z检验结果显著(0.10/0.01=10.00，P<0.0001)，表示不接受该方差等于零的原假设，说明在考虑性别的影响后，不同的青少年个体每年体重指数增长速度还有明显的差异，与表4中对应的第二水平成长率的方差相比，数值没有变化，提示需进一步探索除性别以外还有哪些因素影响成长率。协方差为0.08，表示初始体重指数越大的青少年，其体重指数增长速度越快。

模型拟合度指标与无条件成长曲线模型相应指标相比除BIC略有增高外其它两项有明显降低，提示加入第二水平性别变量的条件成长曲线模型比无条件成长曲线模型对数据的拟合更好。在后续模型中其它可能的第二水平变量被继续加入进行拟合。

5.5 加入其它第二水平变量的条件成长曲线模型

前述加入第二水平性别变量的条件成长曲线模型的结果提示除性别以外还有其它因素影响初始体重指数以及体重指数的增长速度，为此本文进一步探讨了随访对象在中学时的交友频率(CONTACT)和就医情况(TREATMENT)对体重指数初始值和成长率的影响。交友频率问题为“在过去的一周，你有多少次与朋友在一起？”选项包括完全没有、一到两次、三到四次、五次及以上。就医问题为“在过去的一年里，是否有你认为你应该就医却没有就医的情况发生？”选项包括有、没有。对于回归模型中自变量的赋值编码，一般建议最小值从零开始，例如对就医问题，赋值为0=没有和1=有；对交友频率问题，赋值为0=完全没有、1=一到两次、2=三到四次、3=五次及以上，这样使得截距的估计值有实际意义，以利于解释。这样赋值使得在包括交友频率自变量的条件成长曲线模型中，截距表示在过去的一周完全没有与朋友在一起的青少年的平均体重指数。如果将交友频率选项赋值成最小为1或其它大于零的整数，则估计出的截距仍表明交友频率选项为零时的平均体重指数，由于交友频率选项的赋值范围不包括零，此时截距就没有实际意义。赋值为零的选项一般称为参照组，任何有意义的应答选项均可以被指定为参照组。例如交友频率的参照组可以指定是一到两次的选项，将其赋值为0、完全没有项赋值为-1、三到四次项赋值为1、五次及以上项赋值为2，则此时截距表示在过去的一周有一到两次与朋友在一起的青少年的平均体重指数。如果第一水平或第二水平自变量是连续型变量，则可以将其最小值、平均值、最大值或其它有实际意义的值定为零，使得对应的截距有实际意义。这种方法在数据分析中称为对变量的中心化(centering)处理。例如本研究中的随访对象为青少年，其年龄都在十岁以上。如果成长曲线模型中将年龄作为自变量并直接采用调查数据中的实际年龄，则所得截距表示1994-1995学年青少年年龄为零时的体重指数，这显然与本研究的对象不相符。若将在1994-1995学年调查对象中的最低年龄或平均年龄定为参照年龄，即每个青少年的年龄都减去该最低年龄或平均年龄，进行中心化处理后再纳入模型，则所得截距表示1994-1995学年时，青少年年龄为最低年龄或平均年龄时的平均体重指数。具体如何进行中心化处理，在处理具体数据时应具体情况具体分析。

本文用成长曲线模型分析了交友频率与就医因素及性别对研究开始时体重指数以及每年体重指数增长速度的影响。与无条件成长曲线模型一样，第一水平模型仍然是方程(10)，显示体重指数与时间之间的关系。

BMIij=π0i+π1iTIMEij+εij

(10)

在第二水平方程中，方程(16)描述了性别、交友频率、就医与截距之间的关系；方程(17)描述了性别、交友频率、就医与时间斜率之间的关系。

π0i=γ00+γ01GENDERi+γ02CONTACTi+γ03TREATMENTi+ζ0i

(16)

π1i=γ10+γ11GENDERi+γ12CONTACTi+γ13TREATMENTi+ζ1i

(17)

将方程(16)和(17)代入方程(10)后得到方程(18)。时间交互效应表明跨水平变量性别、交友、就医分别对体重指数增长率的影响。

BMIij=γ00+γ01GENDERi+γ10TIMEij+γ02CONTACTi+γ03TREATMENTi+γ11TIMEij*GENDERi+γ12TIMEij*CONTACTi+γ13TIMEij*TREATMENTi+(ζ0i+ζ1i*TIMEij+εij)

(18)

分析后结果表明交友频率与就医对初始体重指数水平有影响，但是对每年体重指数增长速度无显著影响，因为γ12和γ13的t检验结果均不显著即接受这两个系数与零无差异的原假设(P值均大于0.05)，故方程(18)可以简化成方程(19)。

BMIij=γ00+γ01GENDERi+γ10TIMEij+γ02CONTACTi+γ03TREATMENTi+γ11TIME*GENDERi+(ζ0i+ζ1iTIMEij+εij)

(19)

方程(19)中有显著意义的主要参数估计结果见表6。

表6中截距估计值表示在1994-1995学年，在过去的一周完全没有与朋友在一起以及在过去的一年里没有发生应该就医却没有就医的美国女性青少年的平均体重指数为22.65。时间斜率估计值显示在调整了性别、交友频率、就医因素影响后的美国女性青少年的体重指数每年平均增长0.44。截距(22.65/0.14=162.14，P<0.0001)和时间斜率(0.44/0.01=44.00，P<0.0001)的t检验结果均显著，表示截距和时间斜率均拒绝了等于零的原假设。性别的斜率表明当交友频率与就医情况相同的条件下，男性青少年比女性青少年在1994-1995学年时的平均体重指数高0.30，t检验结果显著(0.30/0.11=2.72，P=0.008)，表示该斜率等于零的原假设被拒绝。交友频率的斜率估计值是负数，表明在其它因素的影响保持相同的情况下，在1994-1995学年，与朋友一周内接触越频繁的学生其平均体重指数越低。就医频率的斜率估计值显示在其它因素的影响保持相同的情况下，在1994-1995学年，应就医而未就医的青少年比应就医而就医的青少年平均体重指数高0.47。时间与性别交互效应的系数仍然为-0.03，显示在其它因素的影响保持相同的情况下，男性青少年每年的体重指数增长率为0.44-0.03=0.41，较女生慢，t检验结果显著(-0.03/0.01=-3.00，P=0.002)，表示时间与性别交互效应等于零的原假设被拒绝。

表6 加入其它第二水平变量的条件成长曲线模型结果

表6中第一水平的方差Z检验结果显著(4.23/0.06=67.50，P<0.0001)，表示不接受该方差等于零的原假设，说明每个青少年在不同时间的体重指数观测值围绕个体体重指数平均值波动，个体内差异有显著性，再次提示还有其它第一水平随时间变化的自变量可以引入方程来解释剩余的个体内体重指数的变异。第二水平截距即初始值的方差Z检验结果显著(4.14/0.04=103.50，P<0.0001)，表示不接受该方差等于零的原假设，说明性别、交友频率、就医无法完全解释1994-1995学年体重指数的变异，表明还有其它第二水平不随时间变化的自变量可以引入方程。初始值的方差4.14比表5中初始值的方差4.15稍小，说明交友频率和就医解释1994-1995学年时美国青少年平均体重指数的差异的能力非常有限。第二水平成长率方差Z检验结果显著(0.29/0.01=29.00，P<0.0001)，表示该方差不等于零的原假设被拒绝，说明体重指数增长率还有除性别、交友频率、就医以外的其它潜在因素的影响。表6中成长率的均值与方差与表5中相比没有差别，表明还需要寻找其它影响成长率的因素，这是下一步需继续研究的方向。表6中模型拟合度指标与表5中相比具有较小的指标值，表明加入性别、交友频率和就医的成长曲线模型比仅加入性别的条件成长曲线模型对数据的拟合更好，意味着最终的模型更好地概括了数据中体重指数潜藏的变化规律，支持目前的研究方向。

综合方程(8)、(9)、(12)、(19)的结果可发现在美国青少年至成年人健康纵向队列研究中有明显的聚类效应，即青少年个体之间的体重指数差异对其与时间的关系有明显的干扰效应，在用混合效应模型控制了干扰效应后，在调查初始时即1994-1995学年青少年个体体重指数有高有低，其中女性、交友频繁者、应就医而就医者在调查初始时比男性、交友不频繁者、应就医而没就医者的平均体重指数要低。青少年的体重指数从1994年到2008年有很明显的增长，初始时的体重指数差异与后来14年的体重指数增长呈正相关，没有发现初始时的交友频率和就医情况对体重指数增长率有显著影响。女性青少年尽管在开始时平均体重指数比男性青少年低，但后来居上，每年平均增长速度明显快于男性青少年，快速增长的体重指数可能会增加了女性心脑血管疾病、肿瘤、心理疾患的危险性。另外不是所有男性和女性青少年的体重指数增长速度都一样，个体之间差别很明显。下一步的研究需要找到除性别、交友频率、就医情况外还有哪些因素可以影响体重指数增长率，从而为下一步针对加速体重指数增长的因素进行有效干预提供建议，进而降低人群中体重超重和肥胖症的比例，最终达到降低人群中心脑血管疾病、肿瘤、心理疾患的危险性的目标。

6 总结

过去社会科学及医学研究者认为纵向研究设计收集数据不仅费时费事，成本也较高。随着科技的进步，人类收集数据的能力越来越强，对研究质量的要求也越来越高，近二十年来在社会科学和医学领域里使用纵向设计的研究也越来越多。本文概括性地介绍了目前分析纵向数据的一个主要方法即成长曲线模型的基本原理及应用。使用美国青少年至成年人健康纵向调查数据探讨了美国青少年体重指数从1994年到2008年间的变化趋势以及部分潜在因素即性别、交友频率、就医情况对体重指数变化的影响，详细介绍了成长曲线模型中每个参数的实际意义，结合实例重点强调如何最大限度提取参数给出的信息以回答实际问题，涉及了部分相关的数据分析技巧。依次报告了无条件均值模型、无条件成长曲线模型以及条件成长曲线模型的结果，这些中间过渡模型的结果对于问题的探索和理解大有裨益。但在撰写最终报告中可以省略模型的探索过程，只报告最终模型结果。另外，本文在分析美国青少年至成年人健康纵向研究样本时没有考虑复杂抽样设计和样本的加权值对统计检验结果的影响，故分析结果不一定能代表美国青少年的真正体重指数变化情况，上述结论仅供学习讨论成长曲线模型之用。由于篇幅限制，本文仅介绍了普通成长曲线模型，无法顾及其中许多细节及其它演化的模型，例如如何检验模型的线性关系和正态分布的前提假设、非连续性成长率、非线性模型、针对分类和计数因变量的成长曲线模型、meta分析等。有兴趣的读者在对成长曲线模型有了初步了解后可阅读相关的专著以对成长曲线模型有更全面深入的理解，加强合理运用成长曲线模型的能力。

(本文受湖北省卫生技术评估研究中心和《医学与社会》编辑部邀请撰写，编校过程中编辑部张治国和宋芳给予了大力协助，在此一并表示感谢。)