R geepack：使用GEE得到过大的估计

Question

我正在使用geepack进行R语言中的逻辑边际模型估计，使用geeglm()函数。但是我得到了垃圾估计值，它们大约比期望值高16个数量级。然而，p值似乎与预期相似。这意味着响应基本上变成了一个阶梯函数。请参见附图 以下是生成该图的代码：

require(geepack)
data = read.csv(url("http://folk.uio.no/mariujon/data.csv"))
fit = geeglm(moden ~ 1 + power, id = defacto, data=data, corstr = "exchangeable", family=binomial)
summary(fit)
plot(moden ~ power, data=data)
x = 0:2500
y = predict(fit, newdata=data.frame(power = x), type="response" )
lines(x,y)

这是回归表格：

Call:
geeglm(formula = moden ~ 1 + power, family = binomial, data = data, 
    id = defacto, corstr = "exchangeable")

 Coefficients:
             Estimate   Std.err  Wald Pr(>|W|)    
(Intercept) -7.38e+15  1.47e+15  25.1  5.4e-07 ***
power        2.05e+13  1.60e+12 164.4  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Estimated Scale Parameters:
            Estimate  Std.err
(Intercept) 1.03e+15 1.65e+37

Correlation: Structure = exchangeable  Link = identity 

Estimated Correlation Parameters:
      Estimate  Std.err
alpha    0.196 3.15e+21
Number of clusters:   3   Maximum cluster size: 381

希望能得到帮助。谢谢！ 祝好， Marius

- Mikkel Rev

---

你需要一些正则化或收缩组件。你可以使用广义线性混合模型+贝叶斯先验来处理固定效应(MCMCglmm或blme软件包)，但它将拟合条件模型而不是边际模型...我不知道如何在GEE框架中实现收缩，或者是否已经有人做过这个。 - Ben Bolker

我有一种边际逻辑方法，(Intercept) 的系数为-0.664，power 的系数为0.003。是否有兴趣写出来？ - swihart

@swihart：当然可以。 - Mikkel Rev

出于好奇，什么是数据应用？我很感兴趣，因为我通常在有很多簇但每个簇只有几个观测值的情况下工作，而这里的一个簇有3个簇和381个观测值。 - swihart

@swihart 有一个生物学应用程序。在一项实验中，数百个个体被养育在完全相同的三种环境中。我们想研究个体在给定身体质量指数的情况下成熟的概率。但我们预计环境会引起相关性。 - Mikkel Rev

Answer 1

我将提供三个过程，每个过程都是边际随机截距模型（MRIM）。这些MRIM具有边际逻辑解释的系数，并且比GEE的系数小：

| Model | (Intercept) |  power |  LogL  |
|-------|-------------|--------|--------| 
| `L_N` |       -1.050| 0.00267|  -270.1|
| `LLB` |       -0.668| 0.00343|  -273.8|
| `LPN` |       -1.178| 0.00569|  -266.4|

相对于不考虑任何相关性的glm，以供参考：

| Model | (Intercept) |  power |  LogL  |
|-------|-------------|--------|--------| 
|  strt |       -0.207| 0.00216|  -317.1|

一个边缘化的随机截距模型（MRIM）值得探索，因为您需要具有可交换相关结构的边缘模型来处理聚类数据，而这是MRIM所展示的结构类型。

代码（尤其是带有注释的R脚本）和文献的PDF文件在GITHUB仓库中。我将在下面详细介绍代码和文献。

MRIM的概念自1999年以来就存在了，在GITHUB仓库中有一些背景阅读材料。我建议首先阅读Swihart et al 2014，因为它回顾了其他论文。

按时间顺序--

• L_N Heagerty (1999): 该方法使用一个正态分布的随机截距逻辑回归模型。关键在于，随机截距模型中的预测器非线性参数化，并带有边际系数，以使得所得到的边际模型具有边际逻辑解释。其代码是lnMLE R软件包（不在CRAN上，但可以在Patrick Heagerty的网站这里找到）。该方法在代码中被标记为L_N，表示边际上的logit（L），在条件比例尺上没有解释（_），且随机截距服从正态分布（N）。

• LLB Wang & Louis (2003): 该方法使用一个桥式分布的随机截距逻辑回归模型。与Heagerty 1999不同的是，这个技巧是一种特殊的随机效应分布（桥式分布），它允许随机截距模型和所得到的边际模型都具有逻辑解释。其代码是使用gnlmix4MMM.R（在repo中），该代码使用rmutil和repeated R软件包实现。该方法在代码中被标记为LLB，表示边际上的logit（L），在条件比例尺上也是logit（L），且随机截距服从桥式分布（B）。

• LPN Caffo and Griswold (2006): 该方法使用一个正态分布的随机截距probit模型，而Heagerty 1999使用了一个logit随机截距模型。这种替换使计算变得更容易，并仍然产生一个边际logit模型。其代码是使用gnlmix4MMM.R（在repo中），该代码使用rmutil和repeated R软件包实现。该方法在代码中被标记为LPN，表示边际上的logit（L），在条件比例尺上是probit（P），且随机截距服从正态分布（N）。

• Griswold et al (2013): 另一篇评论/实用介绍。

• Swihart et al 2014: 这是一篇关于Heagerty 1999和Wang & Louis 2003以及其他人的评论文章，它推广了MRIM方法。最有趣的推广之一是允许边际和条件模型中的logistic CDF（等效地，logit link）改为近似于logistic CDF的稳定分布。其代码是使用gnlmix4MMM.R（在repo中），该代码使用rmutil和repeated R软件包实现。我在带注释的R脚本中用SSS表示这个方法，表示边际上是稳定的（S），在条件比例尺上也是稳定的（S），且随机截距服从稳定分布（S）。它包含在R脚本中，但在SO的这篇文章中没有详细介绍。

准备

#code from OP Question: edit `data` to `d` 
require(geepack)
d = read.csv(url("http://folk.uio.no/mariujon/data.csv"))
fit = geeglm(moden ~ 1 + power, id = defacto, data=d, corstr = "exchangeable", family=binomial)
summary(fit)
plot(moden ~ power, data=d)
x = 0:2500
y = predict(fit, newdata=data.frame(power = x), type="response" )
lines(x,y)
#get some starting values from glm():
strt <- coef(glm(moden ~  power, family = binomial, data=d))
strt
#I'm so sorry but these methods use attach()
attach(d)

L_N Heagerty（1999年）

# marginally specifies a logit link and has a nonlinear conditional model
# the following code will not run if lnMLE is not successfully installed.  
# See https://faculty.washington.edu/heagerty/Software/LDA/MLV/
library(lnMLE)
L_N <- logit.normal.mle(meanmodel = moden ~ power,
                        logSigma= ~1,
                        id=defacto,
                        model="marginal",
                        data=d,
                        beta=strt,
                        r=10) 
print.logit.normal.mle(L_N)

准备 LLB 和 LPN

library("gnlm")
library("repeated")
source("gnlmix4MMM.R") ## see ?gnlmix; in GITHUB repo 
y <- cbind(d$moden,(1-d$moden))

LLB Wang和Louis（2003年）

LLB  <- gnlmix4MMM(y = y,
                   distribution = "binomial",
                   mixture = "normal",
                   random = "rand",
                   nest = defacto,
                   mu = ~ 1/(1+exp(-(a0 + a1*power)*sqrt(1+3/pi/pi*exp(pmix)) - sqrt(1+3/pi/pi*exp(pmix))*log(sin(pi*pnorm(rand/sqrt(exp(pmix)))/sqrt(1+3/pi/pi*exp(pmix)))/sin(pi*(1-pnorm(rand/sqrt(exp(pmix))))/sqrt(1+3/pi/pi*exp(pmix)))))),
                   pmu = c(strt, log(1)),
                   pmix = log(1))

print("code: 1 -best 2-ok 3,4,5 - problem")
LLB$code
print("coefficients")
LLB$coeff
print("se")
LLB$se

LPN Caffo和Griswold（2006年）

LPN  <- gnlmix4MMM(y = y,
                   distribution = "binomial",
                   mixture = "normal",
                   random = "rand",
                   nest = defacto,
                   mu = ~pnorm(qnorm(1/(1+exp(-a0 - a1*power)))*sqrt(1+exp(pmix)) + rand),
                   pmu = c(strt, log(1)),
                   pmix = log(1))

print("code: 1 -best 2-ok 3,4,5 - problem")
LPN$code
print("coefficients")
LPN$coeff
print("se")
LPN$se

三种方法得出的系数：

rbind("L_N"=L_N$beta, "LLB" = LLB$coefficients[1:2], "LPN"=LPN$coefficients[1:2])

3个模型的最大对数似然：

rbind("L_N"=L_N$logL, "LLB" = -LLB$maxlike, "LPN"=-LPN$maxlike)

- swihart