13.2.2.1 数据描述和数据处理
本节将用中国居民健康与营养调查（China Health and Nutrition Survey，CHNS）数据。本章我们将通过CHNS数据观察我国近20年城乡间儿童体重指数的变化。1991年作为基线调查，随后在1993年、1997年、2000年、2004年、2006年、2009年、2011年进行追踪调查，共有8个观察时点。本章使用的数据集名称为exe13，模型构建中所涉及的变量如下。
结局变量：
bmi：体重指数，每次调查均会询问儿童的身高和体重，根据BMI＝体重/身高2公式计算，该变量为连续变量。
水平1解释变量：在发展模型中，水平1解释变量是随时间而变化的研究对象个体内测量，在本研究中包括：
time：调查时间，对1991年、1993年、1997年、2000年、2004年、2006年、2009年、2011年和2015年分别编码为：0、2、6、9、13、15、18、20。
fat：脂肪，每日摄入量，这是一个随时间变化而变化的个体内协变量，为连续变量。
水平2解释变量：在发展模型中，水平2解释变量是不随时间而变化的研究对象个体间测量，数据包括：
region：基线调查时的居住地（1＝城市，0＝农村）
age：连续变量，为基线调查时的年龄。
gender：性别（1＝男，0＝女）；
weight：基线调查时是否超重（1＝超重；0＝不超重）
本节，我们只使用以上变量来介绍和阐述发展模型的概念、特征及如何实际运用发展模型来分析纵向数据。在分析模型时需要长数据格式，即每个研究对象有与各次观察相对应的多个记录。如果是宽数据，我们需要先进行数据格式的转换，本章使用的数据集为长数据，如表13-4所示，故此不必转换。
表13-4 实例分析的数据集（节选）
续表
我们将讨论线性发展模型，从最简单的模型入手，仅在模型中加入一个水平1变量（time），并将水平1截距设定为随机系数。随后，将在此基础模型上逐步扩展并完善该模型。
13.2.2.2 截距测量随时间变化的形式
在最初建立模型时，最好先用图的方式观察一下整个研究期间的结局测量的平均发展趋势（average growth trend）。本节研究的结局测量BMI[1]的平均发展趋势线图（图13-1）显示，研究样本的平均BMI在1991—2011年内有上升趋势，且呈直线线性上升趋势，因此可考虑一次线性时间模型。
13.2.2.3 随机截距发展模型
线性回归模型最简单的扩展是用随机截距替代其固定截距，从而有以下随机截距发展模型（random intercept growth model）。
水平1（个体内）模型为：
图13-1 1991—2011年儿童BMI趋势图
水平2（个体间）模型为：
组合模型为：
个体内模型表示个体j的第i次结局测量值由结局测量初始水平（β0j）和变化率（β1j）决定。个体间模型将个体j的初始结局水平分解成两个部分：总体结局测量平均初始水平（γ00）和个体特定效应（u0j）。通过加入个体特异性（u0j），结局测量的初始水平便因人而异，随个体不同而不同。在此模型中，时间变量time的斜率被设定为固定效应，即结局测量随时间推移的变化率是一固定系数γ10，不因研究对象的不同而变化。换言之，个体研究对象可以从不同的结局水平开始，但其随时间推移的变化率相同，即个体发展趋势与模型估计的总体发展趋势平行。
个体特异性效应u0j是随机效应，其被假设为服从均数为0、方差为的正态分布总体中随机抽取的。误差项eij也被假设服从均数为0、方差为σ2的正态条件独立分布。由于由个体间变异引起的噪音已从未解释的随机残差中去除，所以称为条件独立分布。但是两个不同观察水平的随机残差项u0j和eij是相互独立的。
该随机截距模型与传统线性模型的唯一区别是，模型中的截距不再是固定的系数，而是一个随机系数。模型中水平2残差项u0j代表各个体发展趋势线与总体水平发展趋势线
或平均发展轨迹之间的差异，而该差异的方差
代表所有个体发展趋势与平均发展趋势的离散程度。如果
为0，则u0j＝0，因为E（u0j）＝0。那么就会只有一条发展或轨迹线，即平均发展或轨迹线。
SAS程序：
proc mixed plots(maxpoints＝none)data＝exe13 covtest noclprint；
class idind；
model bmi＝time/s ddfm＝kr；
random int/subject＝idind；
run；
上述SAS程序中，语句random int/subject＝idind要求SAS拟合随机截距模型。
SAS结果：
SAS部分结果输出如下：
图13-2 迭代历史
图13-4 模型拟合
图13-5 固定效应估计值
SAS结果解释：
图13-2显示模型估计经过四次迭代收敛。图13-3显示的是残差方差/协方差估计：
，σ2＝4.1297（p＜0.0001），二者均呈统计显著。
的显著性表明研究对象的结局测量（BMI）的初始水平显著不同；σ2显著性表明在模型设定了随机截距后，仍然存在显著的个体内变异（within-subject variation）。组内相关系数ICC＝4.5464/（4.5464＋4.1297）＝0.52，表明一半的总变异是有研究对象个体间异质性（between-subject heterogeneity）引起的。
由于该模型估计使用的是REML法（SAS PROC MIXED程序的默认估计法）（bell，2013），SAS输出中所示的模型拟合统计量是用来比较带有相同固定效应但不同随机效应或不同残差方差结构的模型，限制性－2LL（restricted－2LL）可用于LR检验来比较该类嵌套模型。信息标准统计量AIC，AICC和BIC可用于嵌套模型或非嵌套模型比较（Richard，2011）。但是－2LL和信息标准统计量都不能回答某个单一模型是否拟合数据，它们仅用于模型比较，说明哪个模型相对来说拟合数据更好。
图13-5显示固定效应的参数估计值。截距
，表示模型估计的结局测量BMI的总体平均初始值为16.856；线性时间效应
，表示结局测量BMI每年的平均变化率为0.3042。
13.2.2.4 随机截距-斜率发展模型
随机截距模型中，结局测量的个体发展趋势线有不同的截距，但有相同的斜率，即研究对象个体的初始结局测量水平不同，但每个个体的结局测量随时间的变化率都相同。但实际研究中，通常情况是，结局测量发展不仅初始水平因人而异，且其随时间的变化率也不尽相同。因此，更为符合实际情况的模型是随机截距和随机斜率发展模型，即随机截距－斜率发展模型（random intercept-slope growth model）。该模型如下。
水平1（个体内）模型：
水平2（个体间）模型为：
组合模型为：
式中u0j代表第j个个体的结局测量初始水平偏离模型估计的总体平均初始结局水平程度，u1j代表第j个个体的结局变化率偏离估计的总体平均结局变化率程度。假设u0j和u1j为二元正态分布（bivariate normal distribution），即N～（0，G），其具有如下方差/协方差矩阵：
其中，
和
分别代表结局测量初始水平和变化率变异的方差，
表示结局测量水平和结局测量变化率之间的协方差。
SAS程序：
以下SAS程序拟合随机截距－斜率发展模型。
proc mixed plots(maxpoints＝none)data＝exe13 covtest noclprint；
class idind；
model bmi＝time/s ddfm＝kr；
random int time/subject＝idind G type＝un；
run；
SAS结果解释：
以上SAS程序的random语句包括了int（代表intercept，截距）和时间变量time，因而模型截距和变量time的斜率均被处理为随机回归系数，即个体发展趋势线的结局测量初始水平和变化率均因人而已。SAS程序的random语句中type＝un选项设定模型的G矩阵包含水平2残差有非结构方差/协方差。选项ddfm＝kr要求SAS用Kenward－Roger法计算固定效应检验的分母自由度。
SAS结果：
SAS部分结果输出如下：
图13-6 协方差参数估计值
图13-7 模型拟合
图13-8 零模型似然比检验
图13-9 固定效应估计值
SAS结果解释：(www.xing528.com)
与随机结局发展模型一样，图13-9显示该模型有两个固定效应：截距
（p＜0.0001），表示模型估计的结局测量BMI的总体平均初始值为16.6989。线性时间效应
（P＜0.0001），表示结局测量BMI每年的平均变化率为0.3053。图13-6显示有3个随机效应：随机截距的方差估计
，表明结局测量的初始水平在研究对象之间有显著差异。时间变量的随机斜率的方差估计
（p＜0.0001），表明结局测量随时间的变化率在研究对象之间有显著差异；该截距和斜率系数间的协方差为
，表明研究对象的BMI初始水平越高，则其结局测量时间的推移的变化率就越小。模型中增加随机斜率后，水平1残差方差从4.1297下降到3.0795，减少了25.43%。图13-7显示该模型与数据的拟合优于随机截距模型，因为－2LL从64783.4下降到63729.8，差值为1053.6，表明拟合改善统计显著（df＝2，p＜0.0001）。
13.2.2.5 城乡效应的评估
为了评估城乡对结局测量的效应，需要将水平1随机截距β0j和斜率β1j在水平2方程中处理为城乡变量（region，1－城市，0－农村）的函数，模型表达如下。
水平1（个体内）模型为：
bmiij＝β0j＋β1j timeij＋eij
水平2（个体间）模型为：
β0j＝γ00＋γ01 regionj＋u0j
β1j＝γ10＋γ11 regionj＋u1j
组合模型为：
模型中γ00为农村的儿童BMI平均初始水平，γ01为城市儿童（region＝1）和对照组（region＝0）之间儿童BMI初始水平的平均差异。γ10代表农村儿童BMI平均变化率，γ01代表农村和城市之间儿童BMI平均变化率的差异。
SAS程序：
运行SAS程序如下：
proc mixed plots(maxpoints＝none)data＝exe13 covtest noclprint；
class idind；
model bmi＝region｜time/s ddfm＝kr；
random int time/subject＝idind G type＝un；
run；
SAS结果解释：
model bmi＝region｜time同时设定变量time和region的主效应和交互作用。等效SAS语句为model bmi＝region time region*time。
SAS结果：
SAS部分输出结果如下：
图13-10 协方差参数估计
图13-11 模型拟合
图13-12 零模型似然比检验
图13-13 固定效应估计值
SAS结果解释：
图13-13显示，
，表示模型估计的农村儿童BMI的平均初始值为16.6451。
，表示城市儿童比农村儿童BMI初始水平平均高0.2097（p＜0.0001）；
，表示农村儿童BMI每年平均变化率为0.3050。
，表明region变量与time变量的交互作用没有统计学差异，交互作用大于0表示随一个单位时间的增加，在城市的儿童，相对于农村儿童而言，他们的BMI平均结局值多增长0.001，但是这个额外的增长没有统计学差异。
图13-11显示，与随机截距－斜率模型结果相比，该模型数据的拟合得到显著改善（－2LL之差为63716.3－63710.1＝6.2，df＝2，p＜0.05），[2]所以须将变量region放入模型。
图13-10显示随机效应仍然统计显著，但加入region变量后，随机截距和斜率的方差基本没有变化，说明该变量不能解释随机截距和随机斜率的变异。
13.2.2.6 在模型中控制个体协变量
实际研究中，为了区分地区效应和其他“噪音”，应将理论上有关的个体变量作为协变量加入模型。本研究中，应加入age，gender，weight等变量，模型如下。
水平1（个体内）模型为：
bmiij＝β0j＋β1j timeij＋eij
水平2（个体间）模型为：
组合模型为：
其中，假设所有的协变量对结局测量的影响都不随时间而变化，也就是说，协变量与时间变量time之间不存在交互作用。
SAS程序：
运行SAS程序如下：
proc mixed plots(maxpoints＝none)data＝exe13 method＝reml covtest noclprint；
class idind；
model bmi＝region｜time age gender weight/s ddfm＝kr；
random int time/subject＝idind G type＝un；
run；
SAS结果：
SAS输出的部分结果如下：
图13-14 协方差参数估计
图13-15 模型拟合
图13-16 零模型似然比检验
图13-17 固定效应估计值
SAS结果解释：
图13-17显示，变量gender的效应为
，表示男孩的BMI平均初始水平比女孩高0.03，但是性别差异并没有统计显著。变量weight的效应为
＝6.1532（p＜0.0001），表示基线为超重的儿童BMI平均初始水平比基线不超重的儿童高6.15。age为连续变量，由于参照组的解释变量取0值，而age＝0无实际意义，所以此模型对age进行了总均数中心化处理（age的均值约为8.8），此处age的效应为
（p＜0.0001），表示相对于儿童在8.8岁时BMI的初始平均值，年龄每增加一岁BMI会增加0.2913。
13.2.2.7 在模型中纳入时间变化协变量
发展模型的优势之一是能够容易地将时间变化协变量纳入模型。在此例数据中，变量fat（每日脂肪摄入量）是一个时间变化测量，其可能在某种程度上影响不同结局测量BMI。
SAS程序：
proc mixed plots(maxpoints＝none)data＝exe13 method＝reml covtest noclprint；
class idind；
model bmi＝region｜time mean＿age gender weight mean＿fat/s ddfm＝kr；
random int time/subject＝idind G type＝un；
run；
在以上SAS程序model语句自动识别时间恒定和时间变化变量。时间变化变量在个体内和个体间均可取不同值（如此例中的time，fat），是水平1变量。时间恒定变量（如gender等）只有在个体间变异，是水平2变量。
SAS结果：
SAS输出部分结果如下：
图13-18 协方差参数估计
图13-19 模型拟合
图13-20 零模型似然比检验
图13-21 固定效应估计值
SAS结果解释：
图13-21显示，变量fat（此处fat变量进行了总均数中心化处理）对结局测量BMI有显著正效应（0.0054，p＝0.0015），表明儿童每日脂肪摄入量越高，其BMI越高。
免责声明：以上内容源自网络，版权归原作者所有，如有侵犯您的原创版权请告知，我们将尽快删除相关内容。
