[R-es] factores anidados
Olivier Nuñez
onunez en iberstat.es
Jue Sep 10 22:33:29 CEST 2009
Ricardo,
perdona por no haberte contestado antes pero ayer y hoy estuve fuera.
La diferencia que observas entre las dos salidas (la de aov y la lme)
es debida a tu interpretación de estos resultados.
En la salida de lmer
> Random effects:
> Groups Name Variance Std.Dev.
> operador:(muestra:(lote:cc)) (Intercept) 1.5770e-01 3.9712e-01
> muestra:(lote:cc) (Intercept) 0.0000e+00 0.0000e+00
> lote:cc (Intercept) 6.6065e-15 8.1280e-08
> cc (Intercept) 0.0000e+00 0.0000e+00
> Residual 7.3738e-01 8.5871e-01
se te da las estimaciones de las componentes de la varianza:
variabilidad (depurada de la varianza residual) entre niveles de un
factor (cruzado o no).
Mientras que en la salida del aov,
> > summary(aov(y~cc/lote/muestra/operador))
> Df Sum Sq Mean Sq F value Pr(>F)
> cc 1 0.090 0.090 0.1224 0.72737
> cc:lote 2 2.113 1.056 1.4324 0.24478
> cc:lote:muestra 36 33.032 0.918 1.2444 0.20819
> cc:lote:muestra:operador 40 47.935 1.198 1.6252 0.03315 *
> Residuals 80 58.990 0.737
la columna Mean Sq corresponde a una descomposición de la
variabilidad total en componentes ortogonales (en un diseño
equilibrado).
Estas dos descomposiciones no coinciden en general salvo la
estimación de la varianza residuales.
Para simplificar y completar esta aclaración, considero un modelo
sencillo con un solo factor
Y_ij = mu + B_i + eps_ij, donde i=12,..,I y j =1,2,..,J
donde mu es una constante, los B_i son efectos aleatorios con
varianza V_B (efectos de un factor bloque en el diseño) y
los eps_ij son errores aleatorios con varianza V_R (varianza residual).
La esperanza de la varianza estimada en el anova entre los niveles
del factor B es: V_R + J*V_B.
En otras palabras, esta varianza está "contaminada" por la varianza
residual.
Para ilustrar este ejemplo teorico, utilizamos los datos "oats" del
package MASS, se obiene:
> require(MASS)
> require(lme4)
> summary(aov(Y~B,oats))
Df Sum Sq Mean Sq F value Pr(>F)
B 5 15875 3175 5.8031 0.0001683 ***
Residuals 66 36111 547
> lmer(Y~(1|B),oats)
Linear mixed model fit by REML
Formula: Y ~ (1 | B)
Data: oats
AIC BIC logLik deviance REMLdev
668.2 675 -331.1 667.8 662.2
Random effects:
Groups Name Variance Std.Dev.
B (Intercept) 218.99 14.798
Residual 547.13 23.391
Number of obs: 72, groups: B, 6
Y efectivamente, observamos que la "Mean Sq" del factor B es igual a
3175 = 547.13 + J*218.99
donde J=12 es el número de observaciones en cada bloque.
En resumen, si estás interesado en estimar directamente los
componentes de la varianza te aconsejo utilizar la función lmer.
Por otro lado, te aconsejo expresar en tu modelo cuales son los
efectos fijos y los efectos aleatorios con el fin de llevar tu
inferencia de manera más eficientes.
Un saludo. Olivier
-- ____________________________________
Olivier G. Nuñez
Email: onunez en iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es
____________________________________
El 09/09/2009, a las 18:02, Ricardo Bello escribió:
> Ahh, te cuento que resolví el problema a mano, calculando las sumas
> de cuadrados, los cuadrados medios esperados y las
> componentes de la variancia. La tabla de ANOVA es como sigue:
>
> > summary(aov(y~cc/lote/muestra/operador))
> Df Sum Sq Mean Sq F value Pr(>F)
> cc 1 0.090 0.090 0.1224 0.72737
> cc:lote 2 2.113 1.056 1.4324 0.24478
> cc:lote:muestra 36 33.032 0.918 1.2444 0.20819
> cc:lote:muestra:operador 40 47.935 1.198 1.6252 0.03315 *
> Residuals 80 58.990 0.737
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Los resultados anteriores coinciden con los que calculé a mano.
> Lo siguiente son las componentes de la varianza con las funciones
> lme y lmer.
> Si a partir de la tabla ANOVA anterior calculas las componentes de
> la varianza, notarás que no coinciden con las
> obtenidas por medio de lme o lmer, a no ser que no las esté
> calculando como corresponde.
>
> > lmer(y~(1|cc/lote/muestra/operador))
> Linear mixed model fit by REML
> Formula: y ~ (1 | cc/lote/muestra/operador)
> AIC BIC logLik deviance REMLdev
> 448 466.4 -218 432.8 436
> Random effects:
> Groups Name Variance Std.Dev.
> operador:(muestra:(lote:cc)) (Intercept) 1.5770e-01 3.9712e-01
> muestra:(lote:cc) (Intercept) 0.0000e+00 0.0000e+00
> lote:cc (Intercept) 6.6065e-15 8.1280e-08
> cc (Intercept) 0.0000e+00 0.0000e+00
> Residual 7.3738e-01 8.5871e-01
> Number of obs: 160, groups: operador:(muestra:(lote:cc)), 80;
> muestra:(lote:cc), 40; lote:cc, 4; cc, 2
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 5.05125 0.08112 62.27
>
> > lme(y~1,random=~1|cc/lote/muestra/operador)
> Linear mixed-effects model fit by REML
> Data: NULL
> Log-restricted-likelihood: -217.9942
> Fixed: y ~ 1
> (Intercept)
> 5.05125
> Random effects:
> Formula: ~1 | cc
> (Intercept)
> StdDev: 2.181074e-05
> Formula: ~1 | lote %in% cc
> (Intercept)
> StdDev: 2.901981e-05
> Formula: ~1 | muestra %in% lote %in% cc
> (Intercept)
> StdDev: 2.749631e-05
> Formula: ~1 | operador %in% muestra %in% lote %in% cc
> (Intercept) Residual
> StdDev: 0.3971188 0.8587054
> Number of Observations: 160
> Number of Groups:
> cc
> lote %in% cc muestra %in% lote %in% cc
>
> 2
> 4 40
> operador %in% muestra %in% lote %in% cc
> 80
>
>
> espero haber sido claro y desde ya, agradezco tu colaboración.
> ricardo.
>
>
> From: onunez en iberstat.es
> Subject: Re: [R-es] factores anidados
> Date: Wed, 9 Sep 2009 17:40:02 +0200
> To: r_bello en hotmail.com
>
> ..... En otras palabras: ¿Cómo averiguas que es lo que no es
> correcto en los resultados?
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez en iberstat.es
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 09/09/2009, a las 17:27, Ricardo Bello escribió:
>
> Hola Olivier, gracias por tu respuesta pero con la función lme que
> me propones obtengo los mismos resultados que con lmer y éstos, no
> son correctos. Si te sirve mas información te comento que es un
> diseño completamente aleatorizado con cuatro factores, todos ellos
> anidados y aleatorios: operarios dentro de muestras, muestras
> dentro de lotes y lotes dentro de concentraciones. Muchas gracias
> de todos modos.
> ricardo.
>
>
>
>
>
> CC: r-help-es en r-project.org
> From: onunez en iberstat.es
> Subject: Re: [R-es] factores anidados
> Date: Wed, 9 Sep 2009 17:20:12 +0200
> To: r_bello en hotmail.com
>
> Sin más información sobre tu modelo y/o el diseño de tu experimento,
> te aconsejo utilizar lme en vez de lmer:
>
> require(nlme)
>
> lme(y~1, random=~1|cc/ lote/muestra/ operador)
>
> Un saludo. Olivier
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez en iberstat.es
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 09/09/2009, a las 16:46, Ricardo Bello escribió:
>
>
> Tengo la siguiente consulta:
>
> Tengo un modelo lineal con 4 factores anidados y aleatorios
> (concentració n (cc), lote, muestra y operador) en ese orden. Para
> obtener las componentes de la varianza utilizo la función "lmer"
> pero no obtengo los resultados correctos:
>
> lmer(y~(1|cc/ lote/muestra/ operador) ) ; y es la respuesta.
>
> me podrían decir donde está el error??
> gracias.
> ricardo.
>
>
>
>
> _________________________________________________________________
> Revisá tus correos de Hotmail en tu BlackBerry - Clic Aquí
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
> ¡Creá carpetas para todos tus correos! ¡Hotmail te da todo el
> espacio que puedas necesitar!
>
>
> Ingresá a tu Hotmail desde tu Messenger. ¡Windows Live hace tu vida
> más simple!
Más información sobre la lista de distribución R-help-es