1

请考虑以下几点：

当用 GEE 拟合时，geepack我们会收到一个模型，我们可以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 包（v0.2.1）于 2019 年 2 月 8 日创建

4

1 回答 1

3

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

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

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

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

有一种情况适合使用不同的百分比值，这就是百分比是“中介变量”的情况；也就是说，百分比落在治疗和反应之间的因果关系中，因此group被认为影响percent以及反应。请参阅关于中介协变量的小节中关于混乱数据的小插图。

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

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

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

于 2019-02-08T13:47:51.043 回答