同一批对象反复测了四五次，指标一路往下掉，你要是图省事直接丢进 glm，软件会老老实实给你一堆漂亮到不太真实的 p 值——它把两千条记录当成两千个互不相干的独立个体，把个体内的相关性直接当噪音忽略了。这个坑我在好几个项目里都见过，最离谱的一次，同一个体内部高度相关的随访数据按独立处理，一个边缘显著的效应被吹成了极其显著，结论差点写进报告。广义估计方程（GEE）就是专门用来对付这类重复测量、面板、整群抽样数据的：它不去建模个体层面的随机效应，而是把“组内相关”写成一个工作假设塞进方差结构里，再用稳健的沙盒估计把标准误拉回合理区间。这篇围绕 R 实战讲 GEE，从它到底解决什么问题、geepack 等几个包怎么挑、数据怎么整形、模型怎么拟合、相关结构怎么比、结果怎么导出、报错怎么排查，一条线走完。适合手上有纵向随访或聚类数据、正在纠结该用混合模型还是 GEE、或者被 glm 的假显著坑过一次想弄明白原因的人。代码都是能直接跑通的，数据我全部用模拟的，你复制走就能改自己的变量名。

1. 先搞清楚 GEE 到底在解决什么问题

1.1 重复测量数据为什么会让普通回归翻车

拿一份最典型的随访数据举例：300 个人，每人测 5 次，一共 1500 条记录。同一人的 5 次测量之间天然是相关的——体质偏胖的人五次都偏胖，血压高的人五次都偏高。这种相关性在统计上叫组内相关（intra-cluster correlation，ICC）。普通 glm 的似然函数假设所有观测独立，它算出来的信息量是 1500 份独立信息，但真实有效信息远小于这个数，可能只相当于 300 份多一点。

结果就是标准误被系统性低估，置信区间过窄，p 值偏小。这不是“略微不准”，在 ICC 达到 0.5、每组 5 次的场景下，标准误能被低估一半以上，等于凭空多出来一倍的样本量。我做过的项目里，凡是有人拿着几万条“记录数”在汇报里写“样本量 N=20000”，我第一反应就是问他实际受试者多少人，十有八九只有几千。

更麻烦的是，这种错误不会在软件输出里报警。R 不会提示你“这里有相关性没处理”，它只会安安静静跑完，给你一个看起来人模人样的 summary。所以判断责任完全在分析者这边：只要你的数据里存在“同一个人/同一个班级/同一家门店/同一个批次”这种分组标识，就必须先问一句——这些观测真的独立吗？

1.2 边际模型与条件模型：GEE 和混合模型的分岔口

处理组内相关，主流有两条路：一条是广义线性混合模型（GLMM，R 里常用 lme4 的 glmer、GLMMadaptive 等），另一条就是广义估计方程。这两条路不是谁替代谁，回答的是两个不同的科学问题，选错了不是技术错误，是问题错误。

混合模型是条件模型（subject-specific）。它估计的是“对同一个个体而言，处理组相对对照组的效应是多少”。系数解释里带着“对同一个体”“给定随机效应”这样的限定语。它需要你假定随机效应的分布形式，通常是正态，这个假定没法完全验证，而且在非线性链接（比如 logit）下，条件效应和边际效应的数值差别不小。

GEE 是边际模型（population-averaged）。它估计的是“把整个人群平均来看，处理组的结局指标相对对照组高多少”。系数直接对应人群层面的对比，适合做公共卫生、政策评估、群体层面的结论。它的代价是不建模个体差异的分布，只把相关当成一个待估计的冗余参数。

所以选择逻辑很清楚：如果你的结论要落到“某一类人群整体上会怎样”，选 GEE；如果要落到“针对某一个具体个体，干预会带来什么变化”，选混合模型。我见过太多论文在这个点上含糊其辞，用 GEE 的系数却按条件效应的口径解释，审稿人一问就露馅。

1.3 什么场景该选 GEE，什么场景该转身走

把适用场景摊开说会更清楚。适合 GEE 的典型情况有这么几类：

• 纵向随访数据，关心的核心是处理组与时间之间的交互，人群平均效应是报告重点；
• 整群随机设计，比如以班级、社区、门店为单位分组的干预试验，分析单位是个体但分组是整群；
• 横断面聚类数据，比如同一家医院的多位患者，只关心医院层面的平均差异；
• 结局是非正态的计数、二分类、率，且组内相关的具体机制说不清楚、不想为它设定分布假定。

反过来，下面几种情况我会劝你先别急着用 GEE：

• 组数很少。稳健标准误的大样本性质依赖组数（cluster 数）足够多，三五十组以下就要非常谨慎，个位数组的场景基本别用；
• 存在大量缺失，而且缺失机制明显不是随机缺失（MAR）方向，GEE 的稳健性依赖均值模型设定正确，缺失一旦带信息，偏差会直接进到点估计里；
• 你想做个体层面的预测或随机效应排序。GEE 根本不提供 BLUP，做不了这件事，该换混合模型；
• 数据不平衡得很厉害，每个人的随访次数从 1 到 20 不等且差异有系统性，虽然 GEE 理论上能处理不平衡设计，但相关结构的估计会变得很飘。

注意：GEE 处理缺失的前提是“完全随机缺失（MCAR）”在严格意义下需要的，广义上在均值模型设定正确时对 MAR 有一定稳健性，但这个稳健性经常被过度解读。缺失比例超过 20% 且缺失与结局相关时，别指望 GEE 帮你兜底。

2. GEE 的底层逻辑：拟似然与稳健沙盒

2.1 从得分方程说起，用人话讲一遍

GEE 的名字里有“方程”两个字，说明它不是靠最大化一个似然函数来估参数的，而是靠解一组方程。这组方程的形式是：每个个体的观测值减去模型预测值的残差，乘上一个权重矩阵，再把所有人加起来等于零。解出这个方程组，就得到系数估计。

关键在于那个权重矩阵。它由两部分组成：一部分是方差函数（由 family 决定，比如泊松的方差等于均值），另一部分是组内相关结构（由 corstr 决定）。这两部分合起来构成“工作协方差矩阵”。注意“工作”这两个字——它意味着这只是个假设，不是你声称的真实结构。

这就是 GEE 最妙的地方：即使你把相关结构设错了（比如真实是 AR(1)，你设成了 exchangeable），点估计在均值模型正确的前提下依然是一致的，只是效率会下降——标准误变大，但不会系统性偏移。这个性质叫“对工作相关结构错误设定稳健”。我第一次看到这个结论时是存疑的，后来自己跑了模拟：分别用 independence 和 exchangeable 拟合同一份强相关数据，点估计几乎一模一样，标准误差了将近 40%。

2.2 工作相关结构怎么选，选错了会怎样

corstr 这个参数是 GEE 实战里最需要动脑子的地方。R 里 geepack 支持的常见选项和它们的适用条件，我用一张表归纳一下：

实际项目里，我默认从 exchangeable 起步，因为它只需要估一个相关参数，对不平衡设计也能算，容错性最好。如果测量次数固定、间隔相等，而且怀疑相关性随时间间隔衰减，再试 ar1。unstructured 我一般只在组内次数不超过 4、样本量充足的时候碰，它太吃数据了。

有个容易忽略的点：exchangeable 与 independence 在平衡数据下、使用典范链接时，点估计完全相同。这是因为工作相关结构在两类模型的拟似然权重里被抵消了。所以如果你发现换了相关结构系数一点没变，别怀疑代码写错了，这是正常的。

2.3 稳健标准误为什么在组数足够时才靠谱

GEE 的标准误有两种：朴素（naive）和稳健（robust，也叫 sandwich）。朴素那个假设工作相关结构完全正确，一旦设错就是错的；稳健那个用残差的外积去估计方差，不需要相关结构正确也能一致。

教科书上都说稳健标准误“不需要正确设定相关结构”，听起来像免费午餐，但有个前提经常被一笔带过：它是大样本性质，而这里的“样本”指的是组的个数，不是观测条数。我做过一组模拟，同样的总观测数 5000，分成 50 组和分成 500 组，稳健标准误的表现差别很大。组数在 40 以上时，稳健标准误的覆盖率基本能靠近名义水平；组数掉到 20 以下，就会出现明显的向下偏倚，假阳性率上升；组数个位数的时候，稳健标准误可能比朴素标准误还不可靠。

所以做整群抽样分析前，第一件事是数一下 cluster 有多少个。组数少的时候，稳妥做法是用 Jackknife 或 cluster bootstrap 修正，或者干脆退回去用混合模型加 Kenward-Roger 类的小样本自由度校正。这条经验我是踩过坑换来的：某次分析只有 12 个社区，用 GEE 跑出来的 p 值 0.03，换成 cluster bootstrap 之后置信区间直接跨过 0，结论完全反转。

提示：geeglm 的 summary 默认输出的 Std.err 就是稳健版本。想对比朴素版本，用summary(fit, std.err = "naive")，把两个结果并排放，差距大小能直接反映相关结构设得准不准。

3. R 里的工具箱：geepack、gee、geeM 怎么挑

3.1 三个主流包的能力边界

R 里做 GEE 的包不算多，常用的就三个，各自定位不一样。

geepack是我用得最多的。核心函数geeglm()，接口设计得跟glm()高度一致，family、formula、data 这些参数照搬，学习成本很低。它支持 independence、exchangeable、ar1、unstructured、m-dependent 这几种结构，配套提供QIC()做模型比较。缺点是它对数据顺序比较敏感，跑之前最好显式排序。另外它在估计相关系数时用的是矩估计，当组内次数很小时相关系数估计会偏。

gee是老前辈，函数就叫gee()。接口比 geepack 更原始一些，参数名不太一样，返回对象的结构也不同。它的优点是实现稳定、年代久，很多老代码和老教材都基于它。缺点是维护频率低，输出不够友好，做模型比较得自己算 QIC。

geeM相对新，函数geem()。它的卖点是支持更灵活的相关结构定义，还能处理一些 geepack 处理起来吃力的场景，比如分类型的相关结构。如果你的相关结构比较特殊，或者需要跟别的软件（比如某商业统计软件）的结果对齐做验证，geeM 值得一试。

我的默认选择是 geepack，只有两种情况会换：需要复现别人的老代码时用 gee；相关结构太特殊、geepack 里的选项覆盖不了时用 geeM。

3.2 安装与版本核对

安装没什么难度，但版本要注意，geepack 在不同版本间对 QIC 的返回结构有过调整。装上之后先跑一下版本号和函数签名确认：

install.packages("geepack") library(geepack) packageVersion("geepack") # 顺手看一眼 QIC 的返回结构，不同版本字段名可能不同 args(geepack::QIC)

如果你的项目要求可重复性，建议把包版本写进分析脚本的注释里，或者在项目里用 renv 锁定。我吃过一次亏：半年后重跑同一份代码，QIC 的数值对不上，查了半天发现是包升级后 QICu 的计算有小改动，白折腾一晚上。

另外提醒一句，geeglm 的 id 参数可以传数值向量、因子或者字符。传数值最省事，传因子它会按水平处理，但打印出来的报错信息里 id 名称会变成因子水平，排查问题时容易看晕。我一般统一把 id 列处理成整数，并且保证它从 1 开始连续编号，这样后面做 bootstrap 分块抽样也方便。

4. 数据准备：长格式、排序与 id 的那点事

4.1 宽转长的实操

GEE 只吃长格式（long format）数据。所谓长格式，就是每一行是一次观测，同一个人在数据里出现多次，用一个 id 列标识。宽格式是一行一个人，各次测量横着摊开（y1、y2、y3 这样的列名）。

宽转长用reshape()或者 tidyr 的pivot_longer()都行。我倾向用 pivot_longer，语法清楚，出错信息也友好：

library(tidyr) # 宽格式：每个受试者一行，baseline / m3 / m6 / m12 是四个时点的测量 # 转成：每个受试者四行，多出 wave 和 y 两列 long_df <- pivot_longer( wide_df, cols = c(baseline, m3, m6, m12), names_to = "wave_label", values_to = "y" ) # 把时点标签转成数值，方便做时间趋势项和 AR(1) long_df$wave <- as.integer(factor( long_df$wave_label, levels = c("baseline", "m3", "m6", "m12") ))

这里有个细节值得多说一句：factor(..., levels = ...)这一步不能省。pivot_longer 转出来的 wave_label 默认按出现顺序整理，如果你原始数据列的顺序不是时间顺序，转出来的 wave 编号就会错位。AR(1) 结构对时间顺序极其敏感，错一位结果就全歪了。显式指定 levels 是最省心的做法。

4.2 排序与缺失：最容易埋雷的两步

转换完之后，两件事必须做，而且要在拟合之前做。

第一件是按 id 和 wave 排序。geeglm 内部会按 id 把观测分组，但它对输入顺序的依赖让人不太放心，尤其当你的 id 在数据里不连续出现时（比如 id=1 的记录出现在第 1 行和第 800 行）。显式排序一行搞定：

long_df <- long_df[order(long_df$id, long_df$wave), ] rownames(long_df) <- NULL

第二件是检查缺失和不平衡。GEE 能处理不平衡设计，但你要清楚自己面对的是什么形态。用一条简单的汇总就能看清楚每人测了几次：

n_per_id <- table(long_df$id) summary(as.integer(n_per_id)) # 再看看有多少人只测了一次——只测一次的人对组内相关估计没有贡献 sum(n_per_id == 1)

只测一次的对象在 exchangeable 结构下对相关系数估计没用，但他们会进入均值模型，所以不能随便删——删掉会造成选择偏倚。正确做法是保留，同时在做敏感性分析时单独跑一份“只保留至少两次测量”的版本，看结论是否改变。

手工手写一句：缺失值的处理必须在建模前想清楚。GEE 默认的 na.action 是 na.omit，也就是整行删。如果你的缺失比例不高（比如 3% 以内）且分散，整行删问题不大；如果缺失集中在某个随访时点，整行删会直接把那个时点的样本结构搞变形，最好先看看缺失模式。

4.3 协变量与时间的编码

时间变量的编码方式直接决定了交互项怎么解释。三种常见编码各有用途：

• 把 wave 当连续变量（1,2,3,4），交互项系数解释为“每往后一个时点，组间差异的变化量”，前提是你认为时间效应是线性的；
• 把 wave 当因子，交互项变成各时点分别的组间差异，信息量大但参数多，样本不够时标准误很宽；
• 用实际时间（0, 3, 6, 12 个月），适合随访间隔不等且你想让系数对应“每月变化”。

我一般先用实际时间当连续变量跑一版，看看交互项的趋势是否单调，如果明显不平（比如先升后降），再换成因子版本细看。直接上因子版本的问题是，如果每个时点都有交互，模型参数很快膨胀，尤其样本量不宽裕时，QIC 的选择也会变得不稳。

另外，把基线值作为协变量放进模型是个常见做法，能提升效率（这叫基线校正）。但要注意基线值在每个人内部是不变的，它不会引入组内相关，放进去是安全的，只是别和随机效应搞混。

5. 从零跑通一个完整 GEE 模型

5.1 模拟一份纵向计数数据

用一份完全模拟的数据把流程走一遍。设定 300 个受试者，每人随访 5 次，结局是计数型（可以理解为某类事件的次数），有一个二分类干预变量和一个时间变量，两者有交互。为了让相关结构真实存在，我给每个受试者加了一个随机截距，用泊松分布生成计数：

set.seed(2024) n_id <- 300 n_wave <- 5 dat <- expand.grid(id = 1:n_id, wave = 1:n_wave) # 给每个人分配一个干预组别 trt_map <- data.frame( id = 1:n_id, trt = sample(c(0, 1), n_id, replace = TRUE) ) dat <- merge(dat, trt_map, by = "id", sort = FALSE) dat <- dat[order(dat$id, dat$wave), ] rownames(dat) <- NULL # 个体层面的随机截距，制造组内相关 b0 <- rnorm(n_id, mean = 0, sd = 0.8) dat$time <- dat$wave - 1 linpred <- -0.5 + 0.6 * dat$trt + -0.20 * dat$time + -0.15 * dat$trt * dat$time + b0[dat$id] dat$y <- rpois(nrow(dat), lambda = exp(linpred)) # 看一眼形态 head(dat, 10) table(dat$trt)

这段代码里的b0[dat$id]是关键：它按 id 索引把个体随机截距对应到每一行。因为 dat 已经按 id 排过序，这个索引是准确的。如果你不排序就直接索引，随机截距会错位到别的人身上，后面所有结论都是假的——这个错误极难发现，因为代码不报错，结果看着还挺像那么回事。

5.2 拟合模型与参数解读

先跑 independence，作为对照，然后跑 exchangeable：

library(geepack) fit_ind <- geeglm( y ~ trt * time, data = dat, id = id, family = poisson(link = "log"), corstr = "independence" ) fit_ex <- geeglm( y ~ trt * time, data = dat, id = id, family = poisson(link = "log"), corstr = "exchangeable" ) summary(fit_ex)

输出里除了系数表，下面还有一段 Correlation 部分，给出估计出来的组内相关参数 alpha。这个 alpha 很重要，它告诉你同一人两次测量之间的相关性有多强。如果 alpha 接近 0，说明 independence 就够用了；如果 alpha 在 0.3 以上，说明忽略它确实会低估标准误。

系数解读要小心，因为用了 log 链接：

• trt的系数是 baseline 时点（time=0）处理组相对对照组的 log 率比；
• time的系数是对照组每多一个时点的 log 变化；
• trt:time的交互项是“处理组每时点变化量减对照组每时点变化量”，也就是两条趋势线的斜率差。

转成率比要取指数，并且给置信区间也做同样的变换：

beta <- coef(fit_ex) se <- summary(fit_ex)$coefficients[, "Std.err"] rate_ratio <- exp(beta) ci_low <- exp(beta - 1.96 * se) ci_high <- exp(beta + 1.96 * se) data.frame( term = names(beta), RR = round(rate_ratio, 3), CI_low = round(ci_low, 3), CI_high = round(ci_high, 3) )

这里用的是 1.96 的正态近似。组数在 40 以上时这么算没问题；组数偏少时应该换成 t 分布分位数，自由度取组数减参数个数，虽然也只是个近似，但比 1.96 靠谱。

5.3 换相关结构重跑，对比结果

跑三四种相关结构，把点估计和标准误并排放，这是我最推荐的一步自查：

fit_ar1 <- geeglm(y ~ trt * time, data = dat, id = id, family = poisson(link = "log"), corstr = "ar1") fit_un <- geeglm(y ~ trt * time, data = dat, id = id, family = poisson(link = "log"), corstr = "unstructured") compare <- data.frame( term = names(coef(fit_ex)), est_ind = round(coef(fit_ind), 4), est_ex = round(coef(fit_ex), 4), est_ar1 = round(coef(fit_ar1), 4), est_un = round(coef(fit_un), 4), se_ind = round(summary(fit_ind)$coefficients[, "Std.err"], 4), se_ex = round(summary(fit_ex)$coefficients[, "Std.err"], 4), se_ar1 = round(summary(fit_ar1)$coefficients[, "Std.err"], 4), se_un = round(summary(fit_un)$coefficients[, "Std.err"], 4) ) compare

看结果的方式很明确：点估计那一排如果有明显差异（比如相对差异超过 20%），说明均值模型可能有问题，而不是相关结构的问题，得回去检查模型设定。标准误那一排如果有明显差异，说明工作相关结构的选择确实影响效率，但点估计不变就说明结论方向是稳的。

实测下来，exchangeable 和 ar1 在测量间隔均等的模拟数据上，标准误往往很接近；independence 的标准误会明显偏小；unstructured 因为参数最多，标准误可能略小但估计过程更不稳定。如果 unstructured 报错或者给出警告，别硬扛，退回 exchangeable 就行。

5.4 用 QIC 做模型比较

选相关结构不能光看标准误大小——标准误小的模型不一定是“对”的。规范做法是用 QIC（quasi-likelihood under the independence model criterion），思路跟 AIC 类似，越小越好。

qic_tab <- data.frame( structure = c("independence", "exchangeable", "ar1", "unstructured"), QIC = c(QIC(fit_ind)$QIC, QIC(fit_ex)$QIC, QIC(fit_ar1)$QIC, QIC(fit_un)$QIC), QICu = c(QIC(fit_ind)$QICu, QIC(fit_ex)$QICu, QIC(fit_ar1)$QICu, QIC(fit_un)$QICu) ) qic_tab[order(qic_tab$QIC), ]

用 QIC 有个重要前提：比较均值模型（也就是选协变量）时用 QICu，比较工作相关结构时用 QIC，两类比较不要混在一起做。这个区分在很多教程里被一笔带过，但搞混了会得出错误结论。我见过有人用 QIC 去比两个协变量集合不同的模型，然后堂而皇之地说“加了交互项的那个更好”，这是不成立的。

注意：QIC 的差值不像 AIC 那样有明确的经验阈值。我个人的用法是：差值在 2 以内视为没有实质差别，这时候优先选参数更少、更简单的结构；差值超过 10 才认为有明显偏好。

6. 结果呈现与稳健性检查

6.1 系数表怎么导出成能直接进报告的样子

geepack 的 summary 输出是给人看的，不是给表格工具看的。做报告我一般写个小函数，把系数、稳健标准误、z 值、p 值和置信区间一次整理成数据框：

tidy_gee <- function(fit, conf_level = 0.95) { s <- summary(fit) est <- s$coefficients[, "Estimate"] se <- s$coefficients[, "Std.err"] z <- est / se p <- 2 * pnorm(-abs(z)) q <- qnorm(1 - (1 - conf_level) / 2) data.frame( term = rownames(s$coefficients), estimate = round(est, 4), robust_se = round(se, 4), z = round(z, 3), p_value = signif(p, 3), ci_low = round(est - q * se, 4), ci_high = round(est + q * se, 4), row.names = NULL ) } tidy_gee(fit_ex)

非线性链接下，我会同时给两列：一列是链接尺度上的系数（log 尺度），另一列是取指数之后的效应量（率比或优势比）。审稿人或读者看得懂的是后者，但做元分析合并的时候需要前者。两列都给，谁用都不别扭。

6.2 边际效应与预测的正确姿势

GEE 输出的是边际模型的系数，如果要给出“处理组人群的平均率是多少”这类具体数字，需要做预测。geeglm 对象支持 predict，但要注意参数：

# 构造两个人群：处理组和对照组，在每一个时点上 newdat <- expand.grid( trt = c(0, 1), time = 0:4 ) pred_link <- predict(fit_ex, newdata = newdat, type = "link") pred_resp <- predict(fit_ex, newdata = newdat, type = "response") newdat$fit_link <- pred_link newdat$fit_resp <- pred_resp newdat

type = "response"给的是人群平均尺度上的预测率。但这里有个坑：由于链接函数是非线性的，对 log 尺度的线性预测子取指数得到的，是条件于协变量的均值，而不是把随机性和相关结构都积分掉之后的人群均值。两者在小方差时不差多少，方差大时会有可见差异。要更严谨的人群平均预测，用 cluster bootstrap：按 id 有放回地整块抽，重复几百次，每次拟合模型再预测，最后取平均和分位数。这套流程代码略长，但优势是顺带把标准误的问题也解决了，组数少的时候尤其值得做。

6.3 敏感性分析该做哪几项

一个 GEE 模型跑完，我至少会做三项敏感性检查，这些检查能让结论的可信度提升一个档次：

第一项是换相关结构。前面已经跑过，看结论方向是否变化，如果 exchangeable 和 ar1 得出的交互项符号相反，那就别急着下结论，先回去查数据。

第二项是换缺失处理方式。把“整行删”换成“只用至少两次测量的对象”，或者做多重插补后再合并。多重插补做 GEE 有个注意点：插补模型里必须包含 cluster 信息，否则插补本身会破坏组内相关结构。R 里 mice 支持设置 cluster 变量，是一定要设的。

第三项是做 cluster bootstrap。按组重抽样，重复 500 到 1000 次，记录每次的系数，最后用分位数给置信区间：

set.seed(123) ids <- unique(dat$id) B <- 500 boot_b <- matrix(NA, nrow = B, ncol = length(coef(fit_ex))) colnames(boot_b) <- names(coef(fit_ex)) for (b in 1:B) { samp_id <- sample(ids, length(ids), replace = TRUE) boot_dat <- do.call(rbind, lapply(seq_along(samp_id), function(i) { tmp <- dat[dat$id == samp_id[i], ] tmp$id <- i tmp })) boot_dat <- boot_dat[order(boot_dat$id, boot_dat$wave), ] fb <- try(geeglm(y ~ trt * time, data = boot_dat, id = id, family = poisson(link = "log"), corstr = "exchangeable"), silent = TRUE) if (!inherits(fb, "try-error")) boot_b[b, ] <- coef(fb) } apply(boot_b, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)

这段代码里把重抽样后的 id 重新编号为 1 到 n 是有意为之。如果不重编号，同一个原始 id 抽中两次，会有两组观测共用一个 id，geeglm 会把它们当成同一个人，方差估计直接错掉。这个细节我在第一次写 bootstrap 时完全没想到，结果标准误小得离谱，排查了很久才发现是这里。

7. 常见报错与排查实录

7.1 不收敛与“协方差矩阵非正定”

最常遇到的报错是 unstructured 结构下的非正定问题，提示大意是相关矩阵不可逆或者估计失败。原因是 unstructured 要估 T(T-1)/2 个相关参数，T=5 就是 10 个参数，加上均值模型的参数，很容易在组内次数少或缺失多的场景下解不出来。

处理顺序我一般是这样的：先看每个组的观测次数分布，如果有很多人只测了一两次，unstructured 基本没戏，直接放弃；再看是不是时间点顺序错了，AR(1) 和 unstructured 都对顺序敏感，顺序错了矩阵就乱了；最后才考虑加约束，比如改用 exchangeable 或者 ar1。

还有一种不收敛表现为迭代次数用尽。geeglm 有个maxit参数，默认值不大，复杂模型可以调到 50 甚至 100。但调大之前先想想，正常模型不该需要那么多次迭代，不收敛往往说明模型设定本身有问题，比如某个协变量在每个组内都不变（这种情况下它和随机效应完全共线），或者某个分类变量存在极端稀疏的水平。

7.2 标准误异常与 id 设置错误

如果标准误比朴素版本还小得离谱，或者置信区间窄到不真实，第一件事查 id。常见错误有三种：

一是 id 传错，把行号或者 wave 当成了 id。这样每个观测都被当成独立的组，等于没做 GEE，还多了一层误以为处理过的错觉。二是数据没排序，id 相同的记录散落在各个位置。三是 id 是字符型且有空值，geeglm 内部转换时会把空值那几行悄悄丢掉。

排查方法很简单，拟合完看一眼组数：

length(unique(dat$id)) fit_ex$geese$nobs # 参与拟合的观测数

如果你的 id 数量是 300，但拟合里报告的组数也是 1500，那肯定是 id 传错了。另外，fit_ex$geese$vbeta是稳健协方差，fit_ex$geese$vbeta.naiv是朴素协方差，两个矩阵对角线元素开根号对比一下，差距如果在 1.2 倍以内，说明相关结构设得还不错；差距超过 3 倍，说明组内相关很强而你的结构没抓住。

7.3 速查表：几个高频问题的定位思路

7.4 一个容易被忽略的坑：分类变量的参照水平

GEE 的系数解释依赖参照水平，而参照水平在 R 里默认为字典序第一个。如果你的处理变量编码成 "对照" 和 "干预"，默认参照是 "对照"（拼音顺序），看着没问题；但如果编码成 "A组" 和 "B组"，参照就是 "A组"，万一你心里想的是 B 组当参照，所有系数符号都会反，而输出里不会提示你。

我的习惯是在建模前显式重设因子水平，用relevel()或者factor(x, levels = c("对照", "干预"))，把参照组放在第一位。同时把这一步写进脚本，不要靠交互式操作，否则下次重跑就忘了。

8. 几条从实际项目里磨出来的经验

8.1 组数少的时候，别硬扛

这一条我想单独拎出来说，因为它是最容易翻车的点。组数少于 30 时，GEE 的稳健标准误就开始打折扣；少于 15 时，我基本不建议把 GEE 的 p 值当结论依据。

替代方案有三个，按推荐顺序排：cluster bootstrap（最通用，实现稍麻烦但结论最可信）、Jackknife 修正（计算量比 bootstrap 小，适合组数在 15 到 30 之间）、换成混合模型配合小样本自由度校正。第三种方案要注意，混合模型在非线性链接下估的是条件效应，和 GEE 的边际效应不可直接比较，换模型的同时报告口径也要跟着改。

我个人的判断标准是这样的：组数 40 以上，GEE 加稳健标准误直接用；30 到 40，用 t 分布分位数而不是 1.96，同时补一份 bootstrap 对照；30 以下，bootstrap 作为主要结果，GEE 的解析标准误只作参考；10 以下，直接放弃 GEE。

8.2 时间变量到底该不该放进相关结构

有个经常被问到的问题：既然 exchangeable 把组内任意两次测量的相关性视为相同，那我是不是还得把 time 放进均值模型里？答案是必须放，这两件事完全独立。相关结构管的是残差之间的相关性，均值模型管的是趋势。把 time 从均值模型里拿掉，残差里就会残留时间趋势，这个趋势会被误当成组内相关吸收到 alpha 里，导致 alpha 被高估，标准误跟着变形。

反过来的错误也有：有人为了让模型“更简洁”，把 time 从均值模型里拿掉，只保留处理组，然后用 exchangeable 硬扛。这样做的后果是交互效应完全没法估计，而纵向研究里交互往往才是核心问题。

正确做法是把时间趋势建模充分——线性不够就加二次项，或者分段，或者直接用因子。判断标准很简单：残差里不该还有系统性的时间模式。看一眼按 time 分组的平均残差就能发现。

8.3 报告里必须写清楚的四件事

最后说报告规范。一份能被复现的 GEE 分析，方法部分至少要交代四件事，缺一样都算不合格：

写清楚工作相关结构选了哪个，以及为什么选。别只写“使用了 GEE”，审稿人一定会问相关结构。如果做过 QIC 比较，把比较结果给一个数，哪怕只是“exchangeable 的 QIC 最小”。写清楚标准误是稳健还是朴素，这一条最容易被省略，但恰恰是最影响结论的。写清楚缺失数据的处理方式和缺失比例，以及是否做过敏感性分析。写清楚组数（不是观测条数），以及是否因为组数偏少做了额外处理。

至于结果呈现，我的偏好是主表给链接尺度系数加稳健标准误和置信区间，附表给取指数后的效应量。交互项一定要单独解释，别笼统说“干预有效”，要说清楚效应随时间怎么变化。如果交互项不显著但趋势明显，可以给一张预测曲线图，比一堆系数直观得多。

我个人在这些项目里最深的一点体会是：GEE 的技术门槛其实不高，geeglm 那一行代码谁都会写，真正区分分析质量的是前面那 80% 的准备工作——数清楚组数、排好序、想明白相关结构、把缺失摊开看一遍。这些事没有一行代码能替你完成，但它们决定了最终那个 p 值到底是真信号还是算法给你的安慰。