Extraire les variances à effet aléatoire de l'objet modèle lme4 mer
j'ai un objet mer qui a des effets fixes et aléatoires. Comment puis-je extraire les estimations de la variance pour les effets aléatoires? Voici une version simplifiée de ma question.
study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
study
cela donne une sortie longue - pas trop longue dans ce cas. Quoi qu'il en soit, Comment puis-je sélectionner explicitement le
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.18 37.124
Residual 960.46 30.991
une partie de la production? Je veux que les valeurs elles-mêmes.
j'ai pris de longues regarde
str(study)
et il n'y a rien! Également vérifié l' fonctions de l'extracteur dans le paquet lme4 en vain. S'il vous plaît aider!
4 réponses
lmer retourne un objet S4, donc cela devrait fonctionner:
remat <- summary(study)@REmat
print(remat, quote=FALSE)
Qui affiche:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.18 37.124
Residual 960.46 30.991
...En général, vous pouvez regarder le source de la print et summary les méthodes pour "mer" des objets:
class(study) # mer
selectMethod("print", "mer")
selectMethod("summary", "mer")
Certains des autres réponses sont réalisables, mais je prétends que la meilleure réponse est d'utiliser la méthode d'accesseur qui est conçu pour cela -- VarCorr (c'est le même que dans , le nlme paquet).
UPDATE dans les versions récentes de lme4 (version 1.1-7, mais tout ci-dessous est probablement applicable pour les versions >= 1.0), VarCorr est plus flexible qu'avant, et devrait faire tout ce que vous voulez sans jamais recourir à la pêche autour à l'intérieur le modèle ajusté objet.
library(lme4)
study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
VarCorr(study)
## Groups Name Std.Dev.
## Subject (Intercept) 37.124
## Residual 30.991
Par défaut VarCorr() affiche les déviations standard, mais vous pouvez obtenir des variances à la place si vous préférez:
print(VarCorr(study),comp="Variance")
## Groups Name Variance
## Subject (Intercept) 1378.18
## Residual 960.46
(comp=c("Variance","Std.Dev.") imprimer).
Pour plus de flexibilité, vous pouvez utiliser le as.data.frame méthode pour convertir le VarCorr objet, qui donne la variable de regroupement, la ou les variables d'effet, et la variance / covariance ou l'écart-type / corrélations:
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
Enfin, la forme brute de la VarCorr object (que vous ne devriez probablement pas vous embêter si vous n'en avez pas besoin) est une liste de matrices de variance-covariance avec des informations supplémentaires (redondantes) encodant les déviations et les corrélations standard, ainsi que des attributs ("sc") donnant l'écart-type résiduel et précisant si le modèle a un paramètre d'échelle estimé ("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)
un exemple:
> 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
La fin.