某溶栓药物治疗20名急性脑梗死患者的疗效，采用随机、双盲、安慰剂平行对照设计，每组各10例，分别于治疗前及治疗后8周每周进行随访观测，观测指标为神经系统体征评分（MDNS）。
示例来源：杨珉.李晓松等.医学和公共卫生研究常用多水平统计模型.北京:北京大学医学出版社,2007.5.
该示例我们曾在重复测量数据分析系列笔记的《基于Stata再谈多层混合效应模型》及《广义线性混合模型(GLMM)》中进行过分析。本次笔记对该示例采用广义估计方程进行分析，使用软件stata16。
转自个人微信公众号【Memo_Cleon】的统计学习笔记：重复测量数据分析系列_广义估计方程（Stata）。
广义估计方程在stata中的命令是xtgee，菜单操作比较简单，步骤如下：
统计>>广义线性模型>>广义估计方程(GEE)
或
统计>>纵向/面板数据>>广义估计方程(GEE)>>广义估计方程(GEE)
因变量：MDNS；
自变量：age、Trtg、time。其中Trtg实际是分类变量，二分类变量跟按连续变量分析的结果是一致的，因此按连续变量纳入；进入面板设置，在面板ID变量下拉框中选入变量id；分布和关联选择：视数据类型选择合适分布类型和关联函数，本例默认高斯分布和恒等。
【相关系数】选择组内相关结构，即确定模型的作业相关矩阵。选择哪种作业相关矩阵，可以由相关系数矩阵来大体判断，或者通过通过信息准则比如QIC来选择的更优模型来确定。比较遗憾的是在stata的gee后验估计中并未提供模型的信息准则，组内相关结构可选择[无结构]，然后通过后验估计工具来获取相关矩阵来大体判断。好在在广义估计方程中，无论作业相关矩阵指定正确与否，对参数的影响并不大，尤其是采用稳健估计时能够保证检验结果的一致性。本例默认组内的相关结构是可交换。
【SE/Robust】默认常规标准误，本例选择[稳健]，缩放因子选择[使用除数N-P而不是N]
接着是模型一些基本情况。模型迭代1次即收敛，采用高斯分布恒等连接函数，作业相关矩阵为可交换，患者20人，共有180个观测。模型纳入了3个变量，相比不纳入变量具有统计学意义（Wald chi2(3)=97.15,P<0.001）。Trtg、time的主效应均有统计学意义。校正其他因素的影响后，治疗组的MDNS比对照组高11.86分（安慰剂组赋值为0，治疗组赋值1，本例按连续变量纳入），每过1周，MDNS升高2.81分。年龄对MDNS无影响。
我们建立的上述模型并不含有交互项，其前提是两个组的MDNS随时间变化而变化，图示上看两条线（两个组不同时间的MDNS）是平行的。我们可以构建分组的带有拟合线的散点图来判断一下（图形>>二维图）。从图上看，两组的增长趋势不同，治疗组更快一些。当然统计学上的相等并一定是绝对相等，我们接下来纳入交互项来验证一下。
在GEE操作界面将自变量部分纳入交互项，操作如下，注意我们将第一个水平设为参照水平，其他同前：
模型中加入了4个变量，模型具有统计学意义（Wald chi2(4)=707.2,P<0.001）；模型在加入交互项后，各因素的系数已不再表示其主效应，而是简单效应，即一个变量在另一个变量的一个级别上的效应。Trtg的系数-0.32表示校正年龄因素，在治疗前(time=0)，治疗组的MDNS比对照组第0.32分，没有统计学意义，即治疗前两个组没有统计学差异。time的系数则表示校正年龄因素的影响，对照组每过1周，MDNS增加1.28分，具有统计学意义。交互项可以有两种理解，一种是治疗组和对照组之间的MDNS差值在不同的时间点是不同的，每过1周这个差值会增加3.05分；当然也可以理解为治疗组每过1周MDNS的增加值比对照组每过1周MDNS的增加值多3.05分。
本例将组内重复因素time按连续变量纳入，相当于默认了随着时间的延长，MDNS的变化率是相同的，即平均变化率，但事实往往并非如此，如果我们想具体了解每个时间点真实的增加情况，可以把时间因素作为分类变量纳入，时间因素就会按虚拟变量输出结果。另外由于连续变量的“参照水平”是低水平，所以Trtg的系数代表的是治疗前治疗组与对照组的差值，结果往往无统计学意义，在将时间因素time设为分类变量时，不妨将参照水平设为高水平，这样以来Trtg的系数代表的是治疗8周后治疗组与对照组的差值。自变量纳入如下：
结果如下，解读可参照上述结果就不再赘述了。
