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

Я исследую взаимодействие между двумя непрерывными переменными-предикторами с помощью пакета 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))))

Я не хочу возвращать исходную переменную к исходному масштабу.

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


r emmeans
person Paul    schedule 29.05.2020    source источник

Ответы (1)


arrow_upward
3
arrow_downward

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

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

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

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

Преобразование 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() в будущем обновлении.

person Russ Lenth    schedule 29.05.2020
comment
Спасибо за ваш ответ. Что касается второй части, я не думаю, что scale() является признанной трансформацией: Unknown transformation "scale": no transformation done. Я пробовал ваше предложение, и те, что указаны в виньетке, и значения SE, CI, p всегда given on the scale (not the response) scale. Для сравнения, если я использую log() вместо scale(), все будет работать согласно виньетке. - person Paul; 30.05.2020
comment
Ты прав. Я прочитал слово «преобразование», но не посмотрел на то, что вы использовали. Добавлю к ответу. - person Russ Lenth; 30.05.2020
comment
Вы очень любезно выписали эту функцию. Отлично работает, спасибо. Я отредактирую свой ответ, чтобы исправить ошибку кодировки в вызове emmip(). Спасибо еще раз. - person Paul; 30.05.2020