José Trujillo Carmona
2012-Feb-04 19:05 UTC
[R-es] Comparaciones múltiples en ANOVA anidadp
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
José Trujillo Carmona
2012-Feb-04 19:45 UTC
[R-es] Comparaciones múltiples en ANOVA anidadp
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
José Trujillo Carmona
2012-Feb-06 11:27 UTC
[R-es] Comparaciones múltiples en ANOVA anidadp
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
José,
perdona si no interpreto bien tu correo,
¿pero es la diferencia de salida de fit,lmer y AnovaModel.2 que no te
cuadra?
¿Cual es tu valor de referencia para la suma de cuadrados?
Tenemos que
#> 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
#> anova(fit.lmer)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Tratamiento 4 4.4568 1.1142 1.3455
que resultan en realidad muy parecidas.
Ten en cuenta, que los métodos de cálculos de ambos enfoques (aov y
lmer) son distintos.
Además, lmer está calculado con REML (restricted maximum likelihood).
Puedes probar, REML=FALSE para que lmer utilicé el Maximum Likelihood
a secas:
fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento),
data=Medidas, REML=FALSE)
Un saludo. Olivier
-- ____________________________________
Olivier G. Nuñez
Email: onunez@iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es
____________________________________
El 06/02/2012, a las 12:27, José Trujillo Carmona escribió:
> 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@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@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@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@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@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
>
> _______________________________________________
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
--
____________________________________
Olivier G. Nuñez
Email: onunez@unex.es
http://kolmogorov.unex.es/~onunez
Tel : +34 663 03 69 09
Departamento de Matemáticas
Universidad de Extremadura
____________________________________
[[alternative HTML version deleted]]
José Trujillo Carmona
2012-Feb-06 17:44 UTC
[R-es] Comparaciones múltiples en ANOVA anidadp
El valor que no cuadra es el Cuadrado Medio del denominador de F. Para fit.lmer es: 1.3455/1.1142=1,207592892 Y para AnovaModel.2 es: 1.027 (MS Residuals en Error(Individuo)) La diferencia en las estimas es casi un 18% mayor y eso tratándose de un caso donde no hay varianzas adicionales ni el diseño está desequilibrado. La varianza que para el error estima lmer está entre la estimación que aov hace del error puro (1.284) y la que hace de la varianza entre individuos (1.027). La situación no cambia sustancialmente si utilizo REML=F. Los números sí, pero la idea no. Supongo que la diferencia viene, como apuntas, por la utilización de la estimación ML. AOV en cambio utiliza la estimación Mínimo Cuadrática de los parámetros del modelo lineal y a partir de ahí la estimación insesgada de las varianzas total, residual y parciales. El problema es que la estimación del MS del Error es crítico para el test de Tukey. Quizás tenga que olvidarme de la estimación mínimo cuadrática si tengo que habérmelas con modelos desequlibrados y estudiarme la estimación máximo verosímil. Gracias. El 06/02/12 14:33, Olivier Nuñez escribió:> José, > > perdona si no interpreto bien tu correo, > ¿pero es la diferencia de salida de fit,lmer y AnovaModel.2 que no te > cuadra? > ¿Cual es tu valor de referencia para la suma de cuadrados? > > Tenemos que > #> 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 > > #> anova(fit.lmer) > Analysis of Variance Table > Df Sum Sq Mean Sq F value > Tratamiento 4 4.4568 1.1142 1.3455 > > que resultan en realidad muy parecidas. > Ten en cuenta, que los métodos de cálculos de ambos enfoques (aov y > lmer) son distintos. > Además, lmer está calculado con REML (restricted maximum likelihood). > Puedes probar, REML=FALSE para que lmer utilicé el Maximum Likelihood > a secas: > > fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento), > data=Medidas, REML=FALSE) > > Un saludo. Olivier > -- ____________________________________ > > Olivier G. Nuñez > Email: onunez en iberstat.es > Tel : +34 663 03 69 09 > Web: http://www.iberstat.es > > ____________________________________ > > > > > El 06/02/2012, a las 12:27, José Trujillo Carmona escribió: > >> 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 >> >> _______________________________________________ >> 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
José Trujillo Carmona
2012-Feb-07 11:09 UTC
[R-es] Comparaciones múltiples en ANOVA anidadp
Bueno, doy por resuelto el tema con la penúltima sugerencia de Olivier.
Creo que el siguiente ejemplo prueba que la sugerencia es buena:
> VI<-rnorm(15, mean=0, sd=1)
> VRs<-rnorm(45, mean=0, sd=0.3)
> patient <- as.factor(rep(1:15,rep(3,15)))
> Mesures <- as.data.frame(matrix(10+VI[patient]+VRs, ncol=1))
> colnames(Mesures) <- "VR"
> Mesures$trat <- as.factor(rep(1:5,rep(9,5)))
> Mesures$patient <- patient
En este ejemplo ha de haber diferencias significativas entre pacientes
(Componente VI de VR) pero no entre tratamientos.
Efectivamente en el primer método aparece la pseudoreplicación:
> AnovaModel.1 <-aov(VR ~ trat + trat/patient, data=Mesures)
> summary(AnovaModel.1)
Df Sum Sq Mean Sq F value Pr(>F)
trat 4 22.36 5.589 65.23 2.29e-14 ***
trat:patient 10 36.27 3.627 42.34 6.13e-15 ***
Residuals 30 2.57 0.086
---
Signif. codes: 0 ''***'' 0.001 ''**'' 0.01
''*'' 0.05 ''.'' 0.1 '' '' 1
Pero no en el segundo método:
> AnovaModel.2 <- aov(VR ~ trat + Error(patient),data=Mesures)
> summary(AnovaModel.2)
Error: patient
Df Sum Sq Mean Sq F value Pr(>F)
trat 4 22.36 5.589 1.541 0.263
Residuals 10 36.27 3.627
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30 2.57 0.08568
Ninguno de los tres métodos propuestos con lmer incurre en
pseudorreplicación:
> AnovaModel.3 <- lmer(VR ~ trat+(1|trat/patient), data=Mesures)
> AnovaModel.4 <- lmer(VR ~ trat+(1|patient:trat), data=Mesures)
> AnovaModel.5 <- lmer(VR ~ trat+(1|patient:trat), data=Mesures,
REML=FALSE)
> anova(AnovaModel.3)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
trat 4 0.49632 0.12408 1.4482
> anova(AnovaModel.4)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
trat 4 0.5269 0.13173 1.5374
> anova(AnovaModel.5)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
trat 4 0.79108 0.19777 2.3082
Todos son equivalentes. El último parece el más robusto y no incurre
desde luego en pseudorreplicación.
El penúltimo es el más parecido al aov, pero solo es un caso no una
demostración.
Supongo que en todos la diferencia con aov viene dada por utilizar ML en
las estimaciones en lugar de MC.
Utilizar ML para alumnos de Agronomía que tienen que abordar diseño de
experimentos sin más conocimientos que un cuatrimestre de Estadística
Básica se me antoja arduo, muy arduo. No tengo tiempo ni para el cálculo
más elemental de ML que sea comprensible sin sacrificar otros
contenidos; temo convertir las estimaciones en pura magia, sin
explicaciones.
Además el test no ofece el MS del error, aunque es evidente su cálculo a
partir de F y el MS del tratamiento, la salida no es nada elegante.
Gracias a Luciano y a Olivier por su interés y ruego disculpas al resto
de la lista a los que no haya interesado el tema.
El 06/02/12 15:26, Olivier Nuñez escribió:> José,
>
> perdona si no interpreto bien tu correo,
> ¿pero es la diferencia de salida de fit,lmer y AnovaModel.2 que no te
> cuadra?
> ¿Cual es tu valor de referencia para la suma de cuadrados?
>
> Tenemos que
> #> 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
>
> #> anova(fit.lmer)
> Analysis of Variance Table
> Df Sum Sq Mean Sq F value
> Tratamiento 4 4.4568 1.1142 1.3455
>
> que resultan en realidad muy parecidas.
> Ten en cuenta, que los métodos de cálculos de ambos enfoques (aov y
> lmer) son distintos.
> Además, lmer está calculado con REML (restricted maximum likelihood).
> Puedes probar, REML=FALSE para que lmer utilicé el Maximum Likelihood
> a secas:
>
> fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento),
> data=Medidas, REML=FALSE)
>
> Un saludo. Olivier
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez@iberstat.es <mailto:onunez@iberstat.es>
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 06/02/2012, a las 12:27, José Trujillo Carmona escribió:
>
>> 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@iberstat.es <mailto:onunez@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@unex.es <mailto:trujillo@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@unex.es
<mailto:trujillo@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@r-project.org
<mailto:R-help-es@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@r-project.org <mailto:R-help-es@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
>>
>> _______________________________________________
>> R-help-es mailing list
>> R-help-es@r-project.org <mailto:R-help-es@r-project.org>
>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
>
> --
> ____________________________________
>
> Olivier G. Nuñez
> Email: onunez@unex.es <mailto:onunez@unex.es>
> http://kolmogorov.unex.es/~onunez
<http://kolmogorov.unex.es/%7Eonunez>
> Tel : +34 663 03 69 09
> Departamento de Matemáticas
> Universidad de Extremadura
>
> ____________________________________
>
>
>
>
>
--
_____---^---_____
Univ. de Extremadura
Dept. Matemáticas.
Despacho B29
Tf: + 34 924 289 300
Ext. 86823
[[alternative HTML version deleted]]