Предельные средние и уровни достоверности для группы с emmeans и geepack в R

Question

Обратите внимание на следующее:

При подборе GEE с geepack мы получаем модель, в которой мы можем predict с новыми значениями, но база R не поддерживает модели GEE для вычисления доверительных интервалов.Чтобы получить доверительные интервалы, мы можем использовать emmeans::emmeans().

Если переменные в модели являются категоричными и непрерывными, у меня возникают проблемы.

При оценке предельного среднего с помощью emmeans::emmeans() я обнаружил, чтопредельное среднее значение рассчитывается на основе общих данных, а не данных по группе.

Вопрос: как получить расчетное среднее значение по группе, включая доверительные интервалы, из модели GEE в R?

---

Минимальный воспроизводимый пример:

Данные

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 также возвращает другие оцененные средние значения для группы, но доверительные интервалы для GEE не могут быть оценены для базы R.

# 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

^{Создано в 2019-02-08 пакетом Представить (v0.2.1)}

Answer 1

То, что вы считали вычислительной проблемой, оказывается статистической ...

Когда в модели есть ковариаты, обычный подход в пост-hoc анализе заключается в контроле для эти ковариаты.В контексте данного примера мы хотим сравнить средний ответ в разных группах.Однако на ответ также влияет ковариата, percent, и средний процент отличается в каждой группе.Если мы просто вычисляем предельные средние значения для каждой группы, эти средства частично отличаются из-за эффектов percent.

В крайнем примере, представьте ситуацию, в которой группа не имеет никакого значения, а percent.Тогда, если средние значения percent достаточно различаются между группами, тогда мы могли бы иметь статистически разные средние значения, но они были бы разными из-за эффекта percent, а не из-за эффекта group.

По этой причине "справедливое" сравнение получается путем прогнозирования средних значений на уровне таких же процентов, скажем, общего среднего процента в наборе данных.Это метод по умолчанию, используемый в emmeans , и результаты называются скорректированные средства (найдите его в учебнике по дизайну).

В некоторых ситуацияхуместно использовать разные значения процента, и в этом случае процент является «посреднической переменной»;то есть процент попадает в причинно-следственную связь между лечением и ответом, так что, как полагают, group влияет как на 1024 *, так и на ответ.См. виньетка с беспорядочными данными в подразделе о посреднических ковариатах .

Если вы действительно думаете, что percent является посреднической ковариатой, то вы можете получить отдельные проценты, например:

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

Однако, в ситуации, когда percent считается независимым от group, не делайте этого!