Извлечение случайных отклонений эффекта от объекта модели lme4 mer - PullRequest
Новая
       44

Извлечение случайных отклонений эффекта от объекта модели lme4 mer

39 голосов
/ 16 декабря 2011

У меня есть объект mer, который имеет фиксированные и случайные эффекты. Как извлечь оценки дисперсии для случайных эффектов? Вот упрощенная версия моего вопроса. 

study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
study

Это дает длинный вывод - в данном случае не слишком длинный. Во всяком случае, как мне явно выбрать 

Random effects:
Groups   Name        Variance Std.Dev.
Subject  (Intercept) 1378.18  37.124  
Residual              960.46  30.991

часть вывода? Я хочу сами значения.

Я долго смотрел на 

str(study)

и там ничего нет! Также проверил все функции экстрактора в пакете lme4 безрезультатно. Пожалуйста, помогите!

Ответы [ 5 ]

72 голосов
/ 16 декабря 2011

Некоторые другие ответы работоспособны, но я утверждаю, что лучший ответ - использовать метод доступа, разработанный для этого - VarCorr (это то же самое, что и в предшественнике lme4, nlme пакет).

ОБНОВЛЕНИЕ в последних версиях lme4 (версия 1.1-7, но все, что ниже, вероятно, применимо к версиям> = 1.0), VarCorr более гибкочем раньше, и должен делать все, что вы хотите, не прибегая к рыбалке внутри приспособленного модельного объекта.

library(lme4)
study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
VarCorr(study)
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 37.124  
##  Residual             30.991

По умолчанию VarCorr() печатает стандартные отклонения, но вместо этого можно получить отклонения, если вы предпочитаете: 

print(VarCorr(study),comp="Variance")
##  Groups   Name        Variance
##  Subject  (Intercept) 1378.18 
##  Residual              960.46

(comp=c("Variance","Std.Dev.") напечатает оба).

Для большей гибкости вы можете использовать метод as.data.frame для преобразования объекта VarCorr, который дает группирующую переменную, переменную (и) эффекта и дисперсию / ковариацию или стандартное отклонение / корреляции: 

as.data.frame(VarCorr(study))
##        grp        var1 var2      vcov    sdcor
## 1  Subject (Intercept) <NA> 1378.1785 37.12383
## 2 Residual        <NA> <NA>  960.4566 30.99123

Наконец, необработанная форма объекта VarCorr (с которым вам, вероятно, не стоит связываться, если вам не нужно), представляет собой список матриц дисперсии-ковариации с дополнительной (избыточной) информацией, кодирующей стандарт.отклонения и корреляции, а также атрибуты ("sc"), дающие остаточное стандартное отклонение и указывающие, имеет ли модель оценочный масштабный параметр ("useSc"). 

unclass(VarCorr(fm1))
## $Subject
##             (Intercept)      Days
## (Intercept)  612.089748  9.604335
## Days           9.604335 35.071662
## attr(,"stddev")
## (Intercept)        Days 
##   24.740448    5.922133 
## attr(,"correlation")
##             (Intercept)       Days
## (Intercept)  1.00000000 0.06555134
## Days         0.06555134 1.00000000
## 
## attr(,"sc")
## [1] 25.59182
## attr(,"useSc")
## [1] TRUE
##
14 голосов
/ 16 декабря 2011

lmer возвращает объект S4, поэтому это должно работать: 

remat <- summary(study)@REmat
print(remat, quote=FALSE)

Какие отпечатки: 

 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1378.18  37.124  
 Residual              960.46  30.991

... В общем, вы можете посмотреть на источник методов print и summary для объектов "mer": 

class(study) # mer
selectMethod("print", "mer")
selectMethod("summary", "mer")
5 голосов
/ 16 декабря 2011 
> attributes(summary(study))$REmat
 Groups     Name          Variance  Std.Dev.
 "Subject"  "(Intercept)" "1378.18" "37.124"
 "Residual" ""            " 960.46" "30.991"
0 голосов
/ 20 марта 2019

Этот ответ в значительной степени основан на ответе @Ben Bolker, но если люди новички в этом и хотят сами значения, а не просто распечатывают значения (как, кажется, хотел OP), вы можете извлечьследующие значения:

Преобразование объекта VarCorr во фрейм данных. 

re_dat = as.data.frame(VarCorr(study))

Затем доступ к каждому отдельному значению: 

int_vcov = re_dat[1,'vcov']
resid_vcov = re_dat[2,'vcov']

С помощью этого методастроки и столбцы в созданной вами рамке даты) вы можете получить доступ к любым значениям, которые вам нужны.

0 голосов
/ 16 декабря 2011

Попробуйте 

attributes(study)

Например: 

> women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164

> lm1 <- lm(height ~ weight, data=women)
> attributes(lm1)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"

> lm1$coefficients
(Intercept)      weight 
 25.7234557   0.2872492 

> lm1$coefficients[[1]]

[1] 25.72346


> lm1$coefficients[[2]]

[1] 0.2872492

Конец.

...