形式陈述
对同一个人连续记录两次是否出现某事件,两个回答可能相关,却仍可用 logit 描述各次事件概率。广义线性模型 理路 广义线性模型 Generalized linear model · GLM 用链接函数把指数族响应的条件均值连接到协变量线性预测子。 提供链接、均值导数和工作方差;本页接着构造一个允许簇内相关的向量估计方程。只固定工作相关矩阵,不要求给出簇内响应的完整联合似然,也不把相关参数估计混入本页的证明。
固定的是相关规则,不是全部权重
设完整簇资料 W g = ( X g , Y g ) 独立同分布,g = 1 , … , G ,参数维数 p 固定。Y g ∈ R m g ,1 ≤ m g ≤ m ,上界 m 不随 G 增长。为简化记号,X g 包含整个簇的设计、簇大小及已知暴露量等输入。以下所有条件均值都条件于这份完整设计:
E ( Y g ∣ X g ) = μ g ( β 0 ) , D g ( β ) = ∂ μ g ( β ) ∂ β T ∈ R m g × p . 典型均值为 μ g i ( β ) = ℓ − 1 ( x g i T β + o g i ) ;o g i 是已知 offset。设
A g ( β ) = diag { v g i ( β ) } , V g ( β ) = A g ( β ) 1 / 2 R g A g ( β ) 1 / 2 , 其中各工作方差 v g i > 0 是已规定的设计与参数函数;R g = R ( X g ) 是预先指定、对称正定且对角线为1的工作相关规则。这里使用同一个已固定规则,以保留完整簇的 IID 设定;不从本批 Y 中调出 R 或新的分散参数。R g 不随 β 变化,但 A g 通常随 β 变化,故整个 V g 并未被冻结。
令 r g = Y g − μ g 。簇得分、总方程和平均方程分别为
(1) s g ( β ) = D g T V g − 1 r g , U G ( β ) = ∑ g s g ( β ) , Ψ G ( β ) = U G ( β ) / G . GEE 估计量是 Ψ G ( β ) = 0 的已指定根,或具有已报告残差的近似根;这正是Z-估计的求根接口 理路 估计方程与 Z 估计 Z-estimation · Estimating equation 区分求根和极小化,给出多根定位、近似残差与局部 Jacobian 线性化的完整接口。 。输入必须包括簇标签、设计与响应、链接、工作方差、每簇 R g 、初值、停止容差与最大迭代数。输出不能只有系数,还应包括拟合均值、平均得分残差和下面两种协方差。
在返回的 β ^ 处计算求和矩阵
(2) H = ∑ g D g T V g − 1 D g , M = ∑ g s g s g T , C ^ = H − 1 M H − 1 . C ^ 已是 β ^ 的协方差估计,不能再除以 G 或总行数。按sandwich 的平均约定 理路 Sandwich 协方差估计 Sandwich covariance estimation 从得分与敏感度的样本矩阵算出可复核的协方差、对比标准误和异方差稳健推断。 ,两侧敏感度是 − H / G ,中间矩阵是 M / G ,最后再除 G ,恰好得到式(2)。H − 1 则是工作模型协方差;要使用它,还需保证 B = Q ;真实条件协方差等于工作协方差是一个充分条件,所得公式通常也只是渐近公式。
为什么错误的工作相关不一定移动目标
在真参数处,D g 和 V g 都由 X g 决定。只要相应期望存在,
(3) E { s g ( β 0 ) ∣ X g } = D g T V g − 1 { E ( Y g ∣ X g ) − μ g ( β 0 ) } = 0. 这个等式不使用 R g 等于真实相关矩阵。它解释均值中心为什么可以保留,却没有单独证明总体根唯一、算法选对根或估计量一致;这些仍是后面的独立条件。式(3)也不是说拟合后的 s g ( β ^ ) 条件均值为零,因为 β ^ 已经使用了响应。
条件于完整 X g 很重要。仅知 E ( Y g i ∣ x g i ) = μ g i ,不能随意把它升级为 E ( Y g i ∣ X g ) = μ g i 。例如第二次协变量包含第一次结局时,完整设计已经透露第一次结局的信息;而 V g − 1 还会把不同时间的残差混合。这样的反馈资料需要另外验证矩条件,不能靠“边际均值写对了”一句话跳过。
期望敏感度与观测 Jacobian 的差别
记 K g ( β ) = D g ( β ) T V g ( β ) − 1 、h g = K g D g 。逐坐标求导得到
(4) ∂ j s g = ( ∂ j K g ) r g − K g D g , ⋅ j . 第一项包含均值二阶导数以及工作方差的导数。由式(3)使用的完整条件均值,在 β 0 处对式(4)取条件期望,第一项消失。因此在允许微分与期望交换的可积条件下,
(5) − E { ∂ β s g ( β 0 ) } = E { h g ( β 0 ) } =: Q . 有限样本中却有 − ∂ β U G = H − T G ,其中 T G 的第 j 列是 ∑ g ( ∂ j K g ) r g 。即使 U G ( β ^ ) = 0 ,这些新的残差加权和也未必为零。H 是把期望敏感度结构代入当前参数后求和的 bread,并不等于每次样本的观测负 Jacobian;后者甚至可以不对称。两者在下述条件下除以 G 都趋于 Q ,这才解释其一阶等价。
一套可以逐项检查的充分条件
为了把“均值正确”与真正的抽样结论接起来,固定一个紧凸参数集 Θ ,β 0 在其内部;所有函数联合可测。除前述 IID、有界簇大小与式(3)之外,要求:
在包含 Θ 的开集,μ g 二次连续可微,A g 一次连续可微并保持正对角元,因而 s g 连续可微。写 J g = ∂ β s g ,存在随机界 L s , L J ,使 sup Θ ‖ J g ‖ ≤ L s ,且 ‖ J g ( b ) − J g ( c ) ‖ ≤ L J ‖ b − c ‖ ;E L s 2 < ∞ 、E L J < ∞ 、E ‖ s g ( β 0 ) ‖ 2 < ∞ 。
E ‖ h g ( β 0 ) ‖ < ∞ ,且 ‖ h g ( b ) − h g ( c ) ‖ ≤ L h ‖ b − c ‖ ,其中 E L h < ∞ 。Q = E h g ( β 0 ) 正定。
总体方程 Ψ ( b ) = E s g ( b ) 在 Θ 上只有分离的目标根:对每个 ε > 0 ,inf b ∈ Θ : ‖ b − β 0 ‖ ≥ ε ‖ Ψ ( b ) ‖ > 0 ,空集的下确界约定为正无穷。求解器返回 Θ 中的可测选择,并满足 ‖ Ψ G ( β ^ ) ‖ = o P ( G − 1 / 2 ) 。
这些是充分条件,不声称必要。比如设计有界、暴露量上下有正的有限界、参数限制在固定紧集、R g 最小特征值统一远离零,且 E ‖ Y g ‖ 2 < ∞ 时,本页的 logit 与 log-link 公式在紧集上具有所需光滑包络;但仍需另外核对总体分离和 Q 正定。不能把一次迭代的可逆 H 当成全部总体条件。
记 B = E { s g ( β 0 ) s g ( β 0 ) T } 。则
(6) G ( β ^ − β 0 ) = Q − 1 1 G ∑ g s g ( β 0 ) + o P ( 1 ) ⇒ N p ( 0 , Q − 1 B Q − 1 ) , (7) G C ^ ⟶ P Q − 1 B Q − 1 . 证明并不另外发明一套 GEE 渐近理论。先把紧参数集上的一致大数律 理路 参数化一致大数律 Parametric uniform law of large numbers 用紧参数集和可积随机 Lipschitz 界,把逐点大数律提升为同一样本上的全参数控制。 逐分量用于 s g :L s 给随机 Lipschitz 界,s g ( β 0 ) 给可积锚点。结合总体分离和趋零残差,得到一致定位。再将一致大数律用于 J g ,其锚点由 L s 控制,Lipschitz 界为 L J ;支配收敛同时允许交换微分与期望,使总体 Jacobian 在真值处等于 − Q 。已有Z 定理的积分 Jacobian 展开 理路 估计方程与 Z 估计 Z-estimation · Estimating equation 区分求根和极小化,给出多根定位、近似残差与局部 Jacobian 线性化的完整接口。 于是给出式(6),其中独立样本量是簇数 G 。
协方差的插件也需检查:‖ s g ( β ^ ) − s g ( β 0 ) ‖ ≤ L s ‖ β ^ − β 0 ‖ 。记 d = ‖ β ^ − β 0 ‖ ,外积差的平均至多为
2 d P G { L s ‖ s g ( β 0 ) ‖ } + d 2 P G L s 2 = o P ( 1 ) . Cauchy–Schwarz 与所给二阶矩使两个平均稳定,所以 M / G → P B ;同样,H / G 与 P G h g ( β 0 ) 的差不超过 d P G L h = o P ( 1 ) ,故 H / G → P Q 。最后用矩阵求逆连续性得到式(7)。这避免把共同使用拟合参数的残差误当作 IID 样本直接套大数律。
若真实条件协方差是 Σ g ,则
B = E ( D g T V g − 1 Σ g V g − 1 D g ) . 若进一步假定 Σ g = V g ( β 0 ) ,则 B = Q ,使式(6)的方差化为 Q − 1 ,相应插件是 H − 1 。对满足 c T Q − 1 B Q − 1 c > 0 的固定对比,c T C ^ c 支持渐近正态学生化。有限样本若 H 奇异或对比方差为零,应报告推断失败;为陈述渐近随机序列所作的任意有限值延拓,不是程序的成功输出。
直觉
每簇先把“实际回答减去预测均值”变成残差向量,再用工作协方差重新组合这些残差,最后由 D g T 投到参数方向。工作相关改变同一人的两次回答怎样共同推动系数,也可能改变一次样本的根。只要这个组合在给定完整设计后不偷看响应,真参数处的平均推力仍为零。
图片加载失败 GEE:均值中心、工作相关与两种协方差 均值正确保护中心,很多独立簇提供重复抽样,sandwich 校准中心附近的波动。这三件事缺一不可。图中 H − 1 的捷径需要额外的方差依据,例如真实协方差等于工作协方差;它不是因为矩阵更容易计算就自动成立。
如何实际求根
在当前 β 计算 U G , H ,解 H δ = U G ,再作
(8) β new = β + α δ , 0 < α ≤ 1. 正号来自 U G ( β + δ ) ≈ U G ( β ) − H δ 。这是期望敏感度 scoring,不是直接使用观测 Jacobian 的 Newton 步。固定 R 之后,每一步仍须更新 μ , D , A 。原始论文的印刷式(8)带负号,与其残差 Y − μ 的定义不相合;这里按所声明的得分通过线性化确定正号,而不照抄印刷符号。
可以缩小 α ,检查平均得分范数是否下降、均值是否有效、矩阵是否正定。但 GEE 未必是某个联合对数似然的梯度,故不能把这一步称为联合似然必然上升;残差回溯也不保证任意初值都能到达所需根。若方向未通过回溯、达到最大迭代数、发生溢出、二元均值数值饱和或 H 奇异,应明确退出并保留诊断。
停止时同时报告 ‖ U G ‖ / G 、最后步长或 ‖ H − 1 U G ‖ 、矩阵条件诊断及根选择规则。固定的 10 − 12 是当前算术精度,不是沿无限样本序列的 o P ( G − 1 / 2 ) 定理;小裸残差也不能排除分离方向上曲率一起消失。
若每簇使用稠密线性代数,一次迭代可按 O { ∑ g ( m g 3 + m g 2 p + m g p 2 ) + p 3 } 计费,再乘迭代次数。用 V g = A g 1 / 2 R g A g 1 / 2 可预先分解固定 R g ,将 m g 3 的成本移到初始化;求解所需缩放仍每次更新。最终还要形成 M ,额外约 O ( G p 2 ) 并完成夹乘。通过分解后的线性方程实现计算,避免显式求每个 V g − 1 。
例子与边界
两次二元回答,从第一步到两个最终根
四个独立簇各有两个时点。按簇依次给出 ( x g , y g ) :
( ( 0 , 1 ) , ( 0 , 0 ) ) , ( ( 0 , 2 ) , ( 0 , 1 ) ) , ( ( 1 , 2 ) , ( 1 , 1 ) ) , ( ( 0 , 1 ) , ( 1 , 0 ) ) . 模型为 logit μ g i = β 0 + β 1 x g i ,含截距的设计记作 Z g 。于是 A g = diag { μ g i ( 1 − μ g i ) } 、D g = A g Z g 。两次计算都从 ( 0 , 0 ) 开始,分别固定
R 0 = I 2 , R 1 / 4 = ( 1 1 / 4 1 / 4 1 ) . 在初值处 μ = 1 / 2 、A = I / 4 ,所以 s g = Z g T R − 1 ( y g − 1 2 1 ) 、h g = Z g T R − 1 Z g / 4 。直接相加得
U 0 = ( 0 , 3 / 2 ) T , H 0 = ( 2 7 / 4 7 / 4 11 / 4 ) , β 0 ( 1 ) = ( − 14 / 13 , 16 / 13 ) T , U 1 / 4 = ( 0 , 22 / 15 ) T , H 1 / 4 = ( 8 / 5 7 / 5 7 / 5 8 / 3 ) , β 1 / 4 ( 1 ) = ( − 154 / 173 , 176 / 173 ) T . 继续求根而不是把第一步当最终答案,得到
β ^ 0 ≈ ( − 1.23704824 , 1.43812136 ) , β ^ 1 / 4 ≈ ( − 1.04765658 , 1.16809866 ) . 两种工作相关使用同一个总体均值目标,并不要求同一份样本的根相等。对应稳健协方差为
C ^ 0 ≈ ( 1.53049323 − 1.21735680 − 1.21735680 1.34609971 ) , C ^ 1 / 4 ≈ ( 1.32366165 − 1.01834385 − 1.01834385 1.13841632 ) . 在 ρ = 1 / 4 的根上,两个求和矩阵是
H ≈ ( 1.33445209 1.14259095 1.14259095 2.02354509 ) , M ≈ ( 0.73794071 0.57102441 0.57102441 1.68057811 ) . 工作协方差却为
H − 1 ≈ ( 1.45076765 − 0.81917324 − 0.81917324 0.95672686 ) . 它与稳健协方差的对角元有的更大、有的更小;“稳健”没有逐坐标放大的承诺。有限差分求得的观测负 Jacobian 约为 ( 1.28560951 1.18450604 1.13656778 2.01941187 ) ,既不等于 H ,也不对称。ρ = 0 的典范 logit 得分简化为 ∑ Z g T ( Y g − μ g ) ,该特例才恰有观测负 Jacobian 等于 H 。
四簇足以检查这些算术,不能证明95%区间覆盖可靠;也不能用两张样本协方差宣称增大 ρ 总能提高效率。若 R 接近奇异,权重会放大某些残差方向,数值和矩条件都应重新核验。
计数率的迁移只换均值和导数
对已知正暴露量 e g i ,设
log μ g i = log e g i + x g i T β , μ g i = e g i exp ( x g i T β ) . 选 Poisson 型工作方差 A g = diag ( μ g ) ,则 D g = A g X g ,式(1)、(2)以及完整簇条件均值的责任都不变。e β j 是其他输入固定时的边际率比,不是同一个人的随机效应条件率比。真实方差可以不等于均值;只要式(3)和定理其余条件成立,工作方差错误由 B 而非 Q 反映。
使用上一例的四组 x ,将响应改为 ( 0 , 2 ) , ( 1 , 3 ) , ( 2 , 5 ) , ( 2 , 1 ) ,暴露量依次取 ( 1 , 2 ) , ( 2 , 1 ) , ( 1 , 2 ) , ( 2 , 1 ) ,固定 ρ = 1 / 4 。得到 β ^ ≈ ( − 0.47831681 , 0.71402244 ) ,边际率比约为2.04219;稳健协方差约为 ( 0.07520345 − 0.04306183 − 0.04306183 0.02559158 ) 。暴露量不是需要另外估计的回归系数;漏掉 offset 会改变均值模型,而不是仅仅改变标准误。
正权重也能把正确均值拉偏
令 Y ∼ Bernoulli ( 1 / 2 ) ,本来要估计均值 θ ,却偷偷使用响应相关权重 w = 1 + Y 。则
E { w ( Y − 1 / 2 ) } = 1 / 4 ≠ 0 , E { w ( Y − θ ) } = 1 − 3 2 θ . 总体加权根为 2 / 3 ,样本加权根 ∑ ( 1 + Y i ) Y i / ∑ ( 1 + Y i ) 也收敛到 2 / 3 。所有权重都为正,失败仍发生:权重无法从条件期望外提出。围绕这个错误中心计算 sandwich,不会恢复原来的均值 1 / 2 。这与“给定设计后固定的工作相关允许误设”是不同情况。
不能随名称一起带入的保证
从同一批响应估计 R 、分散参数或权重,需要额外论证,例如给出 nuisance 展开、收敛率与导数控制,或先证明某个共同尺度在根与 sandwich 中完全抵消。本页没有覆盖一般估计权重;固定规则的式(3)不能独自承担这项推广。
失访、只保留完整个案或其他依赖响应的选择会改变观测得分的条件均值。不能把“允许相关误设”读成任意缺失机制都安全。
错误链接或遗漏均值结构,可能使区间围绕另一个方程目标;robust 协方差不修复均值目标。边际 logit 系数也不等于条件随机截距模型的系数。
IID 簇是本页独立单位。更多簇内行不能替代更多独立簇;跨簇依赖、高维增长及无界簇大小均需另证。精确根处 ∑ g s g = 0 ,因此式(2)的秩至多 min ( p , G − 1 ) ,与线性回归的少簇秩障碍 理路 聚类稳健协方差估计 Cluster-robust covariance estimation · CR0 covariance 把独立群的总得分作为方差单位,计算回归 CR0 并识别少群、错误分组和秩不足的边界。 相同。
推论与应用
分析重复响应时,先把均值模型与工作相关分别声明,再用独立簇组织求根和不确定性计算。最短学习路径是:会 GLM 的均值导数,理解 Z 的根定位,再用 sandwich 的定标,最后在完整簇上应用这三者。
完成重复响应的独立选读验收 ,可以将两种二元根、计数 offset 与权重偏差放在同一份报告中。标准库复算器 可直接运行,输出每簇得分、拟合均值、H , M 、两种协方差、求根轨迹及观测导数诊断;它验证算术,不自动认证科学假设或有限样本覆盖。
自测与答案
为什么把 R 固定后还要更新方差?答案:固定的是相关规则;Bernoulli 方差 μ ( 1 − μ ) 和计数工作方差 μ 随 β 变化。冻结全部 V 是另一条迭代描述。
只知道 E s g ( β 0 ) = 0 ,可直接说算法一致吗?答案:不可以,还需总体分离、统一逼近与正确根选择;均值零只把真参数放进候选根集合。
两个参数、两个簇,H 可逆,能反演完整稳健协方差吗?答案:精确根处其秩至多1,不能据此构造通常二维 Wald 椭球;数值残差带来的微小第二特征值不是真实信息。
参考资料