[R-es] Comparaciones múltiples en ANOVA anidadp
José Trujillo Carmona
trujillo en unex.es
Dom Feb 5 18:31:52 CET 2012
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
Más información sobre la lista de distribución R-help-es