[R-es] Comparaciones múltiples en ANOVA anidadp

Mar Feb 7 14:33:03 CET 2012

José,

dos ultimas precisiones:

1) lmer con REML=TRUE ha de coincidir con el método de mínimos  
cuadrados.

2) Sin embargo, la salida de la tabla ANOVA que proporciona la  
función anova() aplicada a lmer no da las MS estándares,
sino más bien, una MS "estandarizadas". En todo caso, las F han de  
coincidir (o ser muy próxima).

Para verlo, una pequeña simulación similar a la tuya, con efectos  
individuales y efectos fijos de ciertos tratamientos:

mu <- 100 #media global
error <- rnorm(45, 0, 1) #error de medidas con desviación típica  
igual a 1
T.effects <- c(50,-50,0,0,0) #efectos de los tratamientos; sólo los  
dos primero tratamientos tienen efectos
I.effects <- rnorm(15,0,3) #efectos individuales aleatorio desviación  
típica igual a 3
VR <- mu + rep(T.effects,each=9) + rep(I.effects,each=3) + error  
#simulación de las medidas
Tratamiento <- gl(5, 9) #niveles de los tratamientos
Individuo <- gl(15, 3) #etiqueta del individuo
Medidas <-data.frame(VR=VR,Tratamiento=Tratamiento,Individuo=Individuo)

La salida de lmer me da

 > options(contrasts=c(factor="contr.sum", ordered="contr.poly"))
 > fit.lmer <- lmer(VR ~ Tratamiento+(1|Individuo:Tratamiento),  
data=Medidas,REML=FALSE)
 > summary(fit.lmer)

Linear mixed model fit by REML
Formula: VR ~ Tratamiento + (1 | Individuo:Tratamiento)
    Data: Medidas
    AIC   BIC logLik deviance REMLdev
  172.3 184.9 -79.13    169.6   158.3
Random effects:
  Groups                Name        Variance Std.Dev.
  Individuo:Tratamiento (Intercept) 9.0434   3.00723
  Residual                          0.9108   0.95436
Number of obs: 45, groups: Individuo:Tratamiento, 15

Fixed effects:
              Estimate Std. Error t value
(Intercept)   98.9494     0.7894  125.35
Tratamiento1  49.6228     1.5787   31.43
Tratamiento2 -51.6139     1.5787  -32.69
Tratamiento3   1.1821     1.5787    0.75
Tratamiento4  -0.1764     1.5787   -0.11

Correlation (omitida)

Como puedes ver, la varianza de los efectos aleatorios y los efectos  
fijos están relativamente bien estimados.
La residual teóricamente igual a 1 es aquí estimada por 0.91.

SI aplico la función anova() a este modelo obtengo:

 > anova(fit.lmer)
Analysis of Variance Table
             Df Sum Sq Mean Sq F value
Tratamiento  4 1496.9  374.23  410.88

Sin embargo, no se trata aquí de MS(Tratamiento) sino del ratio de MS 
(Tratamiento)/ (MS( Individuo:Tratamiento)/MS(Residual)),
de manera que, la F correspondiente se obtiene dividiendo por la MS 
(Residual).
Se puede confirmar (puede ser pedagógico) con la salida de aov:

 > fit.aov <- aov(VR ~ Tratamiento + Error 
(Individuo),data=Medidas);summary(fit.aov)

Error: Individuo
             Df Sum Sq Mean Sq F value    Pr(>F)
Tratamiento  4  46159   11540  411.53 4.792e-11 ***
Residuals   10    280      28
---
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 27.324 0.91081

Tenemos que 11540/(28/.91) = 375 que coincide aproximadamente con la  
MS que nos da 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 07/02/2012, a las 12:09, José Trujillo Carmona escribió:

> 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 en iberstat.es <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto:R-help-es en r-project.org>
>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>
>>
>>
>> -- 
>> ____________________________________
>>
>> Olivier G. Nuñez
>> Email: onunez en unex.es <mailto:onunez en 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]]
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es