Получение матриц случайных эффектов из смешанной модели

ProgramBox

Получение матриц случайных эффектов из смешанной модели

Post author:admin
Запись опубликована:11 декабря, 2022
Post category:Вопросы по программированию

#r #regression #lme4 #mixed-models #nlme

#r #регрессия #lme4 #смешанные модели #nlme

Вопрос:

В моем приведенном ниже коде мне было интересно, как я могу получить эквивалент out и Ts из lme() объекта в library(nlme) ?

 dat <- read.csv("https://raw.githubusercontent.com/rnorouzian/v/main/mv.l.csv")

library(lme4)

x <- lmer(value ~0   name  (1| School/Student), data = dat,
          control = lmerControl(check.nobs.vs.nRE= "ignore"))

    lwr <- getME(x, "lower")
    theta <- getME(x, "theta")

    out = any(theta[lwr == 0] < 1e-4)  # find this from `x1` object below
     Ts = getME(x, "Tlist")            # find this from `x1` object below


# Fitting the above model using `library(nlme)`:

library(nlme)

x1 <- lme(value ~0   name, random = ~1| School/Student, data = dat)

1. Обратите внимание, что для MWE часто лучше использовать один из встроенных наборов данных, а не связываться с удаленным набором данных без каких-либо объяснений структуры этого набора данных.

2. Кроме того, ваш MWE не должен требовать отключения проверок работоспособности (например check.nobs.vs.nRE ), если только вопрос не касается этих проверок работоспособности

Ответ №1:

Я бы настоятельно рекомендовал прочитать документы! lme4 и nlme используют по своей сути разные методы для подгонки смешанных моделей — lme4 использует штрафную формулировку наименьших квадратов, основанную на более низком коэффициенте Холецкого (тета), а nlme4 использует обобщенную формулировку наименьших квадратов, которая при желании может быть сохранена как коэффициент Холецкого, но их документы дают вам информацию, чтобы получить то, что вам нужноиз внутреннего представления. После этого вам нужно выполнить математические вычисления для преобразования между представлениями.

Если вы это сделаете ?lme , тогда есть строка

Смотрите lmeObject компоненты подгонки

Затем вы делаете ?lmeObject , где находите две многообещающие записи:

apVar
приблизительная ковариационная матрица для коэффициентов дисперсии-ковариации. Если apVar = FALSE в управляющих значениях, используемых при вызове lme , этот компонент NULL .

modelStruct объект, наследуемый от class lmeStruct , представляющий список компонентов модели со смешанными эффектами, таких как reStruct , corStruct , и varFunc объекты.

Ну, на самом деле нам нужны не коэффициенты переменных, а скорее матрицы случайных эффектов. Итак, мы можем посмотреть reStruct . Это намного более гибко в nlme, чем lme4, но обычно это просто матрицы случайных эффектов. Чтобы сделать что-либо сопоставимое с lme4, вам необходимо преобразовать их в их более низкий коэффициент Холецкого. Вот пример использования sleepstudy данных:

 > library("nlme")
> library("lme4")
> 
> data("sleepstudy")
> m_nlme <- lme(fixed=Reaction ~ 1   Days,
                random=~ 1   Days | Subject,
                data=sleepstudy,
                method = "ML")
> m_lme4 <- lmer(Reaction ~ 1   Days   (1   Days|Subject),
                 data=sleepstudy,
                 REML=FALSE)
> 
> re_lme4 <- getME(m_lme4, "Tlist")$Subject
> print(re_lme4)
           [,1]      [,2]
[1,] 0.92919061 0.0000000
[2,] 0.01816575 0.2226432
> 
> re_nlme <- m_nlme$modelStruct$reStruct$Subject
> # entire covariance matrix
> print(re_nlme)
Positive definite matrix structure of class pdLogChol representing
            (Intercept)       Days
(Intercept)  0.86344433 0.01688228
Days         0.01688228 0.04990040
> # get the lower cholesky factor
> re_nlme <- t(chol(re_nlme)) # could also use pdMatrix(re_nlme, TRUE)
> print(re_nlme)
            (Intercept)      Days
(Intercept)  0.92921705 0.0000000
Days         0.01816829 0.2226439

Тета-вектор для lme4 — это просто представление нижнего треугольника нижнего коэффициента Холецкого для данной группирующей переменной. (Для моделей с несколькими группирующими переменными вы просто объединяете их вместе.) Нижний коэффициент Холецкого ограничен, чтобы не иметь записей меньше нуля по диагонали (потому что это соответствовало бы отрицательной дисперсии), и в противном случае не ограничен. Другими словами, диагональные записи получают нижнюю границу при 0, все остальные записи получают нижнюю границу при -Inf .

Итак, в lme4:

 > re_lme4[lower.tri(re_lme4,diag = TRUE)]
[1] 0.92919061 0.01816575 0.22264321
> getME(m_lme4, "theta")
     Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
              0.92919061               0.01816575               0.22264321 
> getME(m_lme4, "lower")
[1]    0 -Inf    0

Мы можем реализовать это для nlme (не самый эффективный способ, но он показывает, как все устроено):

 > lowerbd <- function(x){
    dd <- diag(0, nrow=nrow(x))
    dd[lower.tri(dd)] <- -Inf
    dd[lower.tri(dd, diag=TRUE)]
  }
> lowerbd(re_nlme)
[1]    0 -Inf    0
> lowerbd(re_lme4)
[1]    0 -Inf    0

Обратите внимание, что это одно место, где nlme на самом деле более мощный, чем lme4: весь pdMatrix набор ограничений может создавать разные нижние границы для разных записей (а также, например, ограничивать связь между записями).