2011-12-15 7 views
36

Tengo un objeto mer que tiene efectos fijos y aleatorios. ¿Cómo extraigo las estimaciones de varianza para los efectos aleatorios? Aquí hay una versión simplificada de mi pregunta.Extraer variaciones de efectos aleatorios del objeto de modelo lme4 mer

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

Esto proporciona una salida larga, no demasiado larga en este caso. De todos modos, ¿cómo selecciono explícitamente el

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

parte de la salida? Quiero los valores en sí mismos.

yo he tomado largas miradas en

str(study) 

y no hay nada allí! También verifiqué todas las funciones de extractor en el paquete lme4 en vano. ¡Por favor ayuda!

Respuesta

12

lmer devuelve un objeto S4, así que esto debería funcionar:

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

que imprime:

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

... En general, se puede ver en la fuente de la print y summary métodos para objetos "mer":

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

Si quiere los valores, entonces VarCorr() es mucho más eficiente. Eche un vistazo a la publicación de Ben Bolker – Thierry

+1

esto está algo desactualizado ahora (aunque la pregunta original se refiere a "mer objects", que están asociados por definición con pre-1.0 'lme4' - la clase ahora se llama' merMod' –

-6

Trate

attributes(study) 

Como un ejemplo:

> 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 

El extremo.

+2

Err, su código utiliza 'lm()', las preguntas eran sobre 'lmer()' que * no * es lo mismo. –

+0

Así es, esta era una forma general de acceder a los objetos R. – abcde123483

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

Podría estar equivocado, ya que no parece haber atributos 'REmat' en' (resumen (estudio)) 'más – blazej

62

Algunas de las otras respuestas son viables, pero afirmo que la mejor respuesta es usar el método de acceso diseñado para esto - VarCorr (este es el mismo que en el predecesor lme4, el paquete nlme).

ACTUALIZACIÓN en las versiones recientes de lme4 (versión 1.1-7, pero todo a continuación es probablemente aplicable a versiones> = 1.0), VarCorr es más flexible que antes, y se debe hacer todo lo que desee sin tener que recurrir a la pesca en torno dentro del objeto modelo ajustado.

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

Por defecto desviaciones VarCorr() impresiones estándar, pero se puede obtener variaciones en cambio, si lo prefiere:

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

(comp=c("Variance","Std.Dev.") imprimirá ambos).

Para mayor flexibilidad, se puede utilizar el método de as.data.frame para convertir el objeto VarCorr, que da a la variable de agrupación, variable (s) efecto, y la varianza/covarianza o desviación/correlaciones estándar:

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 

Finalmente, la forma cruda del objeto VarCorr (que probablemente no debería meterse con usted si no es necesario) es una lista de matrices de varianza-covarianza con información adicional (redundante) que codifica las desviaciones y correlaciones estándar, así como también como atributos ("sc") dando la desviación estándar residual y especificando si el modelo tiene un parámetro de escala estimado ("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 
## 
+0

VarCorr parece proporcionar solo las desviaciones estándar, no las estimaciones de varianza, que es lo que, en general, las personas desean informar ¿verdad? – user1320502

+3

(1) es bastante fácil cuadrar las desviaciones estándar; (2) 'imprimir (VarCorr (fitted_model)), comp = "Varianza") 'o' como.data.frame (VarCorr (fitted_model)) 'recuperará fácilmente las varianzas; (3) las variaciones de informes frente a las desviaciones estándar dependen del contexto; por lo general, prefiero las variaciones si intento pensar acerca de var descomposición/proporción explicada y desarrolladores estándar si se trata de comparar con la magnitud de los efectos fijos –

+0

Gracias por tu comentario Ben, muy útil! – user1320502

Cuestiones relacionadas