1. 纵向数据的相关性陷阱以及广义估计方程到底在算什么做过随访数据的人大概都有过这种经历同一个人测了四次直接把这四行当成四个独立样本扔进 logistic 回归结果出来一堆显著变量效应量看着还挺漂亮等到审稿意见回来一句你是否考虑了同一个体内观测的相关性前面的功夫基本白费。**广义估计方程Generalized Estimating EquationsGEE**就是专门用来收拾这类数据的工具而 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 三个主流包的差异对比维度geepackgeegeeM主函数geeglm()gee()geem()公式写法与glm一致与glm一致与glm一致默认标准误稳健sandwich模型型naive稳健支持ar1支持需按 id 与时间排序支持支持小样本校正需配合geesmv无内置无内置与emmeans兼容较好一般一般输出对象可update是是是我的默认选择是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 四种相关结构的对比与选型口径结构参数个数4 次随访适用场景主要风险independence0聚类数少、组内相关弱效率损失但推断仍稳exchangeable1纵向数据的默认首选相关性随时间衰减时效率下降ar11等间隔随访、相关性随时间衰减时间编码不连续时报错unstructured6时点少、聚类数充足参数过多导致不收敛m-dependent按阶数只在相邻时点间有相关阶数选择缺乏依据具体到参数个数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。这一节把从系数到结论的转换过程补齐。QICQuasi-likelihood under the Independence model Criterion是 Pan 在 2001 年提出的用来在 GEE 框架下做模型比较思路类似 AIC。它有两个版本QIC 用来比较相关结构固定协变量变动corstrQICu 用来比较协变量组合固定相关结构变动右侧变量。用错场景是常见错误比如拿 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] # 以 visit0 为例的边际概率差 }) quantile(boot_rd, c(0.025, 0.5, 0.975), na.rm TRUE)这段代码里那个重建唯一 id的注释非常关键我自己在这上面浪费过整整一个下午。bootstrap 重抽样会出现同一个人的数据被抽中多次如果直接沿用原来的 idgeeglm会把这些重复行当成同一个体的重复测量聚类结构完全错乱结果是模型能跑通但数值毫无意义。这个坑不会报错只会给你一个看起来正常但实际上错误的区间。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设计矩阵共线或某协变量在个体内不变且与 id 高度关联检查alias()删掉冗余变量NA出现在标准误列该变量在部分聚类内取值恒定sandwich 矩阵奇异合并类别或改用jack标准误迭代不收敛maxit警告相关结构过于复杂或聚类数太少降级到independence/exchangeablear1报时间间隔错误时间列未按 id 排序或编码不连续手动order(id, time)重编码结果与 GLMM 差异巨大边际 vs 条件效应属正常现象用同口径对比不要混用 OR稳健 SE 异常小聚类数太少 30用geesmv做偏差校正第一条值得展开。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 的时候我会优先考虑 GLMMlme4::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 会偏需要改用逆概率加权 GEEipw包或者换用似然的 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 用错地方的代价比不用它更大。