转自个人微信公众号【Memo_Cleon】的统计学习笔记：重复测量数据分析系列：广义线性混合模型（GLMM）。
广义线性混合模型（GLMM）可以看做是广义线性模型和线性混合模型的融合，可以处理不呈正态也不独立的数据。
示例：某溶栓药物治疗20名急性脑梗死患者的疗效，采用随机、双盲、安慰剂平行对照设计，每组各10例，分别于治疗前及治疗后8周每周进行随访观测，观测指标为神经系统体征评分（MDNS）。示例来源：杨珉.李晓松等.医学和公共卫生研究常用多水平统计模型.北京:北京大学医学出版社,2007.5.
此次笔记只是演示广义线性混合效应模型在重复测量数据上的操作步骤，并不是一个完整的案例分析。示例一开始便建立完整的“全模型”，残差方差协方差结构则设定为不同时间点的纵向资料分析常见的一阶自回归，然后逐渐去除掉没有意义的因素。
用广义线性混合模型来分析这个连续型数据的重复测量的示例，本质上就是用广义线性混合模型（GLMM）来实现多层线性混合模型（LMM）而已。为了更好地理解模型参数代表的意义，我们先重新温习一下多层线性混合模型，本例全模型如下：
加入背景协变量age后，结果会有校正。在组合模型里面，截距γ00是Trtgj=0、ageij=0、timeij=0时结局测量MDNS的平均得分值。这里要特别强调一下“变量取值=0”：在向模型中添加变量时，我们往往把分类变量作为因子(Factor)纳入，而连续变量作为协变量(Covariate)纳入。如果把分类变量作为协变量纳入，则会按连续变量处理，取值为分类变量各水平的赋值，同样如果把连续变量作为因子进行分析，则会把连续变量的取值作为其各个水平的赋值。对于按协变量纳入模型的变量而言，变量取值为0和赋值为0是一致的，如本例中的age和time，age=0表示年龄为0岁，time=0表示治疗前。但是对按因子纳入模型的而言，这里的“初始水平”或者说“取值=0的变量水平”指的是变量被设为参照的那个水平，STATA默认低水平为参照水平而且可以指定任意水平为参照，而SPSS一般默认高水平为参照水平，在广义模型里可以通过顺序排序进行调整，如本例的Trtg，在[构建选项]选项卡中分类预测因子默认是按升序排列的，其取值为0的水平表示Trtg=1（治疗组），γ00就是年龄为0岁(age=0)的治疗组(Trtg=1)的研究对象在治疗前(time=0)的MDNS均值，本例改为按降序排列，其取值为0的水平便表示Trtg=0（对照组），截距γ00是Trtgj=0、ageij=0、timeij=0时结局测量MDNS的平均得分值即年龄为0岁的对照组的研究对象在治疗前的MDNS均值，相应的γ01则是治疗组(Trtgj=1)与对照组（Trtgj=0）的年龄为0的研究对象在治疗前(time=0)MDNS的平均差异。在当前模型中我们假定变量age对结局测量MDNS的影响不随时间变化而变化，即对截距可能产生影响而对time的斜率无影响，γ10是对照组(Trtgj=0)研究对象的MDNS平均变化率，γ11则是治疗组(Trtgj=1)与对照组(Trtgj=0)的MDNS的平均变化率差异。
还有一个问题需要注意，对当前模型而言，连续变量age采用的是原始值，age=0是不存在的，所以截距并没有实际意义，因此一般来说连续型变量需要进行中心化处理，中心化处理之后截距γ00代表的就是age取均值时对照组的研究对象在治疗前的MDNS均值，限于篇幅本示例仅演示广义线性混合模型的操作，并没有age中心化的处理。
对于治疗主效应（模型系数），也需要特别说明一下，这关系到结果的正确解读。在临床研究中，多数研究会采用各种方法（如随机化）让基线值无统计学差异，而且基线值常常是在干预之前，此时干预组和对照组都都没有被施加干预因此两组常常无差异。而模型中干预因素系数实际上是初始水平(time=0)的组间差异，个人理解就是用time=0时的单独效应（上图中的γ01：治疗组(Trtgj=1)与对照组（Trtgj=0）的研究对象在治疗前(time=0)MDNS的平均差异），用其来代表干预因素除去交互作用后的效应。所以如果发现固定效应的检测结果治疗因素无统计学意义不要失望，而是应该欣喜。既然基线差异并不能代表干预的效果，我们可以将治疗终点设为参照水平，用治疗结束时治疗组和对照组之间的差异来代表治疗效应，或许这就是为何SPSS默认高水平为参照水平的原因了。当然这个参照水平我们可以修改，一是利用广义混合线性模型[模型选项]中改变显示估计均值的连续变量值，二是直接修改时间的取值，对时间尺度进行重新编码。当然如果时间是按分类变量纳入，默认的就是高水平为参照水平，这个也可以通过广义线性混合模型[构建选项]里面的顺序排序来修改。
操作步骤：
【1】数据录入：具体略。
【2】广义混合线性模型：分析>>混合模型>>广义线性…
①数据结构：将变量id拖到[主体Subjects]上，将变量time拖到[重复测量]；点击更多，重复协方差类型选择一阶自回归。
②字段和效应
目标：选择因变量MDNS。线性模型的目标分布与关系部分（图中红框部分）可以选择不同的数据类型，可以扩展到正态分布以外的数据类型。
固定效应：将变量age、Trtg、time拖到[主(效应)]上，同时选中变量Trtg和time，拖到[双向]上;
随机效应：点击[添加块]打开随机效应块对话框；将变量time拖到[主要]列表框上，此步是建立变量time的随机斜率，即每个个体的MDNS随时间的变化率不同；选中复选框[包含截距]，主体组合选入变量id，随机效应协方差类型使用默认的方差成分，此步是将变量id设为随机变异的来源，即设定截距在不同的个体间是不同的。同时设定随机截距和随机斜率的协方差结构。
③构建选项：分类变量预测因子按降序排列。本例之所以如此，是因为SPSS默认自变量高水平为参照水平，本例安慰剂和治疗组分别赋值0和1，结果是与治疗组相比，安慰剂如何如何。从逻辑上我们想知道，治疗组比安慰剂组有没有改善，即以安慰剂为参照水平，按降序排列后会达到这个效果。另外自由度的估算方法、固定效应及系数检验方法不同，结果可能会有些微的差别。
④模型选项：可以对分类变量进行边际均数比较。本例Trtg选中成对比较，比较治疗结束时(time=8)采用age均值进行校正的结果，多重比较采用默认的LSD法。
如此处设置age=0，time=0，则会得到固定效应系数检验完全一致的结果。如想更好地理解，可参见前面对多层线性混合模型的释义。
【3】结果解读：结果显示基本的个案处理信息和结果缩略图，可以双击缩略图进入模型浏览器查看详细内容。
①模型概要：输出因变量，概率分布，链接函数及信息准则。在纳入不同数量的自变量或选择不同的方差-协方差结构时，可通过信息准则来判定更优的模型。
②数据结构：列出模型的层次结构。本例高层级有20个水平，每个水平重复观测9次。
③预测值和实测值的散点图：可见预测值和实测值存在较好的正向相关，模型拟合良好。
④固定效应检验：默认以图形样式给出各因素的参数关联强度，粗线表示有统计学意义（P<0.05）的变量，可通过左下角的展示样式下拉框将图形样式切换为表格样式，表格样式中有统计学意义的P值带有黄色背景。可以通过横条上的P值来显示相应条件的自变量。
结果显示：
i）模型有统计学意义（F=115.256，P=0.666>0.05），至少有一个变量的系数不为0。
ii）年龄age的主效应无统计学意义（F=2.684，P=0.103>0.05），对整个模型的不产生影响。
iii）变量Trtg主效应无统计学意义（F=0.187，P<0.001）。再次提示此处检测的是治疗前（time=0）治疗组和对照组在治疗前的差值是否等于0，即治疗组和对照组在治疗前有无统计学意义，一般研究都尽量让组间基线无差异，因此其并不能代表治疗组的治疗作用，可以通过改变time编码赋值或者在[模型选项]中的[估计均值]中进行设置需要比较的时间点。
iv）time主效应有统计学意义（F=340.918，P<0.001），随着时间的改变，MDNS的改变不为0。是递增还是递减呢？可以进一步查看后面固定效应的系数。
v）Trtg与time的交互项有统计学意义（F=108.024,P<0.001），随着时间的延长，治疗组和对照组MDNS的改变幅度是不一样的。本例Trtg与time的交互作用有统计学意义，单独分析Trtg与time的主要效应已无多大的实际意义，而是需要进一步分析单独效应（按某个因素不同水平下另外一个因素的效应），而且无论Trtg与time有无统计学意义都应该纳入到模型中。
⑤固定效应的系数估计和检验结果：同固定效应的估计结果视图，可通过左小角的样式下拉框来选图标和表格显示样式，也可以通过横条上的P值来显示相应条件的自变量。
i）截距=98.348(P<0.001)，截距在当前模型中的含义是年龄为0岁的对照组的研究对象在治疗前的MDNS均值，因未进行年龄的中心化处理，其实际意义不大。如将年龄数据进行中心化处理（年龄与平均年龄差值代替原来的年龄值）后，截距的就变得很有意义。
ii）年龄越大，初始值MDNS也会越大，平均每增加1岁MDNS会增加0.239分，但没有统计学差异（P=0.103）。
iii）Trtg的系数值=-1.145（P=0.666>0.05），表明治疗组的MDNS比对照组低1.145，且效应无统计学意义。正如前面模型系数释义部分所言，该差值是治疗前治疗组和对照组的差值，仅能说明基线无差异。用基线的平均差值并不能很好地代表治疗效果，本例在[模型选项]中将治疗终点设为参照水平，交互作用之外的治疗效应可在估计均值的比较结果中进行查看。
iv）Time系数值=1.227（P<0.001）,表明对照组MDNS随时间变化具有统计学意义，每增加1个时间单位，对照组MDNS增加1.227分。
v）交互作用Trtg*time系数=3.161>0，表示时间每增加一个单位，相比对照组，治疗组中研究对象的MDNS将会有更多的增长（多增长3.161分），且这额外的增长具有统计学意义（P<0.001）。每增加1个时间单位，治疗组MDNS增加(1.227+3.161)分。
另外，本例采用默认的自由度计算方法、固定效应及系数检验方法，最终的结果P值结果会与使用线性混合模型的结果略有差异。以Trtg为例，当前模型（t=-0.432，P=0.666）,而采用线性混合模型的结果会是（t=-0.432，P=0.671），原因就是在步骤③构建选项我们采用了默认的自由度的估算方法和固定效应及系数检验方法，改变自由度的估计方法，会得到完全一致的结果。
⑥协方差矩阵：显示当前模型高水平方差协方差矩阵，即G矩阵。当前模型采用的方差成分结构。
⑦协方差参数估计值：显示随机效应的统计结果。
结果默认显示的是残差【AR1结构】的方差和相邻观察间的相关系数。方差为9.099，相邻观察间相关系数为0.317，两者均有统计学意义。点击左下角效果选择框，选择Block1展示随机效应的结果，截距因人而异，截距方差为21.801，且有统计学意义（Z=2.472，P=0.013<0.05），表明截距是随机的，但斜率在不同个体间的变异无统计学意义（Z=0.125，P>0.05），即当前模型按随机截距处理就可以了，并不需要设为随机斜率。
⑧估计显著效应的均值：显示固定效应检验中有统计学意义的因素的不同取值水平下因变量均值的点估计和可信区间。本例因子Trtg无统计学意义，因此未计算估计均值。
⑨估计的总均数
⑩估计的均数：默认显示协变量在均值时因子变量的估计均数。本例在[模型选项]设置显示连续变量time=8，age=均值时的Trtg估计均数，且进行成对比较。结果显示在治疗结束时，治疗组比对照组的MDNS高24.146分，有明显的的统计学差异（P<0.001）。
⑪拟合模型的设定概要：显示当前模型参数设置情况。
当前模型虽然低层级的斜率随组变化，但是这种变化完全由高层级的变量Trtg来解释（斜率变异无统计学意义），背景协变量age对模型无明显影响，可将模型简化为：
操作步骤基本相同，不同的地方如下:
[字段和效应]的固定效应部分，去掉age；
[字段和效应]的随机效应：点击编辑[块]，进入随机效应构建器，删除time后确定。
主要结果如下，感兴趣的可以自己解读一下：
最后还是要再次说明一下，本次笔记重点演示的是广义线性混合模型的操作及模型参数的解读，并不是一个完整的案例分析，比如多次提及的连续变量数据的中心化，以及模型方差协方差结构的比较、交互作用之后的分层分析等都未进行。
转自个人微信公众号【Memo_Cleon】的统计学习笔记：重复测量数据分析系列：广义线性混合模型（GLMM）。
END
