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

Question

请考虑以下几点：

当用 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 进行预测也 returns 每组的其他估计均值，但无法估计基本 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-02-08 创建

Answer 1

您认为是计算问题，结果证明是统计问题...

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

在一个极端的例子中，想象这样一种情况，小组没有任何区别，但 percent 却有区别。然后，如果组间 percent 的均值差异足够大，那么我们可能会有统计上不同的均值，但它们会因为 percent 的影响而不同，而不是因为 group 的影响.

出于这个原因，"fair" 比较是通过预测相同百分比的均值获得的 - 例如，数据集中的总体平均百分比。这是 emmeans 中使用的默认方法，结果称为 adjusted means（在设计教科书中查找）。

在这种情况下，使用不同的百分比值是合适的，在这种情况下，百分比是 "mediating variable;"，也就是说，百分比落在治疗和反应之间的因果路径中，因此 group 被认为会影响 percent 以及响应。见 vignette on messy data, in the subsection on mediating covariates.

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

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

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

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

Marginal means and confidence levels per group with emmeans and geepack in R

r

emmeans