二、广义估计方程
(一)广义估计方程的思想
广义估计方程的计算过程很复杂，但思想却并不难理解。该方法假定在多次测量之间存在一定的相关结构(广义估计方程中叫做作业相关矩阵)。对于重复测量数据而言，最主要的问题就是存在各次测量之间的相关性，从而不能用常规的线性模型等方法。所以广义估计方程思想很简单，就是把这种相关进行校正一下，然后得到校正后的参数估计值，这样就比较可靠了。
(二)广义估计方程中的作业相关矩阵
由于不同时间点观测之间的相关大小存在各种可能性，因此作业相关矩阵也有多种，常见的包括：
(1)独立结构(independence structure)，即不同时间点 上的测量值之间彼此独立，无相关关系。这种结构因为数据完全独立，实际上也无需考虑广义估计方程，直接采用常规的广义线性模型即可。
(2)等相关结构(exchangeable correlation structure)，即 假定任意两次观测之间的相关性是相等的，不随两个时间 点之间的间隔大小而改变。不管是第1次观测与第2次观测，还是第3次观测与第5次观测，相关系数都相等。
(3)一阶相关结构(one‐dependent structure)，表示某时间点的测量值只与其临近时间点的观测存在相关性，而与其他时间点的观测无关。例如，第2次观测只与第1次和第3次有相关，而与第4次无关。
(4)自相关(autocorrelation)，即相关大小与间隔次数有 关，相邻两次观测之间相关较强，间隔越远，相关性越小。例如，第 2次观测与第 1次和第 3次观测相关性较大，与第4次观测的相关性较小。
(5)无结构相关(unstructured correlation)，即假定不同 时间点观测值的相关系数各不相等，不存在前面几种相关 结构的规律。
作业相关矩阵的选择是广义估计方程中很关键的一部分，需要一定的统计学知识来判断。不少研究认为，作业相关矩阵的选择对参数估计结果的影响不大。然而实际数据分析中，指定不同的作业相关矩阵有时确实会产生不同的 参数估计值和标准误(尽管这种情况很少见)。尽管一般差别不大，但笔者仍建议， 尽量指定最为合适的作业相关矩阵，以获得最可靠的估计结果。
如何选择合适的作业相关矩阵，建议结合以下两种方式综合考虑：
(1)根据不同时间点观测值的相关系数矩阵考虑。简单来说，先计算各次相关系数，大致观察一下相关系数情况，然后进行判断。
如果任意两次的相关系数差不多，可考虑等相关；
如果相关系 数出现随时间间隔而规律性减小的趋势，可考虑自相关；
如果无明显的规律，可考虑无结构相关。
理论上，指定无结构相关最为稳妥，可以满足任意情形的相关系数矩阵，但它需要估计的参数也最多。例如，对于 5次重复测量，如果指定等相关，只需要估计 1个参数即可(只有 1个相关系数)；而无结构相关则需要估计任意两个时间点的相关系数，即 10个参数，估计参数过多容易导致统计学效能(power)的降 低。因此，实际分析中需要综合考虑，根据相关系数矩阵的 提示选择较为合理的作业相关矩阵。
(2)结合QIC指标(quasi‐likelihood under the independence model criterion)选择。QIC类似于广义线性模型的拟合优度指标 AIC，只是最大似然值换成了准似然值。
对QIC不理解也无所谓，关键知道，其值越小表示选择的作业相关矩阵越合适。与AIC指标类似，QIC 指标中也有对变量的惩罚项，即 QIC 值不一定随着模型中 变量的增多而变小，只有模型中含有意义的变量，其值才会变小，提示模型更优；如果纳入无意义的变量，其值反而会 升高，提示模型变差。实际分析时，可以分别指定不同的作业相关矩阵，然后比较各自的QIC值，选择其中较小者。
(三)广义估计方程的用途广义估计方程主要用于重复测量数据的分析，这里的重复测量不仅包括临床试验中较为固定、时间点较少的情形，也包括像生长发育监测、流行病学人群纵向观察等时间点较为灵活或时间点较多的情形。在临床试验的重复测量数据分析中，广义估计方程也可以用于组间比较、时间点的比较、组间趋势变化的分析。在其他纵向观测数据中，广义估计方程可根据研究目的进行灵活分析。
(四)广义估计方程的SAS软件实现
我们仍然采用上一篇文章的数据作为例子。为了方便，我们把上一篇文章的基本数据(表1)和图示(图1)放在下面，免得大家来回翻。
广义估计方程的操作需要先进行一定的探索，确定作业相关矩阵(其实往往很多统计分析都是这样，真正写在文章中的结果都是精华，但其实可能前期我们已经做了非常多的工作，但不可能把所有工作都写在文章里)。
本例中我们分别指定了各种不同的作业相关矩阵，结果均一致，因此本例可任意指定一种作业相关矩 阵，结果不受影响。简单起见，我们指定作业相关矩阵为等相关。
对例 1数据采用基于等相关作业相关矩阵的广义估计方程，首先不加入时间与组别的交互项，先分析时间与组别各自的主效应(主效应是基于所有人 (即不分组)的结果)。SAS程序如下：
data ex2； input id group time y； cards；…… ； proc gee data=ex2；class id time/param=reference ref=first； model y=time group；repeatedsubject=id/within=time type=exch corrw； /*subject 指定个体变量，重复测量数据中通常为个体的id编号；within指定重复测量的变量，通常是时间点变量；type指定作业相关矩阵；corrw指定输出作业相关矩阵*/run；表 4 显示了组别与时间的主效应，结果提示，两组之 间 Y 值评分差异有统计学意义(P=0.002)，治疗后第 3周与 治疗前差异有统计学意义(P=0.005)，治疗后第 4周与治疗 前差异有统计学意义(P<0.001)。
主效应是基于所有人 (即不分组)的结果，因此，表 4 结果对应于上一篇文章重复测量方差分析表3 结果中的总体比较(尽管结果并不完全一致，这很正常)。参数估计值显示了差异情况，例如，group 的参数估计值显 示组间差异为 7.8，即试验组的均值(114.6)与对照组的均 值(106.8)相比高 7.8；time 1 vs 0 的参数估计值显示组间差 异为 1.4，提示第 1 周均值(108.9)比治疗前均值(107.5) 高 1.4。其余time 2 vs 0等的解释以此类推。
如果分析中不加入时间与组别的交互项，相当于假定两条线是平行的，然而实际中这一假定并不一定满足。图1可以看出两条线可能不平行(虽然上一篇文章的重复测量方差分析并不认为两条线不平行，但对于数据分析来说，我们一开始并不知道，都是通过简单图示探索先得到一定认识，然后基于这种认识再深入分析)，因此考虑在分析中纳入时间与组别的交互项，以便观察两组的变化趋势是否有差异。加入交互项的SAS程序如下：
data ex2;
input id group time y;
cards;
……
;
proc gee data=ex2;
class id time/param=reference ref=first;
model y=time group time*group;
/*这里加入了交互项，以反映两条线是否平行*/
repeated subject=id/within=time type=exch corrw;
run;
表5显示了加入组别与时间交互效应的结果。一旦加入交互效应，组别与时间点反映的不再是主效应，而是单独效应(这句话非常关键，一定要牢牢记住)。因此，如果想了解组别与时间点的主效应，可以先不加入交互项。单独效应反映的不是所有人的估计结果，而是某一亚组(如对照组的观测、第1周的观测等)的估计结果。
下面这段结果的解释非常重要，建议一定仔细看。对于想了解交互效应如何解释的朋友，尤为重要。这一段不仅是对广义估计方程的解释，也是对常见其它模型中存在分类变量交互项的解释。
单独效应的结果与变量赋值有很大关系，本例中试验组赋值为1，对照组赋值为0，时间点分别赋值为0~4。因此，表5中group反映的不是所有人两组的差值，而是治疗前这一时间点的两组差值(4.2)；同样，time 1 vs 0反映的也不是所有人在第1周与治疗前的差值，而是对照组第1周与治疗前的差值(1.4)。
交互项的结果对应于重复测量方差分析表3结果中的分组比较。例如，group*time(1 vs 0)的参数估计值为0，它反映了第1周两组差值(4.2)与治疗前两组差值(4.2)的差值，也可以说，反映了试验组第1周-治疗前的值(1.4)与对照组第1周-治疗前的值(1.4)的差值(仔细体会一下这两种说法)，两种说法均可，取决于研究目的侧重说明什么。其他交互项的解释含义以此类推。
(五)广义估计方程分析的注意事项
(1)尽管广义估计方程需要考虑作业相关矩阵的设置，但绝大多数情况下，结果是一致的。建议实际分析中，首先可指定不同的作业相关矩阵，观察分析结果是否一致，如果一致，可以任选其一，否则可根据相关矩阵和QIC综合考虑，选择最合适的作业相关矩阵。
(2)广义估计方程的结果比重复测量方差分析更接近模型的形式，因此不少非统计学专业人员可能对结果的解读存在一定困难，尤其是加入交互项的结果解读，需要仔细体会，否则很容易出现结果的解释错误。
(3)广义估计方程比重复测量方差分析在分析思路上更为灵活，但这同时需要对统计学知识和软件操作的更高要求，因为广义估计方程的结果与自变量赋值有很大关系。例如对时间点赋值0~4，与赋值为1~5，二者给的结果会有不同。这一点其实在所有的模型类都是如此，分类资料的赋值很重要。
(4)广义估计方程对缺失值比重复测量方差分析更为耐受。它是基于完全随机缺失的假设(关于随机缺失等概念参见以前文章，下一篇文章也会再次介绍)，因此完全随机缺失模式对广义估计方程的结果影响不大，此时其参数估计值仍是稳定的，但如果是随机缺失，仍会影响广义估计方程的结果，这种情况下，可 考 虑 加 权 的 广 义 估 计 方 程(Weighted Generalized Estimating Equations)，该法是基于随机缺失的假定，但仅限于失访模式(即一个人在某个时间点缺失后，后面的时间点均无数据)。
