[R-es] Comparaciones múltiples en ANOVA anidadp
Luciano Selzer
luciano.selzer en gmail.com
Dom Feb 5 13:57:55 CET 2012
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
>
Más información sobre la lista de distribución R-help-es