Маржинальные средние и уровни доверия на группу с 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

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

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

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

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

Существует ситуация, когда уместно использовать разные значения процента, и именно в этом случае процент является "посреднической переменной"; то есть проценты попадают в причинно-следственную связь между лечением и реакцией, так что group Считается, что влияет percent а также ответ. Смотрите виньетка с грязными данными в подразделе о посреднических ковариатах.

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

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

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