[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