У меня есть объект 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 безрезультатно. Пожалуйста, помогите!
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")
Некоторые из других ответов работоспособны, но я утверждаю, что лучший ответ заключается в использовании метода доступа, который предназначен для этого - VarCorr (это то же самое, что и в предшественнике lme4, nlme пакет).
UPDATE в последних версиях 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
##
> attributes(summary(study))$REmat
Groups Name Variance Std.Dev.
"Subject" "(Intercept)" "1378.18" "37.124"
"Residual" "" " 960.46" "30.991"
Попробуйте
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
Конец.