Statsmodels GEE 完全指南：用广义估计方程处理面板、聚类与重复测量数据-平芜编程栈

数据分析
数据科学
科研

【免费下载链接】statsmodels

Statsmodels: statistical modeling and econometrics in Python

项目地址：https://gitcode.com/gh_mirrors/st/statsmodels

本文基于仓库文档 docs/source/gee.rst，并结合 statsmodels/genmod/generalized_estimating_equations.py、statsmodels/genmod/cov_struct.py 与 statsmodels/genmod/qif.py 等源码与测试展开。

广义估计方程（Generalized Estimating Equations, GEE）是 statsmodels 中处理面板数据、聚类数据与重复测量数据的核心边际回归模型：它允许同一簇（cluster）内部的观测彼此相关，同时假设不同簇之间相互独立。本文将从理论定位、真实数据集上的实战示例、三类模型类、十种依赖结构、结果对象与推断方法（稳健方差、Score 检验、QIC）以及 QIF 备选方案等维度，完整梳理 GEE 在 statsmodels 中的使用方法与底层实现，读者读完即可在自己的纵向/聚类数据上完成建模、拟合、诊断与模型比较。

一、什么是 GEE：边际模型如何建模“组内相关”

GEE 由 Liang 与 Zeger 于 1986 年提出，是广义线性模型（GLM）在依赖数据（dependent data）上的延伸。它估计的是边际（marginal）回归模型：描述响应变量均值如何随协变量变化，而不对簇内的随机效应做完整建模，簇内相关性通过一个“工作相关结构”（working correlation/covariance structure）来描述。

关键适用场景：

面板数据（panel data）：同一实体（如国家、公司、个体）多个时间点的观测；
聚类数据（cluster data）：同一学校、医院、家庭内的多个观测；
重复测量数据（repeated measures data）：同一受试者多次随访的测量值。

在统计假设上，GEE 只要求：

同一簇内的观测可能相关；
不同簇之间的观测相互独立。

它支持与 GLM 完全相同的单参数指数族分布（one-parameter exponential families），如高斯、泊松、二项、负二项、伽马、逆高斯等。正因为“工作相关结构”即使设定错误，回归系数的点估计依然一致（一致性不依赖相关结构），GEE 在生物统计与计量经济学中被广泛使用。

二、快速上手：癫痫发作数据的 Poisson GEE

文档 docs/source/gee.rst 给出了一个可直接运行的示例：使用 MASS 包的epil（癫痫）数据集，拟合以subject（患者）为簇、簇内采用可交换（Exchangeable）相关结构的 Poisson 回归。

import statsmodels.api as sm import statsmodels.formula.api as smf data = sm.datasets.get_rdataset('epil', package='MASS').data fam = sm.families.Poisson() ind = sm.cov_struct.Exchangeable() mod = smf.gee("y ~ age + trt + base", "subject", data, cov_struct=ind, family=fam) res = mod.fit() print(res.summary())

要点拆解：

代码片段	作用
`smf.gee(formula, groups, data, ...)`	通过公式接口创建 GEE 模型，`groups`指定簇标识列名
`sm.families.Poisson()`	指定分布族（计数数据常用）
`sm.cov_struct.Exchangeable()`	指定簇内工作相关结构：任意两个簇内观测相关性相同
`res.summary()`	输出系数、标准误、z 值与 p 值等回归表

关于本示例与真实癫痫数据，仓库测试 statsmodels/genmod/tests/test_gee.py 中的test_poisson_epil给出了可复现版本：使用GEE.from_formula("y ~ age + trt + base", data["subject"], data, cov_struct=ind, family=fam)拟合，并在cov_type="naive"下断言 GEE 的系数与标准误与普通 Poisson GLM（GLM.from_formula）一致（rtol=1e-6）。这印证了一个重要性质：当相关结构退化为独立时，GEE 的均值参数估计与 GLM 一致。

三、模型类：GEE、NominalGEE 与 OrdinalGEE

模块statsmodels.genmod.generalized_estimating_equations（从statsmodels.api与statsmodels.formula.api均可访问，见 statsmodels/api.py 与 statsmodels/genmod/api.py）提供了三个模型类：

类	适用响应变量	默认分布族	默认相关结构
`GEE`	连续/计数/二值等	`Gaussian`	`Independence`
`OrdinalGEE`	有序分类	`Binomial`	`OrdinalIndependence`
`NominalGEE`	无序分类	`_Multinomial`	`NominalIndependence`

3.1 GEE 构造函数与核心参数

GEE直接继承自GLM（generalized_estimating_equations.py），构造函数签名如下：

GEE(endog, exog, groups, time=None, family=None, cov_struct=None, missing="none", offset=None, exposure=None, dep_data=None, constraint=None, update_dep=True, weights=None, **kwargs)

各参数含义（依据源码 docstring）：

endog/exog：因变量与自变量；
groups：簇标签数组，观测按此分组，簇间要求独立；
time：观测的时间/位置信息，用于依赖距离的相关结构（如 AR、Stationary、Unstructured）；若不提供，默认取簇内等间隔的 0,1,2,... 网格（L653-662）；
family：GLM 分布族实例，缺省为Gaussian；若传入非families.Family子类实例会抛出ValueError（L599-603）；
cov_struct：工作相关结构实例，缺省为Independence（L606-613）；
missing：缺失值处理方式，'none'/'drop'/'raise'；
offset/exposure：线性预测项中的偏移量与暴露量（当链接为 log 时，log(exposure)会加入 offset）；
dep_data：供Nested等结构使用的依赖数据；
constraint：形如(lhs, rhs)的线性等式约束，满足lhs * params = rhs；
update_dep：是否在迭代中更新相关参数；当所有簇都只有 1 个观测（等价于拟合 GLM）时会自动置为False（L686-690）；
weights：簇内观测权重。

3.2 from_formula 公式接口

GEE.from_formula(formula, groups, data, subset=None, time=None, offset=None, exposure=None, *args, **kwargs)是文档示例使用的入口。它的特殊之处在于：groups、time、offset、exposure、dep_data均支持传字符串形式的列名，内部会从data中取值（L771-794）。例如文档示例中"subject"就是数据框中簇标识列名。

3.3 OrdinalGEE 与 NominalGEE：分类响应的重编码

有序与无序分类响应无法直接套用高斯/泊松框架，statsmodels 的做法是将多分类数据重编码为一组二值指示变量再走 GEE 流程：

OrdinalGEE.setup_ordinal将每个有序观测展开为ncut（水平数减 1）个累积指示变量I(y > 阈值)，并额外为每个阈值生成截距列（列名形如I(y>1.0)），要求使用Binomial族（L2468-2589）；
NominalGEE.setup_nominal则用np.kron构造分块对角设计矩阵，为每个非基准水平生成独立系数块（列名形如x1[1.0]），默认使用_Multinomial族（L2798-2927）。

两者对应的相关结构是OrdinalIndependence与NominalIndependence（详见下一节）。仓库中 statsmodels/genmod/tests/gee_categorical_simulation_check.py 专门对这两类模型做仿真检验。

四、依赖（相关）结构：cov_struct 模块全景

GEE 的核心在于“工作相关结构”。所有结构均定义在 statsmodels/genmod/cov_struct.py 中，统一继承自基类CovStruct。基类通过dep_params属性保存当前相关参数，并约定四个核心接口：initialize(model)、update(params)、covariance_matrix(endog_expval, index)、covariance_matrix_solve(...)。

文档 docs/source/gee.rst 列出的十种结构及源码要点如下：

结构类	描述	关键实现细节
`Independence`	独立结构	相关矩阵恒为单位阵，`update`不做任何事（L208-235）
`Exchangeable`	可交换结构	簇内任意两观测相关性相同，`dep_params`为标量相关；用矩估计（残差平方和）更新（L330-413）
`Unstructured`	无结构	每个位置对都有独立相关参数；必须提供整数类型的`time`，否则抛`ValueError`（L260-268）
`Autoregressive`	一阶自回归	相关性 =`dep_params ** 距离`；`dist_func`自定义距离函数，`grid=True`走网格加速实现；默认`time`为簇内索引位置（L792-863）
`Stationary`	平稳结构	相关性是观测间距的任意函数，`max_lag`控制纳入的最大距离；`grid`参数决定用索引位置还是`time`定义距离（L620-789）
`Nested`	嵌套结构	刻画层级分区（如学校→班级→学生）每层独立的随机效应方差；层级由`dep_data`定义，公式形式建议`0 + a + b + ...`（L416-448）
`GlobalOddsRatio`	全局比值比	面向有序/无序分类数据的相关结构，`endog_type`取`"ordinal"`或`"nominal"`，用粗比值比（crude odds ratio）初始化（L1087-1161）
`NominalIndependence`	名义独立	名义分类模型专用（L1354）
`Equivalence`	等价结构	基于配对定义等价关系的相关结构，支持`pairs`/`labels`构造（L1383）

4.1 关于 PSD 投影与数值稳健性

CovStruct.covariance_matrix_solve负责求解形如covmat * soln = rhs的线性方程组。当相关矩阵不是半正定（PSD）时，基类会用cov_nearest将其投影到最近的 PSD 矩阵；投影方法由构造参数cov_nearest_method控制，取值"clipped"（默认）或"nearest"（L46-49、L160-198）。若 20 次投影仍无法分解，会回退为对角矩阵并发出ConvergenceWarning。Exchangeable与Stationary(grid=True)等结构重写了covariance_matrix_solve，利用结构的解析性质避免显式构造矩阵，大幅提升大簇计算效率。

4.2 结构选择的实践建议

簇内相关性大致恒定 →Exchangeable；
时间序列式重复测量、相关性随间隔衰减 →Autoregressive/Stationary；
观测位置不完全等间隔 → 为Autoregressive提供time并自定义dist_func；
数据存在多级嵌套 →Nested+dep_data；
有序/无序分类 →OrdinalIndependence/NominalIndependence/GlobalOddsRatio；
不确定相关形式时，Independence仍能给出一致的均值参数估计（配合稳健方差）。

五、分布族与链接函数：与 GLM 完全一致

GEE 支持的分布族与 GLM 相同，文档明确指出“当前实现的分布族与 GLM 相同”，完整列表见 docs/source/glm.rst 的家族清单，实现位于 statsmodels/genmod/families/，包括：

Gaussian（高斯）
Binomial（二项）
Poisson（泊松）
NegativeBinomial（负二项）
Gamma（伽马）
InverseGaussian（逆高斯）
Tweedie（Tweedie，含var_power参数）

链接函数同样复用 GLM 的链接体系。并非每个链接都适用于每个分布族，可用链接列表可通过以下方式查询（文档原句）：

>>> sm.families.family.<familyname>.links

例如sm.families.Poisson().links会列出泊松族可用的链接（如log、identity、sqrt）。源码中GEE.__init__会检查传入族实例的链接是否在该族的safe_links内，若不满足会发出DomainWarning（L525-537）。

六、拟合与推断：fit 参数、稳健方差与偏差校正

6.1 fit 方法与收敛控制

GEE.fit的签名（L1295-1307）：

fit(maxiter=60, ctol=1e-6, start_params=None, params_niter=1, first_dep_update=0, cov_type="robust", ddof_scale=None, scaling_factor=1.0, scale=None)

maxiter：最大迭代次数，默认 60；达到上限未收敛会发IterationLimitWarning；
ctol：收敛阈值，基于回归参数更新向量的 L2 范数（del_params = sqrt(sum(score**2))）；
params_niter/first_dep_update：控制相关参数更新的频率与起始时机——只有同时满足itr % params_niter == 0且itr >= first_dep_update才更新依赖结构，且至少完成一次相关参数更新后才允许提前收敛（L1360-1372）；
cov_type：可选"robust"（默认，三明治稳健协方差）、"naive"（模型假定方差）、"bias_reduced"（Mancl–DeRouen 小样本偏差校正）；
ddof_scale：估计尺度参数时从样本量中减去的自由度，默认取exog列数；
scale：若为数值则直接固定尺度参数；None时对Binomial、Poisson、NegativeBinomial等族默认取 1.0（estimate_scale见 L981-1026）。

拟合循环的核心是交替执行两步：_update_mean_params求解当前相关结构下的拟得分方程更新均值参数，_update_assoc基于残差更新相关参数。二者在源码 L1081-1136 与 L1337-1372 中实现。

6.2 三明治协方差：robust / naive / bias_reduced

_covmat（L1166-1233）同时计算两个协方差：

naive（model-based）：cov_naive = B^{-1} * scale，仅在相关结构正确设定时有效；
robust（sandwich）：cov_robust = B^{-1} * C * B^{-1}，其中C是各簇得分向量的外积之和，即使工作相关结构设定错误依然渐进有效，这也是 GEE 的默认选择。

_bc_covmat（L1237+）实现 Mancl 与 DeRouen（2001）提出的偏差校正三明治估计，用于小簇数场景，对应cov_type="bias_reduced"。

GEEResults中还提供了standard_errors(cov_type=...)便捷方法，可在三种协方差间任意切换（L1946-1977）。

6.3 假设检验：Score 检验与模型比较

GEE.compare_score_test(submodel)：给定一个已拟合的子模型，对大模型做 Score 检验，返回{"statistic", "p-value", "df"}字典，无需对大模型调用 fit（L834-979）；
GEEResults.score_test()：针对带线性约束（constraint）拟合的模型返回 Score 检验结果；与compare_score_test互为补充——后者适合比较两个显式模型，前者支持任意线性等式约束（L1984-2013）。

两个方法均以 Guo 与 Pan（2002）关于 GEE Score 检验小样本表现的研究为依据。

6.4 模型选择：QIC 与 QICu

GEE 没有传统似然，模型比较使用 Pan（2001）提出的拟信息准则 QIC。GEE.qic（L1798-1883）通过数值积分 Wedderburn 拟似然，返回三个量：

ql：拟似然值；
qic：可同时比较均值结构与相关结构；
qicu：简化的 QIC，只能比较均值结构。

结果对象的res.qic(scale=..., n_step=1000)是便捷入口（L2052-2076）。源码特别提示两点（L1827-1842）：

数值拟似然与其它软件用解析式算出的 QIC 绝对值不同，只有同一数据上不同模型的 QIC 之差才有意义；
当尺度参数未知时，比较模型应使用同一个scale 估计值。

6.5 正则化 GEE

GEE.fit_regularized(pen_wt, scad_param=3.7, maxiter=100, ddof_scale=None, update_assoc=5, ctol=1e-5, ztol=1e-3, eps=1e-6, scale=None)提供 SCAD 惩罚的稀疏估计，面向高维纵向数据（Wang, Zhou & Qu, 2012），要求使用正则链接（L1538-1661）。测试 statsmodels/genmod/tests/test_gee.py 中的test_regularized_poisson与test_regularized_gaussian覆盖了该路径。

6.6 边际效应与其他结果方法

GEEResults继承自GLMResults，因此具备残差（resid、resid_split、resid_centered）、get_margeff边际效应、plot_added_variable/plot_partial_residuals/plot_ceres_residuals诊断图等能力；GEEMargins（L3258）封装边际效应的推断结果。簇级残差按簇拆分的方法resid_split/resid_centered_split由test_gee_results_resid_split验证。

七、备选方案：二次推断函数 QIF

模块 statsmodels/genmod/qif.py 提供 QIF（Quadratic Inference Functions）作为 GEE 的高效替代。QIF 的优势在于：对相关结构设定更稳健、可能更高效，并提供不同于 GEE 的模型选择与推断途径（L117-148）。其理论依据为 Qu, Lindsay & Li（2000）发表在Biometrika的论文。

使用方式与 GEE 对称：

from statsmodels.genmod.qif import QIF from statsmodels.genmod.families import Poisson from statsmodels.genmod import cov_struct model = QIF(endog, exog, groups, family=Poisson(), cov_struct=cov_struct.Exchangeable()) res = model.fit()

QIF 配套三种协方差结构（继承自QIFCovariance）：QIFIndependence、QIFExchangeable、QIFAutoregressive，结果类为QIFResults（含aic、bic、fittedvalues、summary）。模型类的fit(maxiter=100, start_params=None, tol=1e-6, gtol=1e-4, ddof_scale=None)与from_formula(formula, groups, data, ...)均已实现；测试见 statsmodels/genmod/tests/test_qif.py，其test_qif_numdiff用数值微分验证了 QIF 目标函数的梯度。

八、模块参考速查表

按 docs/source/gee.rst 的 Module Reference 整理：

类别	符号	所在源码文件
模型类	`GEE`、`NominalGEE`、`OrdinalGEE`	statsmodels/genmod/generalized_estimating_equations.py
模型类	`QIF`	statsmodels/genmod/qif.py
结果类	`GEEResults`、`GEEMargins`	generalized_estimating_equations.py
结果类	`QIFResults`	statsmodels/genmod/qif.py
依赖结构	`CovStruct`、`Autoregressive`、`Equivalence`、`Exchangeable`、`GlobalOddsRatio`、`Independence`、`Nested`、`NominalIndependence`、`Stationary`、`Unstructured`	statsmodels/genmod/cov_struct.py
分布族	同 GLM（Gaussian/Binomial/Poisson/NegativeBinomial/Gamma/InverseGaussian/Tweedie）	statsmodels/genmod/families/
链接函数	同 GLM，可用`sm.families.family.<name>.links`查询	statsmodels/genmod/families/links.py

九、进一步阅读

核心实现：statsmodels/genmod/generalized_estimating_equations.py（3530 行，含 GEE/OrdinalGEE/NominalGEE 完整实现）
依赖结构：statsmodels/genmod/cov_struct.py
QIF 实现：statsmodels/genmod/qif.py
测试与验证：statsmodels/genmod/tests/test_gee.py、statsmodels/genmod/tests/test_qif.py、statsmodels/genmod/tests/gee_simulation_check.py（仿真恢复真实参数检验）
分布族与链接：statsmodels/genmod/generalized_linear_model.py 及 docs/source/glm.rst

十、参考文献（与文档一致）

Liang KY, Zeger SL. "Longitudinal data analysis using generalized linear models".Biometrika(1986) 73(1): 13-22.
Zeger SL, Liang KY. "Longitudinal Data Analysis for Discrete and Continuous Outcomes".BiometricsVol. 42, No. 1 (Mar., 1986), pp. 121-130.
Rotnitzky A, Jewell NP (1990). "Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data",Biometrika, 77, 485-497.
Guo X, Pan W (2002). "Small sample performance of the score test in GEE".
Mancl LA, DeRouen TA (2001). "A covariance estimator for GEE with improved small-sample properties".Biometrics57(1): 126-134.
Pan W (2001). "Akaike's information criterion in generalized estimating equations".Biometrics57(1).
Qu A, Lindsay B, Li B (2000). "Improving Generalized Estimating Equations using Quadratic Inference Functions",Biometrika87:4.
Heagerty PJ, Zeger SL (1996). "Marginal Regression Models for Clustered Ordinal Measurements".JASA91(435).
Wang L, Zhou J, Qu A (2012). "Penalized generalized estimating equations for high-dimensional longitudinal data analysis".Biometrics68(2): 353-360.