Как получить ковариационную матрицу для случайных эффектов (BLUP/условные режимы) из lme4

Итак, я подогнал линейную смешанную модель с двумя случайными перехватами в R:

Y = X beta  + Z b + e_i,

где b ~ MVN (0, Sigma); X и Z — матрицы модели фиксированных и случайных эффектов соответственно, а beta и b — параметры фиксированных эффектов и случайных эффектов BLUP/условные режимы.

Я хотел бы получить в свои руки базовую ковариационную матрицу b, которая не кажется тривиальной вещью в пакете lme4. Вы можете получить только отклонения по VarCorr, а не фактическую матрицу корреляции.

Согласно одному из виньеток пакета (стр. 2) , можно рассчитать ковариацию бета: e_i * lambda * t(lambda). И все эти компоненты вы можете извлечь из вывода lme4.

Мне было интересно, если это путь? Или у вас есть другие предложения?


r matrix lme4 covariance
person pa_ka    schedule 13.06.2016    source источник
comment
Добро пожаловать в StackOverflow. Пожалуйста, уточните свои математические обозначения (Xbeta на самом деле X * beta? Возможно, вам также следует указать, что такое beta, b, sigma; хотя я не эксперт, и для некоторых эти обозначения могут быть очевидными). Помните, что чем четче вы напишете свой вопрос, тем вероятнее и быстрее получите соответствующий ответ.   -  person YakovL    schedule 13.06.2016
comment
Да, Xbeta должна была быть X*beta. Бета — вектор фиксированных эффектов матрицы плана X, b — вектор случайных эффектов, а sigma — ковариационная матрица дисперсии b. Спасибо за подсказку, буду иметь в виду.   -  person pa_ka    schedule 13.06.2016

Ответы (1)


arrow_upward
5
arrow_downward

От 1_:

Если «condVar» имеет значение «TRUE», каждый из фреймов данных имеет атрибут, называемый «postVar», который представляет собой трехмерный массив с симметричными гранями; каждая грань содержит матрицу дисперсии-ковариации для определенного уровня фактора группировки. (Имя этого атрибута является историческим артефактом и может быть изменено на «condVar» в какой-то момент в будущем.)

Установите пример:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
rr <- ranef(fm1,condVar=TRUE)

Получите матрицу дисперсии-ковариации среди значений b для перехвата

pv <- attr(rr[[1]],"postVar")
str(pv)
##num [1:2, 1:2, 1:18] 145.71 -21.44 -21.44 5.31 145.71 ...

Итак, это массив 2x2x18; каждый срез представляет собой матрицу дисперсии-ковариации между условным пересечением и наклоном для конкретного субъекта (по определению, пересечения и наклоны для каждого субъекта не зависят от пересечений и наклонов для всех других субъектов).

Чтобы преобразовать это в матрицу дисперсии-ковариации (см. getMethod("image",sig="dgTMatrix") ...)

library(Matrix)
vc <- bdiag(  ## make a block-diagonal matrix
            lapply(
                ## split 3d array into a list of sub-matrices
                split(pv,slice.index(pv,3)),
                   ## ... put them back into 2x2 matrices
                      matrix,2)) 
image(vc,sub="",xlab="",ylab="",useRaster=TRUE)

введите описание изображения здесь

person Ben Bolker    schedule 13.06.2016
comment
Большое спасибо! Я проверю это :) - person pa_ka; 13.06.2016
comment
StackOverflow не рекомендует использовать комментарии для выражения благодарности; если этот ответ был полезен, вы можете проголосовать за него (если у вас достаточно репутации), и в любом случае, если он удовлетворительно отвечает на ваш вопрос, вам предлагается нажать на галочку, чтобы принять его. - person Ben Bolker; 13.06.2016
comment
Извините, это может быть некоторое время. Но я пытался понять ваше слово each slice is the variance-covariance matrix among the conditional intercept and slope for a particular subject. Так что же такое каждый срез? Вы имели в виду ковариацию b, как спрашивал автор? - person Nan; 19.02.2019
comment
Причина, по которой я спросил, заключается в том, что я пытался использовать ваш код, чтобы получить оценку матрицы ковариации b, но результат сильно отличается от фактической ковариации b, которую мы использовали для генерации данных. Большое спасибо! - person Nan; 19.02.2019
comment
И в чем разница между получением ковариационной матрицы с использованием вашего кода выше и чтением из Std.Dev. в выводе lmer? Согласно заголовку stats.stackexchange.com/questions/24452/, они должны означать одно и то же, но в моих экспериментах они не совпадают. . - person Nan; 19.02.2019
comment
Я публикую вопрос по этому вопросу, вот ссылка stackoverflow.com/questions/54770082/ Спасибо за любую помощь! - person Nan; 19.02.2019