Дисперсия неизвестного контраста от lm или coxph

Question

Со вчерашнего дня я пользуюсь своим вопросом

Предположим, у нас есть три процедуры, и мы хотим получить все парные различия. По умолчанию в R используется контрастность и отображается только 2. IE 2 против 1 и 3 против 1. Чтобы получить 3 против 2, нам нужно вычесть бета-коэффициенты (3 против 1 - 2 против 1).

Итак, у нас есть оценка сейчас, но есть ли способ получить дисперсию этой оценки? Или мы должны снова запустить регрессию с другой контрольной группой, чтобы получить ее?

Answer 1

Мы можем использовать glht функция от multcomp пакет для проведения пост-специального тестирования с возможностью настройки для множественного тестирования. (см. http://www.ats.ucla.edu/stat/r/faq/testing_contrasts.htm для большего количества примеров)

Маленький пример

library(multcomp)
data(mtcars)
mtcars$cyl <- factor(mtcars$cyl, levels=c(4, 6, 8))

# run model
m <- lm(mpg ~ cyl, data=mtcars)
coef(summary(m))
#               Estimate Std. Error   t value     Pr(>|t|)
# (Intercept)  26.663636  0.9718008 27.437347 2.688358e-22
# cyl6         -6.920779  1.5583482 -4.441099 1.194696e-04
# cyl8        -11.563636  1.2986235 -8.904534 8.568209e-10
coef(m)[2] - coef(m)[3]
#     cyl6 
# 4.642857 

# use glht function to estimate other contrasts
# define the linfct matrix to test specific contrast
mf <- glht(m, linfct = matrix(c(0, 1, -1), 1))
summary(mf)
#    Simultaneous Tests for General Linear Hypotheses
# 
# Fit: lm(formula = mpg ~ cyl, data = mtcars)
# 
# Linear Hypotheses:
#        Estimate Std. Error t value Pr(>|t|)   
# 1 == 0    4.643      1.492   3.112  0.00415 **
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)

# look at estimate by changing the reference level manually
mtcars$cyl2 <- factor(mtcars$cyl, levels=c(6, 4, 8))
m2 <- lm(mpg ~ cyl2, data=mtcars)
coef(summary(m2))
#              Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 19.742857   1.218217 16.206357 4.493289e-16
# cyl24        6.920779   1.558348  4.441099 1.194696e-04
# cyl28       -4.642857   1.492005 -3.111825 4.152209e-03

Мы также можем получить ответ, зная теорию, стоящую за ним. Возможно, этот ответ лучше подходит для резюме, но я решил опубликовать его здесь.

Ответ сводится к знанию статистики, стоящей за этим. У нас есть разница двух случайных величин, которые мы хотим получить дисперсию.

В математическом плане мы хотели бы var(beta_{2vs1} -beta_{3vs1}), Зная нашу элементарную статистику, мы знаем, что это равно var(beta_{2vs1}) +var(beta_{3vs1})-2*cov(beta_{2vs1},beta_{3vs1})

v <- vcov(m)
sqrt(v[2,2] + v[3,3] - 2*v[2,3])
# [1] 1.492005