[R-es] Comparaciones múltiples en ANOVA anidadp
José Trujillo Carmona
trujillo en unex.es
Lun Feb 6 12:27:27 CET 2012
Gracias Olivier.
Pero el análisis de la varianza sigue sin ser el que cabría esperar de
los Cuadrados Medios del modelo mixto:
Cuadrados medios esperados (Textos ya citados):
Suma (media(x)_i - media(x)_T)^2 / (a-1) -> var(epsilon) + n sigma^2_B +
n b (suma alfa_i)^2 / (a-1)
Suma (media(x)_ij - media(x)_i)^2 / (ba-a) -> var(epsilon) + n sigma^2_B
Suma (x_ijl - media(x)_ij)^2 / (abn-ab) -> var(epsilon)
Las ejecuciones son cada vez distintas por la ejecución de rnorm()
#> AnovaModel.2 <- aov(VR ~ Tratamiento + Error(Individuo),data=Medidas)
#> AnovaModel.3 <- lmer(VR ~ Tratamiento+(1|Tratamiento/Individuo),
data=Medidas)
#> fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento),
data=Medidas)
#> summary(AnovaModel.2)
Error: Individuo
Df Sum Sq Mean Sq F value Pr(>F)
Tratamiento 4 4.457 1.114 1.085 0.415
Residuals 10 10.270 1.027
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30 22.85 0.7618
#> > anova(AnovaModel.3)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Tratamiento 4 1.2155 0.30388 0.367
#> anova(fit.lmer)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Tratamiento 4 4.4568 1.1142 1.3455
En realidad el valor de F se aleja aún más del valor que estoy buscando.
¿Quizás es por la utilización de estimas de Máxima Verosimilitud?
¿Se le pueden pedir estimas Minimo Cuadráticas a lmer?
Gracias.
UN PROBLEMA que surgiría con la propuesta de Luciano: promediar los
individuos, es que si se pierden datos tendría que ponderar sus medias o
algo parecido que tendría que calcular previamente al análisis. El
problema se complicaría. Supongo que el modelo tendría este problema
El 06/02/12 01:17, Olivier Nuñez escribió:
> José,
>
> la expansión de (1|Tratamiento/Individuo) es (1|Individuo:Tratamiento)
> + (1|Tratamiento)
>
> SI como creo, sólo te interesa estimar los efectos fijos de
> "Tratamiento", te sugiero escribir
>
> fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento),
> data=Medidas);anova(fit.lmer)
>
> Un saludo. Olivier
>
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez en iberstat.es
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 05/02/2012, a las 18:31, José Trujillo Carmona escribió:
>
>> Sí, promediar cada individuo sería una solución. La solución parece
>> un poco chapuza, pero es perfectamente correcta.
>>
>> Como expliqué no es la opción utilizada en los textos consultados y
>> me parece extraño que no pueda escribir el modelo completo en R.
>>
>> Respecto a que en (1|Tratamiento/Individuo) Tratamiento sea
>> aleatorio, no sé muy bien lo que quieres decir.
>>
>> El resultado de "VR ~ (1|Tratamiento) + (1|Tratamiento/Individuo)",
>> donde especifico que Tratamiento es aleatorio no es el mismo que en
>> "VR ~ Tratamiento + (1|Tratamiento/Individuo)".
>>
>> Yo intento decir que Individuo sí es un factor aleatorio y que está
>> dentro de Tratamiento; si esto no se dice con
>> (1|Tratamiento/Individuo) es precisamente el motivo de mi pregunta
>> ¿Cómo lo digo con lmer o con lme? Con lme solo me interesa en el caso
>> de que después pueda utilizar glht con el modelo de lme.
>>
>>
>>
>>
>>
>> El 05/02/12 13:57, Luciano Selzer escribió:
>>> José,
>>> 1. de todas formas si cada individuo se uso solo en un tratamiento no
>>> es necesario poner Tratamiento/Individuo. Al poner esa forma estás
>>> especificando que Tratamiento es un factor aleatorio
>>> 2. El método anova.lmer va a obviar los factores aleatorios, no es que
>>> no lo este considerando sino que de esa forma no va a aparecer en el
>>> analisis.
>>> 3.Si no es de tu interés estimar la variabilidad dentro de cada
>>> individuo puedes promediar sus mediciones y listo. Así podes usar aov
>>> sin necesidad de usar el término de error.
>>>
>>> Un Saludo
>>> Luciano
>>>
>>>
>>>
>>> El día 5 de febrero de 2012 07:58, José Trujillo Carmona
>>> <trujillo en unex.es> escribió:
>>>
>>>> Es que lo individuos son únicos para cada tratamiento.
>>>>
>>>> Los tres individuos de cada tratamiento no reciben los otros
>>>> tratamientos.
>>>>
>>>> Por eso no puedo utilizar simplemente:
>>>>
>>>> Anova.Model<- aov(VR ~ Tratamiento + Individuo, data=Medidas)
>>>>
>>>> Sino que debería señalar que la variabilidad entre Tratamientos la
>>>> dan los
>>>> individuos, que son distintos en cada tratamiento (función 'Error' en
>>>> 'aov').
>>>>
>>>>
>>>> En los datos tengo:
>>>>
>>>>
>>>>
>>>>> xtabs(~Individuo+Tratamiento, data=Medidas)
>>>>>
>>>> Tratamiento
>>>> Individuo 1 2 3 4 5
>>>> 1 3 0 0 0 0
>>>> 2 3 0 0 0 0
>>>> 3 3 0 0 0 0
>>>> 4 0 3 0 0 0
>>>> 5 0 3 0 0 0
>>>> 6 0 3 0 0 0
>>>> 7 0 0 3 0 0
>>>> 8 0 0 3 0 0
>>>> 9 0 0 3 0 0
>>>> 10 0 0 0 3 0
>>>> 11 0 0 0 3 0
>>>> 12 0 0 0 3 0
>>>> 13 0 0 0 0 3
>>>> 14 0 0 0 0 3
>>>> 15 0 0 0 0 3
>>>>
>>>> Los individuos están "anidados" respecto de los tratamientos y la
>>>> variabilidad entre tratamientos, aunque no haya diferencias
>>>> significativas
>>>> entre los tratamientos, al menos ha de ser tan grande como la
>>>> variabilidad
>>>> entre individuos. Precisamente porque al no estar repetidos los
>>>> individuos
>>>> en los tratamientos no puedo descontar la variabilidad de los
>>>> individuos en
>>>> los tratamientos.
>>>>
>>>> Las observaciones, al estar repetidas las de cada individuo, pueden
>>>> variar
>>>> menos que los individuos y si hay diferencias entre individuos, hay
>>>> que
>>>> evitar en el contraste que se confunda con las posibles diferencias
>>>> entre
>>>> tratamientos.
>>>>
>>>> Por otra parte el análisis que tu propones incurre en
>>>> pseudorreplicación
>>>> extrema: en el Anova obvia por completo a los individuos como si
>>>> fuese un
>>>> factor inexistente:
>>>>
>>>>
>>>>> AnovaModel.5<-lmer(VR ~ Tratamiento + (1|Individuo), data = Medidas)
>>>>>
>>>>
>>>>> anova(AnovaModel.5)
>>>>>
>>>> Analysis of Variance Table
>>>> Df Sum Sq Mean Sq F value
>>>> Tratamiento 4 0.57051 0.14263 0.245
>>>>
>>>> Si hago el análisis con un único factor, como si todas las
>>>> observaciones
>>>> fuesen independientes:
>>>>
>>>>
>>>>> AnovaModel.6<-aov(VR ~ Tratamiento, data=Medidas)
>>>>>
>>>>
>>>>> summary(AnovaModel.6)
>>>>>
>>>> Df Sum Sq Mean Sq F value Pr(>F)
>>>> Tratamiento 4 0.571 0.1426 0.245 0.911
>>>> Residuals 40 23.287 0.5822
>>>>
>>>> Es el mismo análisis de la varianza (la misma F).
>>>>
>>>>
>>>> Gracias Luciano por el interés.
>>>>
>>>>
>>>> El 05/02/12 04:14, Luciano Selzer escribió:
>>>>
>>>>> Hola José, puedo estár equivocado porque son las 12 de la noche ...
>>>>>
>>>>> Pero creo que el modelo equivalente en lmer es
>>>>> lmer(VR ~ Tratamiento + (1|Individuo), data = Medidas)
>>>>>
>>>>> El otro modelo implica que tus individuos son únicos para cada
>>>>> tratamiento.
>>>>>
>>>>> Un saludo
>>>>>
>>>>> Luciano
>>>>>
>>>>>
>>>>>
>>>>> El día 4 de febrero de 2012 16:45, José Trujillo Carmona
>>>>> <trujillo en unex.es> escribió:
>>>>>
>>>>>
>>>>>> 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
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-help-es mailing list
>>>>>> R-help-es en r-project.org
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>>>>>
>>>>>>
>>>>
>>>> --
>>>> _____---^---_____
>>>>
>>>> Univ. de Extremadura
>>>> Dept. Matemáticas.
>>>> Despacho B29
>>>> Tf: + 34 924 289 300
>>>> Ext. 86823
>>>>
>>>>
>>
>> --
>> _____---^---_____
>>
>> Univ. de Extremadura
>> Dept. Matemáticas.
>> Despacho B29
>> Tf: + 34 924 289 300
>> Ext. 86823
>>
>> _______________________________________________
>> R-help-es mailing list
>> R-help-es en r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-help-es
--
_____---^---_____
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