R中使用emmeans和geepack的每组边际均值和置信水平

Posted 2023-02-18

【中文标题】R中使用emmeans和geepack的每组边际均值和置信水平

【英文标题】：Marginal means and confidence levels per group with emmeans and geepack in R

【发布时间】：2019-07-02 13:40:41

【问题描述】：

请考虑以下几点：

当使用geepack 拟合 GEE 时，我们会收到一个模型，我们可以使用新值 predict，但基础 R 不支持 GEE 模型来计算置信区间。要获得置信区间，我们可以使用emmeans::emmeans()。

如果模型中的变量是分类的和连续的，我就会遇到问题。

在使用emmeans::emmeans() 估计边际均值时，我发现边际均值是根据整体数据而不是每组数据计算的。

问题：如何从 R 中的 GEE 模型获得每组的估计平均值，包括置信区间？

---

最小的可重现示例：

数据

library("dplyr")
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library("emmeans")
#> Warning: package 'emmeans' was built under R version 3.5.2
library("geepack")

# Adding a grouping variable
pigs.group <- emmeans::pigs %>% mutate(group = c(rep("a", 20), rep("b", 9)))

拟合模型

# Fitting the model
fit <- geepack::geeglm(conc ~ as.numeric(percent) + factor(group),
                       id = source, data = pigs.group)

# Model results
fit
#> 
#> Call:
#> geepack::geeglm(formula = conc ~ as.numeric(percent) + factor(group), 
#>     data = pigs.group, id = source)
#> 
#> Coefficients:
#>         (Intercept) as.numeric(percent)      factor(group)b 
#>           20.498948            1.049322           10.703857 
#> 
#> Degrees of Freedom: 29 Total (i.e. Null);  26 Residual
#> 
#> Scale Link:                   identity
#> Estimated Scale Parameters:  [1] 36.67949
#> 
#> Correlation:  Structure = independence  
#> Number of clusters:   3   Maximum cluster size: 10

使用emmeans::emmeans() 计算边际均值和 LCL/UCL。但是，percent 的组均值在两个组中均为 12.9。这是percent 的整体观察平均值，而不是组平均值。

# Calculating marginal means per group.
# Note that 'percent' is the same for both groups
emmeans::emmeans(fit, "percent", by = "group")
#> group = a:
#>  percent emmean    SE  df asymp.LCL asymp.UCL
#>     12.9   34.1 3.252 Inf      27.7      40.4
#> 
#> group = b:
#>  percent emmean    SE  df asymp.LCL asymp.UCL
#>     12.9   44.8 0.327 Inf      44.1      45.4
#> 
#> Covariance estimate used: vbeta 
#> Confidence level used: 0.95

# Creating new data with acutal means per group
new.dat <- pigs.group %>%
        group_by(group) %>%
        summarise(percent = mean(percent))

# These are the actual group means
new.dat
#> # A tibble: 2 x 2
#>   group percent
#>   <chr>   <dbl>
#> 1 a        13.2
#> 2 b        12.3

使用predict 进行预测还会返回每组的其他估计均值，但无法估计基数 R 中的 GEE 的置信区间。

# Prediction with new data
# These should be the marginal means but how to get the confidence interval?
predict(fit, newdata = new.dat)
#>        1        2 
#> 34.35000 44.14444

^{由reprex package (v0.2.1) 于 2019 年 2 月 8 日创建}

【参考方案1】：

您认为是计算问题的结果是统计问题...

当模型中有协变量时，事后分析中的常用方法是控制 for 这些协变量。在给定示例的上下文中，我们想要比较不同组中的平均响应。但是，响应也受协变量 percent 的影响，并且每个组的平均百分比不同。如果我们只计算每个组的边际均值，这些均值的不同部分是因为percent 的影响。

在一个极端的例子中，想象一下这样一种情况，该组没有任何区别，但percent 确实如此。那么如果percent 的平均值在各组之间有足够的差异，那么我们可以有统计上不同的平均值，但它们会因为percent 的影响而不同，而不是因为group 的影响。

因此，通过预测相同百分比的均值（例如，数据集中的总体平均百分比）来获得“公平”比较。这是emmeans中使用的默认方法，其结果称为adjusted mean（在设计教科书中查找）。

有一种情况适合使用不同的百分比值，即百分比是“中介变量”的情况；也就是说，百分比落在治疗和反应之间的因果关系中，因此group 被认为会影响percent 以及反应。请参阅vignette on messy data, in the subsection on mediating covariates。

如果你真的认为percent 是一个中介协变量，那么你可以像这样得到单独的百分比：

 emmeans(model, "group", cov.reduce = percent ~ group)

但是，在percent 被视为独立于group 的情况下，不要这样做！

【讨论】：

我不希望每个协变量的边际均值。我想要每个组的平均值的边际平均值（这里：group）。我用你的方法和ref_grid()：summary(ref_grid(fit, at = list(percent = new.dat$percent)))尝试了以下方法。这返回了一个带有预测的边际平均值和 SE 的表。但是没有 LCL/UCL。

试试emmeans(fit, “group”)。

emmeans(fit, "group") 返回与emmeans::emmeans(fit, "percent", by = "group") 相同的边际均值（请参阅上面示例中后者的结果）。这正是促使我写这篇文章的原因。我原以为emmeans(fit, "group") 会使用观察到的（百分比）平均值每组，但由于我收到与emmeans::emmeans(fit, "percent", by = "group") 相同的结果，我得出结论，无论组如何，结果都是基于总体平均值.