Estrarre le varianze degli effetti casuali dall'oggetto modello lme4 mer

Question

Ho un oggetto mer che ha effetti fissi e casuali. Come estrarre le stime della varianza per gli effetti casuali? Ecco una versione semplificata della mia domanda.

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

Questo dà un output lungo, non troppo lungo in questo caso. In ogni caso, come faccio a selezionare esplicitamente lo

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

parte dell'uscita? Voglio i valori stessi.

mi hanno preso sguardi lunghi a

str(study) 

e non c'è niente lì! Inoltre ha controllato tutte le funzioni estrattore nel pacchetto lme4 senza alcun risultato. Per favore aiuto!

Answer 1

lmer restituisce un oggetto S4, quindi questo dovrebbe funzionare:

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

che stampa:

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

... In generale, si può guardare il sorgente della print e metodi summary per gli oggetti "mer":

class(study) # mer 
selectMethod("print", "mer") 
selectMethod("summary", "mer") 

Answer 2

Prova

attributes(study) 

Per fare un esempio:

> 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 fine.

Answer 3

> attributes(summary(study))$REmat 
Groups  Name   Variance Std.Dev. 
"Subject" "(Intercept)" "1378.18" "37.124" 
"Residual" ""   " 960.46" "30.991" 

Answer 4

Alcune delle altre risposte sono realizzabili, ma sostengo che la risposta migliore è utilizzare il metodo accessor progettato per questo - VarCorr (questo è lo stesso del predecessore lme4, il pacchetto nlme).

UPDATE nelle versioni recenti di lme4 (versione 1.1-7, ma tutto sotto è probabilmente applicabile alle versioni> = 1.0), VarCorr è più flessibile rispetto a prima, e dovrebbe fare tutto quello che vuoi, senza mai ricorrere alla pesca intorno all'interno dell'oggetto modello montato.

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

Per impostazione predefinita deviazioni VarCorr() stampa standard, ma è possibile ottenere variazioni, invece, se si preferisce:

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

(comp=c("Variance","Std.Dev.") stamperà entrambi).

Per una maggiore flessibilità, è possibile utilizzare il metodo as.data.frame per convertire l'oggetto VarCorr, che dà la variabile di raggruppamento, variabili effetto (s), e la varianza/covarianza o standard di deviazione/correlazioni:

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 

Infine, la forma grezza dell'oggetto VarCorr (che probabilmente non dovresti rovinare con te se non devi) è un elenco di matrici di varianza-covarianza con informazioni aggiuntive (ridondanti) che codificano le deviazioni standard e le correlazioni, nonché come attributi ("sc") dando la deviazione standard residua e specificando se il modello ha un parametro di scala stimato ("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 
##