计算边际效应:ggeffect与ggemmeans结果为何存在差异?
示例
library(glmmTMB) library(ggeffects) ## 零膨胀负二项混合模型 (m <- glmmTMB(count ~ spp + mined + (1|site), ziformula=~spp + mined, family=nbinom2, data=Salamanders, na.action = "na.fail")) summary(m) ggemmeans(m, terms="spp") spp | Predicted | 95% CI -------------------------------- GP | 1.11 | [0.66, 1.86] PR | 0.42 | [0.11, 1.59] DM | 1.32 | [0.81, 2.13] EC-A | 0.75 | [0.37, 1.53] EC-L | 1.81 | [1.09, 3.00] DES-L | 2.00 | [1.25, 3.21] DF | 0.99 | [0.61, 1.62] ggeffects::ggeffect(m, terms="spp") spp | Predicted | 95% CI -------------------------------- GP | 1.14 | [0.69, 1.90] PR | 0.44 | [0.12, 1.63] DM | 1.36 | [0.85, 2.18] EC-A | 0.78 | [0.39, 1.57] EC-L | 1.87 | [1.13, 3.07] DES-L | 2.06 | [1.30, 3.28] DF | 1.02 | [0.63, 1.65]
问题
为何ggeffect与ggemmeans计算输出的边际效应结果存在差异?该差异是否仅源于emmeans、effects两个包的内部计算逻辑差异?此外,是否有相关参考资料可指导从零计算本例这类零膨胀负二项混合模型的边际效应?
回答
两个函数输出结果有差异,核心是非焦点预测变量的边际化规则不一致,绝非只是emmeans和effects两个包的底层计算逻辑差异这么简单,具体区别可以拆成三点:
- 非焦点分类变量的权重规则不同:
ggemmeans()调用emmeans计算时,对模型中的非焦点分类协变量(本例是mined变量)默认给每个水平分配相等权重,再计算加权平均预测值;ggeffect()调用effects计算时,默认按照非焦点分类变量各水平在原始样本中的实际占比分配权重,再做加权平均。这是本例结果出现小幅偏差的最主要原因。 - 随机效应的处理逻辑不同:
ggemmeans()默认把随机效应项设为0,直接输出总体层面的平均预测值,不考虑不同分组的随机截距差异;ggeffect()默认会将随机效应纳入预测过程,按观测层面的随机效应分布做边际化,当随机效应方差较大时,这部分带来的结果差异会非常明显。 - 零膨胀组件的实现细节不同:针对零膨胀类模型,两个包在计数部分、零膨胀部分的协变量积分顺序、链接函数反变换的平均步骤上存在细微实现差异,会进一步拉大结果的偏差。
如果要从零计算这类零膨胀负二项混合模型的边际效应,可以按以下逻辑推进,相关参考内容都可以在所用R包的内置文档中找到,不需要跳转外部站点:
- 先明确定义你要计算的边际效应类型:是平均边际效应(AME,即对所有样本的协变量取值逐一做预测,再对所有预测值取平均),还是均值处边际效应(MEM,即将所有非焦点协变量固定在均值/参考水平,再做预测),两种定义对应的计算逻辑和结果差异很大。
- 牢记零膨胀负二项模型的响应尺度期望值公式:
E[Y] = (1 - p) * mu,其中p是零膨胀组件预测的“结构零”概率,mu是负二项计数组件预测的计数期望值,两个值都需要先通过模型估计的系数计算线性预测值,再通过对应链接函数(零膨胀部分一般是logit链接,计数部分是log链接)反变换到响应尺度。 - 按你选定的权重规则,对所有非焦点变量做加权平均/积分:如果是混合模型要算总体水平的边际效应,还需要对随机效应的正态分布做数值积分,不能直接把随机效应固定为0就结束计算。
- 置信区间计算推荐用参数Bootstrap法或者Delta方法,注意一定要先把预测值反变换到响应尺度再计算区间,不能在线性预测尺度算完区间再反变换,否则会带来区间偏差。
相关计算的完整逻辑可以直接查阅glmmTMB包内置的模型预测说明文档、ggeffects包手册中不同函数边际化规则的对比章节,以及广义线性混合模型相关教材中关于边际效应计算的部分。
内容的提问来源于stack exchange,提问作者Constantin
相关产品推荐
相关产品推荐

