Первые линейные контрасты: Кенвард Роджерс или Саттертуэйт DF и SE

В R я ищу способ оценить доверительные интервалы для линейных контрастов для lmer моделей, которые используют либо степени свободы Кенварда-Роджера, либо Саттертуэйта и SE.

Например, я могу вычислить CI для параметра с фиксированным эффектом в смешанной модели, такой как SAS с R, используя t-значение (с df из KR) и SE.

mod<-lmerTest::lmer(y~time1+treatment+time1:treatment+(1|PersonID),data=data)
lmerTest::summary(mod,ddf = "Kenward-Roger")

Этот вывод:

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      49.0768     1.0435 56.4700  47.029  < 2e-16 ***
time1             5.8224     0.5963 48.0000   9.764 5.51e-13 ***
treatment         1.6819     1.4758 56.4700   1.140   0.2592    
time1:treatment   2.0425     0.8433 48.0000   2.422   0.0193 * 

Разрешает ЭК для time1, например:

5.8224+abs(qt(0.05/2, 48))*0.5963 #7.021342
5.8224-abs(qt(0.05/2, 48))*0.5963 #4.623458

Я хотел бы сделать то же самое для линейного контраста фиксированных коэффициентов. Это p-значение, но выход SE отсутствует.

pbkrtest::KRmodcomp(mod,matrix(c(0,0,1,0),nrow = 1)) 

         stat     ndf     ddf F.scaling p.value
Ftest  1.2989  1.0000 56.4670         1  0.2592

Можно ли каким-либо образом получить SE или CI из более линейных контрастов, использующих этот тип df?


person B_Miner    schedule 21.08.2015    source источник


Ответы (1)


Для этого у вас есть как минимум два варианта: использовать пакет lsmeans или сделать это вручную (используя функции vcovAdj.lmerMod и pbkrtest::get_Lb_ddf). Лично я выбираю более поздний вариант, если тестируемый контраст не очень «простой», потому что я нахожу синтаксис в lsmeans немного сложным.

Для примера возьмем следующую модель:

library(pbkrtest)
library(lme4)
library(nlme) # for the 'Orthodont' data

# 'age' is a numeric variable, while 'Sex' and 'Subject' are factors
model <- lmer(distance ~ age : Sex + (1 | Subject), data = Orthodont)

Linear mixed model fit by REML ['lmerMod']
Formula: distance ~ age:Sex + (1 | Subject)
…
Fixed Effects:
    (Intercept)    age:SexMale  age:SexFemale  
        16.7611         0.7555         0.5215 

из которого мы хотели бы получить статистику о разнице между коэффициентами возраста у мужчин и женщин (т. е. age:SexMale - age:SexFemale).

Использование lsmeans:

library(lsmeans)
# Evaluate the contrast at a value of 'age' set to 1,
# so that the resulting value is equal to the regression coefficient
lsm = lsmeans(model, pairwise ~ age : Sex, at = list(age = 1))$contrast

производит:

contrast           estimate         SE    df t.ratio p.value
1,Male - 1,Female 0.2340135 0.06113276 42.64   3.828  0.0004

В качестве альтернативы, выполняя расчет вручную:

# Specify the contrasts: age:SexMale - age:SexFemale
# Must have the same order as the fixed effects in the model
K = c("(Intercept)" = 0, "age:SexMale" = 1, "age:SexFemale" = -1)

# Retrieve the adjusted variance-covariance matrix, to calculate the SE
V = pbkrtest::vcovAdj.lmerMod(model, 0)

# Point estimate, SE and df
point_est = sum(K * fixef(model))
SE = sqrt(sum(K * (V %*% K)))
df = pbkrtest::get_Lb_ddf(model, K)

alpha = 0.05 # significance level

# Calculate confidence interval for the difference between the 'age' coefficients for males and females
Delta_age_CI = point_est + SE * qt(c(0.5 * alpha, 1 - 0.5 * alpha), df)

приведет к точечной оценке, равной 0.2340135, SE 0.06113276, df 42.63844 и доверительному интервалу [0.1106973, 0.3573297]

person tguzella    schedule 21.08.2015
comment
Какой фантастический ответ! Спасибо большое - очень понятно! Я получил известие от одного из авторов lmerTest, что они рассматривают возможность добавления возможности оценки контрастов в свой следующий выпуск. - person B_Miner; 21.08.2015
comment
Интересно, что вы находите семь операторов для ручного расчета проще, чем с lsmeans. На самом деле, если вы используете возраст как переменную by, все еще проще: lsmeans(model, pairwise ~ Sex | age) - person Russ Lenth; 22.08.2015