[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