Модельные матрицы для моделей со смешанными эффектами

10

В lmerфункции lme4in in Rесть вызов для построения модельной матрицы случайных эффектов , как описано здесь , стр. 7 - 9.Z

Вычисление влечет за собой произведения ХатриРао и / или Кронекера двух матриц, и . ZJiXi

Матрица представляет собой глоток: «Матрица индикаторов группирующих индексов», но, похоже, она представляет собой разреженную матрицу с фиктивным кодированием, чтобы выбрать, какая единица (например, субъекты в повторяющихся измерениях), соответствующие более высоким иерархическим уровням, «включена» для любое наблюдение. матрица , кажется, действует как селектор измерений в нижнем иерархическом уровне, так что сочетание обоих «селекторов» даст матрицу, вида , показанного на бумаге с помощью следующего примера:JiXiZi

(f<-gl(3,2))

[1] 1 1 2 2 3 3
Levels: 1 2 3

(Ji<-t(as(f,Class="sparseMatrix")))

6 x 3 sparse Matrix of class "dgCMatrix"
     1 2 3
[1,] 1 . .
[2,] 1 . .
[3,] . 1 .
[4,] . 1 .
[5,] . . 1
[6,] . . 1

(Xi<-cbind(1,rep.int(c(-1,1),3L)))
     [,1] [,2]
[1,]    1   -1
[2,]    1    1
[3,]    1   -1
[4,]    1    1
[5,]    1   -1
[6,]    1    1

Транспонирование каждой из этих матриц и выполнение умножения Хатри-Рао:

[11......11......11][111111111111]=[11....11......11....11......11....11]

Но это транспонирование этого:Zi

(Zi<-t(KhatriRao(t(Ji),t(Xi))))

6 x 6 sparse Matrix of class "dgCMatrix"

[1,] 1 -1 .  . .  .
[2,] 1  1 .  . .  .
[3,] .  . 1 -1 .  .
[4,] .  . 1  1 .  .
[5,] .  . .  . 1 -1
[6,] .  . .  . 1  1

Оказывается, что авторы используют базу данных sleepstudyв lme4, но на самом деле не уточняют матрицы проектирования, поскольку они применимы к этому конкретному исследованию. Поэтому я пытаюсь понять, как составленный код из статьи, приведенной выше, может sleepstudyпривести к более осмысленному примеру.

Для наглядности я сократил набор данных до трех предметов - «309», «330» и «371»:

require(lme4)
sleepstudy <- sleepstudy[sleepstudy$Subject %in% c(309, 330, 371), ]
rownames(sleepstudy) <- NULL

Каждый индивидуум будет демонстрировать очень разные точки пересечения и наклона, если простая регрессия OLS будет рассматриваться индивидуально, что предполагает необходимость модели смешанного эффекта с более высокой иерархией или уровнем единицы, соответствующей субъектам:

    par(bg = 'peachpuff')
    plot(1,type="n", xlim=c(0, 12), ylim=c(200, 360),
             xlab='Days', ylab='Reaction')
    for (i in sleepstudy$Subject){
            fit<-lm(Reaction ~ Days, sleepstudy[sleepstudy$Subject==i,])
            lines(predict(fit), col=i, lwd=3)
            text(x=11, y=predict(fit, data.frame(Days=9)), cex=0.6,labels=i)
        }

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

Вызов регрессии со смешанным эффектом:

fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)

И матрица, извлеченная из функции, дает следующее:

parsedFormula<-lFormula(formula= Reaction~Days+(Days|Subject),data= sleepstudy)
parsedFormula$reTrms

$Ztlist
$Ztlist$`Days | Subject`
6 x 12 sparse Matrix of class "dgCMatrix"

309 1 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . .
309 0 1 2 3 4 5 6 7 8 9 . . . . . . . . . . . . . . . . . . . .
330 . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 . . . . . . . . . .
330 . . . . . . . . . . 0 1 2 3 4 5 6 7 8 9 . . . . . . . . . .
371 . . . . . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1
371 . . . . . . . . . . . . . . . . . . . . 0 1 2 3 4 5 6 7 8 9

Это кажется правильным, но если это так, что за этим стоит линейная алгебра? Я понимаю, что ряды 1людей отбирают. Например, предмет 309включен для базовой линии + девять наблюдений, поэтому он получает четыре 1и так далее. Вторая часть, несомненно, является фактическим измерением: 0для исходного уровня, 1для первого дня лишения сна и т. Д.

Но каковы действительные и и соответствующие им или , в зависимости от того , что уместно?X i Z i = ( J T iX T i ) Z i = ( J T iX T i ) Ji Xi Zi=(JiTXiT) Zi=(JiTXiT)

Здесь есть возможность,

[1111111111,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1111111111,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1111111111]*[11111111110123456789]знак равно

[1111111111,,,,,,,,,,,,,,,,,,,,0123456789,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1111111111,,,,,,,,,,,,,,,,,,,0123456789,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,1111111111,,,,,,,,,,,,,,,,,,,0123456789]

Проблема в том, что это не транспонирование, как, по- lmerвидимому, требует функция, и все еще неясно, каковы правила для создания .Икся

Антони Пареллада
источник
1
Это намного проще, чем ты думаешь. матрица здесь просто (транспонированная) произведение Кронекера единичной матрицы с проектной матрицей. Z
Донни
Спасибо за подсказку. Я буду продолжать работать над осмыслением всех подузлов в скелете линейной алгебры этой функции. Если это произойдет, я пойду дальше и отвечу на свой вопрос, но, хотя я знаю, что это очень просто, соответствие между номенклатурой математических лесов и приложением к любому примеру сбивает с толку.
Антони Пареллада
1
Другим хорошим ресурсом для вас может быть их реализация lme4pureR , которая следует вместе с приведенной выше виньеткой и полностью написана на R. Может быть mkZt()(поиск здесь ) это хорошее начало?
alexforrence

Ответы:

5
  1. Создание матрицы влечет за собой производство 3 уровня ( , и ) каждый из 10 наблюдений или измерений ( ). После кода в оригинальной ссылке в ОП:Jя309330371nrow(sleepstudy[sleepstudy$Subject==309,]) [1] 10

f <- gl(3,10) Ji<-t(as(f,Class="sparseMatrix"))

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

  1. Построение матрицы может быть оказана помощь при использовании функции в качестве ссылки:ИксяgetME

    library(lme4) sleepstudy <- sleepstudy[sleepstudy$Subject %in% c(309, 330, 371), ] rownames(sleepstudy) <- NULL fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)

Xi <- getME(fm1,"mmList")

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

Поскольку нам понадобится транспонирование, а объект Xiне является матрицей, его t(Xi)можно построить так:

t_Xi <- rbind(c(rep(1,30)),c(rep(0:9,3)))

  1. Zя рассчитывается как :Zязнак равно(JяT*ИксяT)

Zi<-t(KhatriRao(t_Ji,t_Xi)):

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

Это соответствует уравнению (6) в оригинальной статье :

Zязнак равно(JяT*ИксяT)Tзнак равно[Jя1TИкся1TJя2TИкся2TJяNTИксяNT]

И чтобы увидеть это, мы можем вместо этого поиграть с усеченными и , представив, что вместо 9 измерений и базовой линии (0), есть только 1 измерение (и базовая линия). Полученные матрицы будут:JяTИксяT

JяTзнак равно[110000001100000011] и .ИксяTзнак равно[111111010101]

А также

JяT*ИксяTзнак равно[(100)(10)(100)(11)(010)(10)(010)(11)(001)(10)(001)(11)]

знак равно[Jя1TИкся1TJя2TИкся2TJя3TИкся3TJя4TИкся4TJя5TИкся5TJя6TИкся6T]

знак равно[110000010000001100000100000011000001] . Какой транспонированная и расширенный приведет к .Zязнак равно[100000110000120000001000001100001200000010000011000012]

  1. Извлечение вектора коэффициентов случайных эффектов можно выполнить с помощью функции:б

b <- getME(fm1,"b")

[1,] -44.1573839
[2,]  -2.4118590
[3,]  32.8633489
[4,]  -0.3998801
[5,]  11.2940350
[6,]   2.8117392

Если мы добавим эти значения к фиксированным эффектам вызова, fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)мы получим перехваты:

205.3016 for 309; 282.3223 for 330; and 260.7530 for 371

и склоны:

2.407141 for 309; 4.419120 for 330; and 7.630739 for 371

значения в соответствии с:

library(lattice)
xyplot(Reaction ~ Days | Subject, groups = Subject, data = sleepstudy, 
       pch=19, lwd=2, type=c('p','r'))

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

  1. Zб можно рассчитать как as.matrix(Zi)%*%b.
Антони Пареллада
источник