[R-es] Comparaciones múltiples en ANOVA anidadp

José Trujillo Carmona trujillo en unex.es
Sab Feb 4 20:45:38 CET 2012


Explico un poco más el problema con lmer:

Si utilizo lmer, todo parece funcionar, pero el análisis de la varianza 
no es consistente con el de aov:

 > AnovaModel.4 <- lmer(VR ~ Tratamiento + (1|Tratamiento/Individuo), 
data=Medidas)

 > Pares <- glht(AnovaModel.3, linfct = mcp(Tratamiento = "Tukey"))
 > summary(Pares)

      Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = VR ~ Tratamiento + (1 | Tratamiento/Individuo),  
data = Medidas)

Linear Hypotheses:
            Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 -0.14933    0.93547  -0.160    1.000
3 - 1 == 0 -0.13078    0.93547  -0.140    1.000
4 - 1 == 0  0.24593    0.93547   0.263    0.999
5 - 1 == 0 -1.08025    0.93547  -1.155    0.777
3 - 2 == 0  0.01855    0.93547   0.020    1.000
4 - 2 == 0  0.39526    0.93547   0.423    0.993
5 - 2 == 0 -0.93091    0.93547  -0.995    0.858
4 - 3 == 0  0.37671    0.93547   0.403    0.994
5 - 3 == 0 -0.94947    0.93547  -1.015    0.849
5 - 4 == 0 -1.32618    0.93547  -1.418    0.616
(Adjusted p values reported -- single-step method)

 > anova(AnovaModel.4)
Analysis of Variance Table
             Df Sum Sq Mean Sq F value
Tratamiento  4 2.4993 0.62482  0.5819

 > AnovaModel.2 <- aov(VR ~ Tratamiento + Error(Individuo),data=Medidas)

 > summary(AnovaModel.2)

Error: Individuo
             Df Sum Sq Mean Sq F value Pr(>F)
Tratamiento  4  9.166  2.2915   3.473 0.0502 .
Residuals   10  6.599  0.6599
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30  36.35   1.212

No comprendo como habría de escribir la función en lmer para que fuese 
equivalente al de aov que creo congruente con los textos:

Bioestaticas Analysis de Zar
Bimetry de Sokal y Rohlf
Bioestadística de Steel y Torrie


Reitero mi agradecimiento por aguantar



El 04/02/12 20:05, José Trujillo Carmona escribió:
> Dispongo de un experimento en el que cinco tratamientos ha sido 
> aplicados a cinco grupos de voluntarios.
>
> En cada grupo había tres personas y a cada persona se le tomaron 3 
> medidas,
>
> En total dispongo de 45 medidas, pero evidentemente no son 
> independientes entre sí. Si no tomo en cuenta que las medidas de la 
> misma persona son más parecidas entre sí (bloques anidados) estaría 
> incurriendo en pseudoreplicación.
>
> Los datos y el análisis se podrían hacer de la siguiente forma:
>
> Generación de datos para ejemplo:
> > Medidas <- as.data.frame(matrix(rnorm(45*1, mean=10, sd=1), ncol=1))
> > colnames(Medidas) <- "VR"
> > Medidas$Tratamiento <- as.factor(rep(1:5,rep(9,5)))
> > Medidas$Individuo <- as.factor(rep(1:15,rep(3,15)))
>
> El análisis de la varianza sugerido podría ser:
>
> > AnovaModel.1 <-aov(VR ~ Tratamiento + Tratamiento/Individuo, 
> data=Medidas)
> > summary(AnovaModel.1)
>                       Df Sum Sq Mean Sq F value Pr(>F)
> Tratamiento            4   1.15  0.2871    0.25  0.907
> Tratamiento:Individuo 10  12.39  1.2395    1.08  0.407
> Residuals             30  34.44  1.1479
>
>
> Pero este análisis incurre en evidente pseudoreplicación: El valor F 
> es el cociente entre el Cuadrado Medio del Tratamiento y los Residuos 
> como si las 45 mediciones fuesen independientes entre sí.
>
> Para tener en cuenta que la variabilidad inducida en la respuesta por 
> los Tratamientos debe ser medida entre individuos y no entre 
> mediciones el análisis procedente podría ser:
>
> > AnovaModel.2 <- aov(VR ~ Tratamiento + Error(Individuo),data=Medidas)
> > summary(AnovaModel.2)
>
> Error: Individuo
>             Df Sum Sq Mean Sq F value Pr(>F)
> Tratamiento  4  1.148  0.2871   0.232  0.914
> Residuals   10 12.395  1.2395
>
> Error: Within
>           Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 30  34.44   1.148
>
> Pero el objeto AnovaModel.2 no es un objeto aov sino aovlist:
> > attr(AnovaModel.2,"class")
> [1] "aovlist" "listof"
>
> Como consecuencia no puedo utilizarlo ni para la función glht ni para 
> la función TukeyHSD.
>
> Puedo "extraer" un aov:
>
> > ls.str(pat="AnovaModel.2")
> AnovaModel.2 : List of 3
>  $ (Intercept):List of 9
>  $ Individuo  :List of 9
>  $ Within     :List of 6
> > AnovaModel.3<-AnovaModel.2$Individuo
> > attr(AnovaModel.3,"class")
> [1] "aov" "lm"
>
> A este modelo ya puedo pedirle límites de confianza adecuados para los 
> parámetros del modelo (diferencias entre Tratamientos):
>
> > confint(AnovaModel.3)
>                      2.5 %    97.5 %
> Tratamiento[T.2] -1.485182 0.8535804
> Tratamiento[T.3] -1.072891 1.2658714
> Tratamiento[T.4] -1.258890 1.0798723
> Tratamiento[T.5] -1.453857 0.8849054
>
> Pero el nuevo modelo sigue siendo inválido para glht o para TuleyHSD:
>
> > Pares <- glht(AnovaModel.3, linfct = mcp(Tratamiento = "Tukey"))
> [1] ERROR:   no 'model.matrix' method for 'model' found!
>
> > TukeyHSD(AnovaModel.3, "Tratamiento")
> [2] ERROR:  this fit does not inherit from "lm"
>
> Para no alargarme en mis indagaciones. Para los análisis de la 
> varianza podría haber utilizado las funciones lme del paquete nlme o 
> lmer del paquete lme4, pero tampoco proporcionan objetos válidos para 
> glht o TukeyHSD.
>
> ¿Alguien sabe como abordar las comparaciones múltiples con medidas 
> repetidas o anova anidado?
>
> Gracias de antemano.
>

-- 
_____---^---_____

Univ. de Extremadura
Dept. Matemáticas.
Despacho B29
Tf: + 34 924 289 300
Ext. 86823



Más información sobre la lista de distribución R-help-es