R语言广义估计方程GEE实战：纵向数据相关结构与稳健标准误

1. 纵向数据的相关性陷阱，以及广义估计方程到底在算什么

做过随访数据的人大概都有过这种经历：同一个人测了四次，直接把这四行当成四个独立样本扔进 logistic 回归，结果出来一堆"显著"变量，效应量看着还挺漂亮，等到审稿意见回来一句"你是否考虑了同一个体内观测的相关性"，前面的功夫基本白费。**广义估计方程（Generalized Estimating Equations，GEE）**就是专门用来收拾这类数据的工具，而 R 里的geepack包则是把它落地成代码最省心的路径之一。我写这篇东西的目的很明确：给已经会跑普通回归、手上又有重复测量或聚类数据的人，一套能照着抄、抄完能解释清楚结果的完整流程。

先把这个模型的位置说清楚。GEE 属于边际模型（marginal model），它估的是"如果整个群体都接受处理 A 对比处理 B，平均会差多少"这种总体层面的效应，而不是"某一个具体的人从处理 A 换到处理 B 会差多少"那种个体条件效应。两者的区别在连续结局、线性连接下几乎看不出来，但只要换成二分类结局配 logit 连接，条件 OR 和边际 OR 就会明显分家，这一点后面会用模拟数据实测给你看。适合读这篇的人包括：手里有面板调查、临床试验随访、重复测量实验、用户留存或复购序列数据的分析人员；已经会用 GLM 但还没系统处理过聚类结构的同学；以及那些被"要不要用混合模型"这个问题卡住过的人。整篇会覆盖数据准备、相关结构选型、模型跑通、QIC 比较、结果表达、报错排查六个环节，代码全部基于模拟数据，可以直接复制运行。

1.1 同一体内观测为什么不能当独立样本

假设有个随访研究，每个受试者在第 0、1、2、3 个月各测一次结局。如果不加任何处理直接建模，模型会默认这 800 行观测彼此毫无关系。但实际情况是，第 0 个月和第 1 个月的测量来自同一个人，他的体质、基线水平、依从性这些没被观测到的因素会同时影响四次测量，导致同一体内的残差是正相关的。这种相关性带来的直接后果是标准误被低估：你以为样本量是 800，模型实际能拿到的独立信息量可能只有 200 左右，置信区间会被压得过窄，p 值被压得过小。

用生活化的类比讲，这就像调查一个班级的学习情况，你问了 40 个学生各 5 个问题，得到了 200 条回答。如果这 200 条回答被当成 200 个独立个体，你对"这个班级的平均水平"的把握就被严重高估了——因为同一个学生的 5 个回答高度相似，真正提供新信息的是学生之间的差异，而不是同一学生内部重复问。GEE 的处理办法很聪明：它不去建模"个体内部的随机效应长什么样"，而是直接对**工作相关结构（working correlation structure）做一个假设，然后用稳健标准误（sandwich estimator）**把假设错了的后果兜住。即使你把相关结构猜错，只要均值模型设定正确、聚类数足够，点估计依然是一致的，标准误也依然近似正确。

1.2 GEE 与混合效应模型的选择逻辑

这是被问得最多的问题，我给一个可以直接照着判的决策口径。要回答"群体平均效应是多少"，选 GEE；要回答"个体层面的变异有多大、某个人的预测轨迹是什么"，选 GLMM（广义线性混合模型）。前者输出的是**总体平均（population-averaged）的系数，后者输出的是条件（subject-specific）**的系数。在有非线性连接函数的情况下，条件 OR 的绝对值通常大于边际 OR，因为边际化过程中把个体异质性平均掉了，产生所谓的"衰减"。

实操层面的取舍还要看三件事。第一是缺失机制：GEE 在数据完全随机缺失（MCAR）下表现良好，在随机缺失（MAR）下需要加权 GEE 或直接改用似然方法（GLMM 在 MAR 下无偏），如果你的数据缺失比例很高又不太可能是完全随机的，GEE 的默认做法风险更大。第二是聚类数量：GEE 的稳健标准误是大样本性质，聚类数少于 30 到 40 时会明显偏小，需要用偏差校正。第三是相关结构本身是否是研究兴趣：如果组内相关系数（ICC）是要报告的核心指标，GLMM 会更自然，GEE 只能通过工作相关结构间接反映。

注意：GEE 的系数解释必须落在"平均"这个词上。写成"用药使患者症状改善的概率提高 X%"是错的，正确表述是"在总体水平上，用药组的结局概率平均比对照组高 X 个百分点"。这个措辞差异在很多评审场合是硬伤。

2. 数据准备与包选型，动手之前先把这几件事定下来

数据格式是 GEE 落地时第一个拦路虎。R 里绝大多数建模函数吃的是长格式（long format），也就是一个受试者有几条观测就占几行，同时需要一个 id 列标识归属、一个时间列标识次序。很多人手里的原始表是宽格式的，每个受试者一行，四次随访摊在y0到y3四列里，直接喂给geeglm会报错。转换用reshape或者tidyr::pivot_longer都行，关键是转完之后要检查 id 和时间列的取值是否符合预期。

选包这件事值得单独说，因为不同包对标准误的默认行为完全相反，这是新手踩坑最密集的地方。R 生态里能做 GEE 的主要有geepack、gee、geeM三个，语法相近但细节差得不少。geepack的geeglm在设计上模仿glm，公式写法和family参数都跟glm一致，学习成本最低，输出默认就是稳健标准误；gee包历史更久，参数名不太一样（corstr、id、Mv），默认输出的是模型型标准误（naive SE），必须显式传robust = TRUE才会给稳健版本，这个默认值坑过太多人。

2.1 长格式数据的构造与必查项

我习惯在做任何建模之前先跑一遍"体检三连"：每个体的观测条数分布、时间点的取值范围、结局的缺失模式。观测条数分布能告诉你数据是否平衡——平衡数据（每个人都有相同的随访次数）处理起来最简单，非平衡数据也不影响 GEE 运行，但会影响你对相关结构的假设是否合理。时间点必须是数值型且从 0 或者 1 开始连续编码，因为ar1相关结构是靠时间间隔算相关系数的，如果你的时间列是 "基线""第1月""第2月" 这种字符串，得先转成数字。

下面这段模拟代码构造了一个典型的随访数据集，200 个受试者、每人 4 次随访、二分类结局，处理组和对照组各半，同时带一个年龄协变量。数据生成时我特意加了个体随机截距，这样同一体内的观测天然相关，正好用来验证 GEE 的工作效果。

set.seed(2024) n_sub <- 200 n_visit <- 4 id <- rep(1:n_sub, each = n_visit) visit <- rep(0:(n_visit - 1), times = n_sub) trt <- rep(rbinom(n_sub, 1, 0.5), each = n_visit) age <- rep(round(rnorm(n_sub, 55, 10)), each = n_visit) # 个体随机截距，制造组内相关 b_i <- rep(rnorm(n_sub, 0, 0.8), each = n_visit) linpred <- -1.2 + 0.35 * trt - 0.15 * visit - 0.12 * trt * visit + 0.02 * (age - 55) + b_i y <- rbinom(n_sub * n_visit, 1, plogis(linpred)) dat <- data.frame(id, visit, trt, age, y) dat$id <- as.numeric(dat$id)

生成完先做检查：

table(table(dat$id)) # 每个人几条观测 range(dat$visit) # 时间编码范围 colSums(is.na(dat)) # 缺失情况 with(dat, tapply(y, id, mean)) |> summary() # 个体层面结局均值分布

最后那一行很关键。如果一个人四次测量全是 0 或全是 1，他对估计组内相关的贡献很有限，但样本里这类"全同个体"太多时，工作相关结构的估计会不稳定。这是实操里不常被提起但确实会影响收敛的细节。

2.2 三个主流包的差异对比

我的默认选择是geepack，理由很实际：它跟glm的接口几乎一致，broom::tidy()能直接收拾结果，MuMIn::QIC()也认它的对象。gee包我只在需要复现老文献结果时用，因为它当年是很多教材的示例包。geeM的优势是底层用 C++ 实现，聚类数上万、观测数几十万的数据跑起来明显更快，做大规模面板数据时值得考虑。

注意：换包对比结果时，一定要确认两边输出的都是稳健标准误。点估计应该一致到小数点后几位，标准误差异巨大基本都是默认设置不同造成的，不是模型本身有问题。

2.3 建模前必须敲定的五个问题

在敲代码之前，我会强迫自己回答五个问题，答不上来就先别跑模型。

第一个是结局类型。二分类用binomial(link = "logit")，计数用poisson(link = "log")，连续且近似正态用gaussian(link = "identity")，这都是常规操作。麻烦的是过度离散的计数数据，GEE 本身不估计尺度参数，泊松假设下如果方差远大于均值，稳健标准误虽然还能用，但效率会掉。这种情况下可以改quasi(variance = "mu", link = "log")，或者老实上负二项。有序多分类结局要用multgee包，geepack不支持。

第二个是聚类单位是什么。听起来很基础，但实操里经常搞错。临床试验里聚类单位是受试者，学校调查里可能是班级，纵向研究里是同一个人。聚类单位定错了，整个模型就废了。

第三个是时间趋势是不是要进模型。随访数据里结局随时间系统性变化非常常见，visit作为一个连续协变量放进去，或者转成因子放进去，看你想不想给时间点之间的差异单独建参数。前者更省自由度，后者更灵活。

第四个是处理与时间有没有交互。如果研究问题本身就关心"处理效应是否随时间衰减"，那交互项必须放，而且是核心参数。如果不关心，放进去只会增加解释成本。

第五个是聚类数量够不够。20 个聚类和 200 个聚类得到的稳健标准误质量完全不是一回事，前者需要偏差校正。

3. 工作相关结构怎么选，以及一次完整的模型跑通记录

相关结构的选择是 GEE 里最容易被神话的一环。很多教程会告诉你"选择与真实结构最接近的那个"，但真实结构你根本不知道。我更愿意用另一个思路：先固定一个通用结构跑主分析，再用多种结构做敏感性检查。如果结论在不同结构下保持一致，说明你的推断对相关结构的假设不敏感，这是比"选对了结构"更强的稳健性证据。

常用的四种结构各有偏好场景。independence假设同一体内观测完全独立，听起来最不合理，但它在三种情况下反而是好选择：聚类数很少（参数估计不稳）、随访次数极少（每人只有 2 次）、或者研究者只想拿到一个稳健的标准误而不关心结构本身。exchangeable假设同一体内任意两次观测的相关系数相同，是纵向数据里最常用的默认选项，因为它对个体内部的"平均相关性"做假设，代价只是多估一个参数。ar1假设相关性随时间间隔衰减，间隔越远相关性越弱，这在等间隔随访里很自然，但它要求数据结构是真正的有序序列。

unstructured给每两个时间点之间的相关系数单独估一个参数。四次随访意味着 6 个相关参数，聚类数只有几十的时候，这些参数会被估得非常飘。我一般只在聚类数上百、随访次数不超过 5 的时候考虑它，用它的主要目的是当作"最灵活的上界"，看看简化结构有没有丢掉重要信息。

3.1 四种相关结构的对比与选型口径

具体到参数个数，exchangeable无论随访多少次都只估一个相关系数，这就是它受欢迎的核心原因；unstructured是 k(k-1)/2 个，4 次随访 6 个，8 次随访直接飙到 28 个，聚类数不够根本撑不住。

有一个判断我用了很多次，效果不错：先把exchangeable和ar1都跑一遍，如果两个模型的点估计几乎重合、标准误差异在 10% 以内，说明数据结构没那么敏感，报exchangeable就行，解释起来最省事。如果差异明显，那就要去看 CIC 或者 QIC 的值，并且老实报告"不同结构下结论的稳定性"。这个做法比死磕"哪个结构是真"务实得多。

3.2 完整代码与逐步解读

数据准备好之后，主模型我一般这么写：

library(geepack) fit_ex <- geeglm( y ~ trt * visit + age, id = id, data = dat, family = binomial(link = "logit"), corstr = "exchangeable", std.err = "san.se" ) summary(fit_ex)

几个参数值得单独拎出来讲。id传的是聚类标识，geeglm内部会据此对数据做重排，但ar1结构下重排之后的时间顺序是否还正确就不能保证了，所以我习惯在建模前手动dat <- dat[order(dat$id, dat$visit), ]排一次。std.err有三个取值，san.se是 sandwich 估计（默认，最常用），jack是刀切法，j1s是一阶刀切，后两个计算量大但小样本表现更好，聚类数在 30 上下时我会考虑换jack。

模型跑完，结果长这样（数值来自上面的模拟数据，你的随机种子不同会有出入）：

Call: geeglm(formula = y ~ trt * visit + age, family = binomial(link = "logit"), data = dat, id = id, corstr = "exchangeable", std.err = "san.se") Estimate Std.err Wald Pr(>|W|) (Intercept) -1.18426 0.42318 7.831 0.00513 trt 0.31247 0.28764 1.180 0.27737 visit -0.14823 0.06172 5.769 0.01628 age 0.01985 0.00712 7.771 0.00530 trt:visit -0.10438 0.08716 1.434 0.23107 Estimated Scale Parameters: Estimate Std.err (Intercept) 1 0 Correlation structure = exchangeable Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.3247 0.0562

解读顺序我建议这么走。先看Estimated Correlation Parameters里的alpha，这是组内相关系数的估计值，0.32 说明同一个体的四次测量之间确实存在中等程度的相关，直接当独立样本处理会出问题——这一条也是说服合作者"该用 GEE"最有力的证据。接着看各系数的 Wald 检验，visit显著为负，说明随时间的推移结局概率下降；trt主效应不显著但方向为正。

这里要特别小心交互项的处理。trt:visit不显著并不意味着"处理没效果"，只意味着"处理效果不随时间变化"这个假设没有被数据推翻。在存在交互项的模型里，trt的系数含义是"在 visit = 0 时处理组与对照组的对数优势比"，不是"处理的总体效应"。很多人在这里读错，写出来的结论必然是错的。要得到处理的总体效应，得做边际预测或者用lincom类的线性组合检验，这个后面第 4 节会展开。

3.3 稳健标准误和模型型标准误的差距有多大

上面那张表输出的是稳健标准误，这是geepack的默认行为。为了让你直观感受"不处理相关性的代价"，我用普通的glm跑同一个公式做对比。glm假定 800 条观测完全独立，它的标准误就是模型型标准误，等价于independence结构下的 naive SE。

fit_glm <- glm(y ~ trt * visit + age, data = dat, family = binomial()) cbind( glm_naive = sqrt(diag(vcov(fit_glm))), gee_robust = sqrt(diag(vcov(fit_ex))) )

我实测下来，visit的模型型标准误通常只有稳健标准误的 60% 到 75%，也就是说 naive 的 p 值会被人为压小一大截。这个差距跟组内相关系数正相关：alpha越大、随访次数越多，差距越夸张。有一次组内相关接近 0.7 的数据，naive SE 只有稳健 SE 的一半，好几个在 GLM 里"极其显著"的变量换成 GEE 之后直接掉到 0.1 开外。

顺带说个容易混淆的点。fit_ex里那个Estimated Scale Parameters的估计值是 1，标准误是 0，这不是错误。二分类结局的方差完全由均值决定（Var = p(1-p)），没有额外的尺度参数可估，所以固定为 1。如果你看到这里不是 1，说明family选错了，比如拿gaussian跑了二分类数据。

4. 模型比较、边际效应与结果呈现

模型选好了不等于能直接写报告。审稿人和业务方想看的是"处理组比对照组高多少"，而不是"对数优势比是 0.31"。这一节把从系数到结论的转换过程补齐。

QIC（Quasi-likelihood under the Independence model Criterion）是 Pan 在 2001 年提出的，用来在 GEE 框架下做模型比较，思路类似 AIC。它有两个版本：QIC 用来比较相关结构（固定协变量，变动corstr），QICu 用来比较协变量组合（固定相关结构，变动右侧变量）。用错场景是常见错误，比如拿 QIC 去选变量，选出来的模型可能过拟合。

library(MuMIn) fit_ind <- geeglm(y ~ trt * visit + age, id = id, data = dat, family = binomial(), corstr = "independence") fit_ar1 <- geeglm(y ~ trt * visit + age, id = id, data = dat, family = binomial(), corstr = "ar1") QIC(fit_ind, fit_ex, fit_ar1)

输出里会给出 QIC、QICu 和 CIC 三个值。我的经验是，三个结构之间的 QIC 差距小于 2 到 4，就可以认为它们在实际意义上等价，这时候优先选参数少的exchangeable。差距很大（比如超过 10）才值得认真考虑换结构，同时要在报告里说明"敏感性分析显示结论对相关结构敏感"。另外 CIC 是 QIC 的一个变体，专门针对相关结构选择，有些文献更推荐它，两个一起看更稳妥。

4.1 从对数优势比到边际概率的转换

处理效应的呈现我一般给两个层面的数字：一个是优势比（OR）及其 95% 置信区间，一个是边际概率差（risk difference）。前者是统计报告的标准格式，后者是临床或业务方真正能理解的东西。

OR 直接指数化就行：

est <- coef(fit_ex) se <- sqrt(diag(vcov(fit_ex))) data.frame( term = names(est), OR = exp(est), lo = exp(est - 1.96 * se), hi = exp(est + 1.96 * se) )

边际概率差麻烦一点，因为geeglm对象的predict不支持se.fit，得自己做。最稳的办法是按聚类重抽样的 bootstrap——注意是按 id 抽，不是按行抽，同一受试者的四条记录必须整组进出，否则你就把组内相关结构人为打散了，估出来的区间会偏窄。

set.seed(123) newdat <- expand.grid(visit = 0:3, trt = 0:1, age = 55) boot_rd <- replicate(1000, { ids <- sample(unique(dat$id), replace = TRUE) idx <- unlist(lapply(ids, function(i) which(dat$id == i))) # 注意：重抽样后要重建唯一 id，否则 geeglm 会把重复个体识别为同一个体 bdat <- dat[idx, ] bdat$id <- rep(seq_along(ids), each = 4) bfit <- try(geeglm(y ~ trt * visit + age, id = id, data = bdat, family = binomial(), corstr = "exchangeable"), silent = TRUE) if (inherits(bfit, "try-error")) return(rep(NA, 8)) p <- predict(bfit, newdata = newdat, type = "response") p[2] - p[1] # 以 visit=0 为例的边际概率差 }) quantile(boot_rd, c(0.025, 0.5, 0.975), na.rm = TRUE)

这段代码里那个"重建唯一 id"的注释非常关键，我自己在这上面浪费过整整一个下午。bootstrap 重抽样会出现同一个人的数据被抽中多次，如果直接沿用原来的 id，geeglm会把这些重复行当成同一个体的重复测量，聚类结构完全错乱，结果是模型能跑通但数值毫无意义。这个坑不会报错，只会给你一个看起来正常但实际上错误的区间。

4.2 交互效应下的结果表达策略

有交互项的时候，我强烈建议用分层预测图代替系数表作为主结果。具体做法是把newdat扩展成处理组×对照组×四个时间点的完整网格，算出每组的预测概率和置信带，然后画两条线。

library(ggplot2) grid <- expand.grid(visit = 0:3, trt = 0:1, age = 55) pred <- predict(fit_ex, newdata = grid, type = "response") grid$prob <- as.numeric(pred) grid$grp <- factor(grid$trt, levels = c(0, 1), labels = c("对照组", "处理组")) ggplot(grid, aes(x = visit, y = prob, color = grp)) + geom_line(linewidth = 1) + geom_point(size = 2) + scale_y_continuous(labels = scales::percent) + labs(x = "随访时间点", y = "预测结局概率", color = NULL)

图的价值在于它把"交互项"这个抽象概念变成了肉眼可见的收敛或发散。如果两条线平行，交互不存在；如果处理组的线更平，说明处理缓冲了时间带来的下降。这种表达方式在给非统计背景的合作者看的时候，沟通效率比系数表高一个数量级。

注意：预测时协变量的取值必须明确交代。上例里年龄固定为 55 岁，报告里要写"以年龄 55 岁为例"，不能笼统说"预测概率"。这是审稿人最常挑的毛病之一。

5. 报错、不收敛与结果对不上的排查实录

GEE 的报错信息普遍不太友好，geeglm有时候只给你一句system is computationally singular或者干脆静默返回一堆NaN。这一节整理我实际踩过的坑，按出现频率排序。

5.1 常见报错与对应处理

第一条值得展开。system is computationally singular最常见的原因是放了一个"个体内恒定"的协变量，比如性别、出生年份这类不随时间变化的变量。理论上它们可以进 GEE（因为它们不是被 id 完全解释），但如果同时放了好几个这类变量，或者某个变量在样本里某一类只有极少数个体，就会撞上共线。排查办法是跑一遍glm看它的alias()结果，GLM 报出来的冗余关系 GEE 基本也会报。

第三条也很典型。unstructured结构加上 5 次随访加 40 个聚类，模型参数数直接逼近可用信息量，必然不收敛。这时候不要死磕，降结构是最务实的选择。我试过强行调maxit = 100让它"跑完"，结果确实是收敛了，但相关系数矩阵出现了大于 1 的值，完全不能用。

5.2 聚类数太少时的偏差校正

这是 GEE 最容易被忽视的软肋。sandwich 估计是大样本性质，聚类数在 30 以下时它会系统性低估标准误，导致假阳性。有个模拟研究做过对比，20 个聚类的情况下，名义 5% 的检验实际第一类错误率能到 10% 以上。

处理办法有两个。一是用geesmv包做 Mancl-DeRouen 偏差校正，或者用刀切法：

# 方案一：geepack 内建的刀切标准误 fit_jack <- geeglm(y ~ trt * visit + age, id = id, data = dat, family = binomial(), corstr = "exchangeable", std.err = "jack") # 方案二：geesmv 的偏差校正 # install.packages("geesmv") library(geesmv) GEE.var.md(y ~ trt * visit + age, id = "id", data = dat, family = binomial(), corstr = "exchangeable")

二是干脆换模型。聚类数少于 20 的时候，我会优先考虑 GLMM（lme4::glmer）或者带小样本校正的clubSandwich包配合geeglm。clubSandwich提供 CR2 类型的方差估计，在小样本下表现明显优于默认 sandwich，用法是在coeftest里传入vcovCR(fit, cluster = dat$id, type = "CR2")。这条路子我在聚类数只有十几个的项目里用过几次，标准误的变化很可观。

5.3 缺失数据：GEE 能吃多少

geeglm对缺失的处理比较直接：na.action默认是na.omit，只要某一行的任何一个变量缺失，整行被丢掉。这里有两个层面的风险要说清。

第一个层面是整条观测被删除。如果缺失发生在协变量上，丢一整行可能就把这个受试者剩下的观测也带走了，有效样本量损失比想象中大。建模前跑一遍colSums(is.na(dat))和sum(complete.cases(dat))，心里有个数。

第二个层面是缺失机制。GEE 的稳健性建立在完全随机缺失（MCAR）的假设上，也就是缺失跟任何已观测或未观测的变量都无关。这在真实数据里几乎不可能成立。数据如果是随机缺失（MAR），也就是缺失概率依赖于已观测变量，标准的 GEE 会偏，需要改用逆概率加权 GEE（ipw包）或者换用似然的 GLMM。

我的实践口径是：缺失比例低于 5%，MCAR 假设随便说说，直接跑；5% 到 20%，跑个敏感性分析，用加权 GEE 对比一下；超过 20%，别自己扛，先跟数据提供方确认缺失原因，再决定用哪种方法。这个分层不是统计学上的硬标准，但在实际项目里能帮你避开大部分麻烦。

6. 几个不太写在文档里但很影响结果的细节

第一个细节关于标准误的选择与聚类数的关系。同一批数据，std.err = "san.se"和std.err = "jack"给出的标准误在小样本下能差 20% 以上。我的经验是聚类数大于 50 用san.se就够了，30 到 50 之间跑一遍jack看看差异，小于 30 直接上偏差校正方法，别指望默认设置能给你可靠推断。

第二个细节关于**geepack输出的 Wald 统计量对应的自由度**。summary.geeglm给的是基于正态近似的 Wald 检验，没有像t分布那样的小样本修正。聚类数少的时候，这个检验偏激进。想看更保守的 p 值，可以用anova(fit_ex, test = "Wald")配合clubSandwich的方差估计，或者直接报置信区间让读者自己判断。我在报告里现在越来越少写具体 p 值，转而给"估计值（95% CI）"这种格式，争议少很多。

第三个细节是模型收敛状态的自检。geeglm不像glmer那样会明确告诉你"收敛了但是奇异拟合"，它只会在真的崩掉时才报警。所以我每次跑完都会检查三件事：fit$converged是否为 TRUE、相关参数估计是否在 (-1, 1) 区间内、fitted是否出现 0 或 1 这种极端值。有一回相关参数估到 0.998，模型没报任何错，但换一批数据一跑结果就翻车——这种"接近奇异"的状态最危险，因为你看不到任何警示。

第四个细节是我自己吃过亏的：不要用update()改corstr。直觉上update(fit_ex, corstr = "ar1")应该没问题，但实际测试中update有时候不会重新走一遍完整的 GEE 拟合流程，结果里的相关参数还是旧的。我现在一律写完整调用，宁可多打几行字。检验方法很简单，把update的结果和手动重跑的结果放在一起比summary()$corr，如果对不上就说明有问题。

第五个细节关于结果的可复现性。GEE 的点估计是确定性的，但只要涉及 bootstrap、jackknife 或者任何重抽样，就一定要在报告或脚本开头固定随机种子。我还会把sessionInfo()的输出连同结果一起存档，因为geepack不同版本之间在相关参数的初始值上做过调整，跨版本重跑可能出现小数点后几位的差异，虽然不影响结论，但在严格复现的场合会被追问。

最后一个心得关于什么时候该放弃 GEE。我见过不少人拿着只有 15 个聚类、12 个时间点、缺失率 35% 的数据硬要跑 GEE，然后在相关结构、标准误校正、缺失处理之间反复折腾。这种数据形态下，GEE 的渐近性质根本没有落脚点，换成 GLMM 或者干脆退回到描述性分析加敏感性讨论，是更负责任的选择。工具的价值在于匹配场景，把 GEE 用错地方的代价，比不用它更大。