如何为geepack::geeglm计算预测标准误？兼论与glmgee的差异

问题背景 #

我想用geepack::geeglm()拟合的GEE模型结果做可视化，不需要针对新个体的预测，但调用predict()并设置se.fit = TRUE时出现错误：Error in scale^2 : non-numeric argument to binary operator。

错误重现示例 #

注：使用自回归相关结构，符合我的实际数据需求

library(geepack)

rm(list = ls())

data(mtcars)
mtcars$am <- as.character(mtcars$am)
mtcars <- mtcars[order(mtcars$cyl, mtcars$disp),]
m1 <- geeglm(mpg ~ disp * am, data = mtcars, id = cyl, corstr = "ar1")
summary(m1)

newd <- expand.grid(disp = seq(min(mtcars$disp), max(mtcars$disp), 1), am = c("0", "1"))

newd$pred <- predict(m1, newdata = newd, type = "response", se.fit = TRUE)

# Error in scale^2 : non-numeric argument to binary operator

glmtoolbox替代实现示例 #

我发现glmtoolbox包的glmgee()函数可以计算预测标准误并生成所需可视化：

library(glmtoolbox)

rm(list = ls())

data(mtcars)
mtcars$am <- as.character(mtcars$am)
mod1 <- mpg ~ disp * am
mtcars <- mtcars[order(mtcars$cyl, mtcars$disp),]
m1 <- glmgee(mod1, id=cyl, family=gaussian(), data=mtcars, corstr="AR-M-dependent")
summary(m1)
newd <- expand.grid(disp = seq(min(mtcars$disp), max(mtcars$disp), 1), am = c("0", "1"))


newd$fit    <- as.data.frame(predict(m1, newdata = newd, type = "response", se.fit = TRUE))$fit
newd$se.fit <- as.data.frame(predict(m1, newdata = newd, type = "response", se.fit = TRUE))$se.fit


ggplot(newd, aes(x = disp, y = fit, color = am)) +
  geom_line() +
  geom_ribbon(aes(ymin = fit - 1.96 * se.fit,
                  ymax = fit + 1.96 * se.fit, fill = am), alpha = 0.2) +
  ylab("fit (mpg)") +
  scale_y_continuous(limits = c(-5,35)) +
  theme_light()

模型结果差异对比 #

但两个函数的模型结果存在明显差异：

• geepack::geeglm()结果：

geeglm(formula = mpg ~ disp * am, data = mtcars, id = cyl, corstr = "ar1")

 Coefficients:
            Estimate  Std.err   Wald Pr(>|W|)    
(Intercept) 26.70389  1.98299  181.3  < 2e-16 ***
disp        -0.03235  0.00571   32.1  1.5e-08 ***
am1          5.04045  0.13026 1497.4  < 2e-16 ***
disp:am1    -0.01908  0.00374   26.0  3.4e-07 ***

• glmtoolbox::glmgee()结果：

Coefficients
            Estimate Std.Error z-value Pr(>|z|)
(Intercept)   26.026     1.691  15.391  < 2e-16
disp          -0.030     0.005  -6.333 2.40e-10
am1            6.207     0.296  20.982  < 2e-16
disp:am1      -0.024     0.006  -4.191 2.78e-05

我已验证，当使用独立相关结构时两者结果一致，推测差异源于相关结构处理或数据排序方式。另外，我尝试用混合模型实现类似功能，但置信区间远宽于预期。

咨询问题 #

1. glmtoolbox::glmgee()与geepack::geeglm()结果差异的原因是什么？
2. 能否为geepack::geeglm()的预测计算标准误？我长期使用该函数，更熟悉其操作。

内容的提问来源于stack exchange，提问作者Avl