Использует ли emmeans lm_robust кластерные стандартные ошибки при вычислении доверительных интервалов для преобразованных переменных результата?

Question

Я исследую взаимодействия между двумя непрерывными переменными-предикторами с помощью пакета emmeans. Я использую lm_robust() из пакета estimatr для выполнения линейной регрессии и получения устойчивых к кластеру стандартных ошибок. Переменная результата центрируется и масштабируется до отклонения единицы SD. Например:

fit <- lm_robust(scale(Y) ~ X1 * X2 + X3 + X4, data = mydata, cluster = school, se_type = 'CR2')

Затем я могу выполнить попарные контрасты или визуализировать линии на трех уровнях X2, используя код, похожий на:

emmip(fit, X2 ~ X1, CIs = TRUE, at = list(X2 = c(mean(X2) - sd(X2),
                                                 mean(X2),
                                                 mean(X2) + sd(X2))))

I не используйте sh, чтобы преобразовать конечную переменную в исходную шкалу.

Мой вопрос заключается в том, использует ли emmeans стандартные ошибки кластерной устойчивости для вычисления доверительных интервалов или значений p, которые он сообщает и зависит ли это поведение от того, находится ли переменная результата в исходном масштабе или преобразована? В коротком примере на сайте создателей пакетов estimatr предполагается, что объекты lm_robust могут использоваться с emmeans, но я не вижу lm_robust в списке поддерживаемых моделей в разделе «Модели. поддерживается виньеткой emmeans page или документацией пакета.

Answer 1

Я считаю, что объекты lm_robust являются расширением lm, поэтому он использует поддержку emmeans для lm. В свою очередь, это будет означать, что оценки получены с помощью coef (модели), а их SE получены с помощью vcov (модели). Поэтому, если vcov () возвращает нужные устойчивые дисперсии, emmeans будет их использовать.

Что касается большинства трансформаций, это будет работать, как описано в виньетке о трансформациях. В частности, указание type = "response" приводит к обратному преобразованию оценок и доверительных интервалов, значения P оставить в покое, а SE должны вычисляться дельта-методом (но не используется в КИ и тестах).

Дополнительная информация

Во-первых, я выяснил, что lm_robust не наследуется от lm; скорее, пакет Estimatr включает собственную поддержку emmeans . Приведено не так много подробностей, но разработчик Estimatr должен полагать, что предоставленное должно быть подходящим.

Преобразование scale() не является встроенным, потому что оно сложно. Сказать, что мы использовали "scale", не так просто, как сказать, что это "log", скажем, потому что для работы с scale() результатами нам нужно знать, что было использовано для центрирования и разделения результатов.

Обходной путь - создать объект, который emmeans() и его родственники должны инвертировать преобразование; и это список функций формы, возвращаемых stats::make.link() или emmeans::make.tran(). Вот функция, которая будет служить этой цели:

make.scaletran = function(y, ...) {
    sy = scale(y, ...)
    if(is.null(m <- attr(sy, "scaled:center")))
       m = 0
    if(is.null(s <- attr(sy, "scaled:scale")))
        s = 1
    list(
        linkfun = function(mu) (mu - m) / s,
        linkinv = function(eta) s * eta + m,
        mu.eta = function(eta) s,
        valideta = function(eta) TRUE,
        name = paste0("scale(", signif(m, 3), ", ", signif(s, 3), ")")
    )
}

Чтобы использовать это, вам необходимо вручную указать преобразование, поскольку оно не определяется автоматически. Вот пример с использованием данных warpbreaks, уже доступных в R:

> warp.lmr = lm_robust(scale(breaks) ~ tension, cluster = wool, 
+     se_type = 'CR2', data = warpbreaks)

> tran = make.scaletran(warpbreaks$breaks)

> emmeans(warp.lmr, "tension", tran = tran)
 tension emmean    SE df lower.CL upper.CL
 L        0.624 0.619 51   -0.618   1.8666
 M       -0.133 0.181 51   -0.497   0.2301
 H       -0.491 0.219 51   -0.930  -0.0517

Results are given on the scale(28.1, 13.2) (not the response) scale. 
Confidence level used: 0.95 

> emmeans(warp.lmr, "tension", tran = tran, type = "response")
 tension response   SE df lower.CL upper.CL
 L           36.4 8.17 51     20.0     52.8
 M           26.4 2.39 51     21.6     31.2
 H           21.7 2.89 51     15.9     27.5

Confidence level used: 0.95 
Intervals are back-transformed from the scale(28.1, 13.2) scale 

Код в OP для вызова emmip() неверен, поскольку он использует спецификации для emmeans(), а не emmip().

Я рассмотрю возможность добавления этой опции преобразования масштаба в emmeans::make.tran() в будущем обновлении.