[R-es] Comparaciones múltiples en ANOVA anidadp
José Trujillo Carmona
trujillo en unex.es
Dom Feb 5 11:58:42 CET 2012
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