Curso R - Sesión 9
Álvaro Cabana
Medias marginales, pruebas pareadas, contrastes, y todo eso.
Por lo general, uno ajusta un modelo (en el mejor de los casos, nuestro diseño nos dice qué modelo formular), y mira las pruebas de hipótesis sobre los parámetros. Casi siempre, además, nos interesan los “efectos” de las variables, y los de sus interacciones. Generalmente utilizamos para esto último pruebas F o “ANOVAs”.
Pero, como casi siempre tenemos más de una variable predictora, y como mínmo una complejidad de 2x2 en el diseño (dos variables categóricas de dos niveles cada una), nos suele interesar contrastar directamente las diferencias entre pares de niveles. Las (a veces mal llamadas) comparaciones post-hoc. En esta sesión veremos algunas formas de calcular estos contrastes pareados.
Las medias marginales y el paquete emmeans
El paquete emmmeans (antes conocido como) lsmeans permite, entre otras cosas, calcular “medias marginales”, esto es, los promedios de la variable dependiente según diferentes niveles de una o más variables predictoras categóricas.
¿Qué? ¿Que qué cornos quiero decir con esto? Ok, basta de cháchara, pongamos un ejemplo. Eso sí, necesitaremos volver los datos de Valentina Paz et al que usamos en el [ejercicio 5] (http://rpubs.com/almadana/ejercicio5) y en el ejercicio 8. Porque reusar, reciclar, reparar, es un gran mantra. En particular, necesitaremos el data frame largo parts3 que luce así:
load("../.RData")
library(plyr)
str(parts3)## 'data.frame': 336 obs. of 6 variables:
## $ codigo : Factor w/ 112 levels "v001 ",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ sexo : Factor w/ 2 levels "H","M": 2 1 2 2 2 1 2 2 1 2 ...
## $ edad : num 22 27 21 21 25 26 20 21 24 31 ...
## $ Beck : num 2 38 13 8 30 6 13 7 7 21 ...
## $ tipo_oferta: Factor w/ 3 levels "justas","injustas",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ enojo : num 0.667 0 0.333 0.833 0.5 ...En sesiones pasadas vimos cómo usando paquetes como plyr podíamos obtener los promedios para los valores de enojo según, sexo, tipo de oferta, o un para la combinación de los niveles de ambas. Algo así como:
mediasEnojo=ddply(parts3,c("sexo","tipo_oferta"),summarize, mediaEnojo = mean(enojo,na.rm=T))
mediasEnojo## sexo tipo_oferta mediaEnojo
## 1 H justas 0.7333333
## 2 H injustas 4.0375000
## 3 H medias 2.2875000
## 4 M justas 0.3028169
## 5 M injustas 4.0187793
## 6 M medias 1.8098592Ahora, ¿qué pasa si quieren comparar las medias de, digamos, ofertas justas e injustas, para los participantes hombres? Un ejemplo típico de comparación pareada, que cuando se calcula a posteriori (cuando no se “planifica” de antemano, o se realiza mansalva todos contra todos, ver más adelante), se suele llamar post-hoc.
En tal caso, podemos utilizar el paquete emmeans para calcular estimaciones un poco más “prolijas” y menos caseras de las medias. Primero, debemos ajustar un modelo lineal relevante: en este caso, queremos el enojo según sexo y tipo de oferta.
enojo.lm1 = lm(enojo~sexo + tipo_oferta,data=parts3)
summary(enojo.lm1)##
## Call:
## lm(formula = enojo ~ sexo + tipo_oferta, data = parts3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2231 -1.0129 -0.3222 0.9627 4.6537
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6556 0.2182 3.004 0.00287 **
## sexoM -0.3090 0.2080 -1.485 0.13838
## tipo_ofertainjustas 3.5676 0.2446 14.585 < 2e-16 ***
## tipo_ofertamedias 1.5240 0.2446 6.231 1.42e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.822 on 329 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.3968, Adjusted R-squared: 0.3913
## F-statistic: 72.15 on 3 and 329 DF, p-value: < 2.2e-16Primero que nada, calculemos una “grilla de referencia” con ref_grid.
library(emmeans)
#Primer paso, calcular la "grilla de referencia". Esto no es siempre necesario, pero es un buen hábito.
enojo.lm1.rg = ref_grid(enojo.lm1)Luego, calculemos las medias marginales para la variable sexo:
enojo.lm1.emm.sexo = emmeans(enojo.lm1.rg,"sexo")bla
summary(enojo.lm1.emm.sexo)## sexo emmean SE df lower.CL upper.CL
## H 2.352778 0.1663461 329 2.025542 2.680014
## M 2.043818 0.1248571 329 1.798199 2.289438
##
## Results are averaged over the levels of: tipo_oferta
## Confidence level used: 0.95Bien! En promedio, los hombres tienen 2.35 de enojo, y las mujeres 2.04. Eso es lo que informa la columna emmean. La gracia es que ahora sabemos el error estándar de la estimación, los grados de libertad, y el inervalo de confianza (lower.CL y upper.CL son los límites inferior y superior del intervalo).
Ahora, ¿cómo hago para incluir el “desglose” del nivel de enojo según sexo y tipo de oferta?
enojo.lm1.emm.sexo.tipoOferta = emmeans(enojo.lm1.rg,c("sexo","tipo_oferta"))
summary(enojo.lm1.emm.sexo.tipoOferta)## sexo tipo_oferta emmean SE df lower.CL upper.CL
## H justas 0.6555806 0.2182066 329 0.22632451 1.0848367
## M justas 0.3466213 0.1885004 329 -0.02419677 0.7174393
## H injustas 4.2231481 0.2182066 329 3.79389207 4.6524042
## M injustas 3.9141888 0.1885004 329 3.54337080 4.2850069
## H medias 2.1796046 0.2182066 329 1.75034853 2.6088607
## M medias 1.8706453 0.1885004 329 1.49982726 2.2414633
##
## Confidence level used: 0.95Comparemos estos resultados a los que dió ddply…
#según emmeans
summary(enojo.lm1.emm.sexo.tipoOferta)[1,]## sexo tipo_oferta emmean SE df lower.CL upper.CL
## 1 H justas 0.6555806 0.2182066 329 0.2263245 1.084837#según ddply
mediasEnojo[1,]## sexo tipo_oferta mediaEnojo
## 1 H justas 0.7333333Eh… es 0.7333, o es 0.655 ?
Qué está pasando? 
Contrastes
OK, dejemos por un momento esa discrepancia… Sigamos adelante, porque, porqué no? Digamos que nos interesa comparar los niveles de enojo, para cada nivel de oferta, entre hombres y mujeres… Una de las tantas formas de hacerlo, es con la función contrast
contrast(enojo.lm1.emm.sexo.tipoOferta,"pairwise")## contrast estimate SE df t.ratio p.value
## H,justas - M,justas 0.3089593 0.2079912 329 1.485 0.6739
## H,justas - H,injustas -3.5675676 0.2446002 329 -14.585 <.0001
## H,justas - M,injustas -3.2586083 0.3210756 329 -10.149 <.0001
## H,justas - H,medias -1.5240240 0.2446002 329 -6.231 <.0001
## H,justas - M,medias -1.2150647 0.3210756 329 -3.784 0.0025
## M,justas - H,injustas -3.8765269 0.3210756 329 -12.074 <.0001
## M,justas - M,injustas -3.5675676 0.2446002 329 -14.585 <.0001
## M,justas - H,medias -1.8329833 0.3210756 329 -5.709 <.0001
## M,justas - M,medias -1.5240240 0.2446002 329 -6.231 <.0001
## H,injustas - M,injustas 0.3089593 0.2079912 329 1.485 0.6739
## H,injustas - H,medias 2.0435435 0.2446002 329 8.355 <.0001
## H,injustas - M,medias 2.3525029 0.3210756 329 7.327 <.0001
## M,injustas - H,medias 1.7345842 0.3210756 329 5.402 <.0001
## M,injustas - M,medias 2.0435435 0.2446002 329 8.355 <.0001
## H,medias - M,medias 0.3089593 0.2079912 329 1.485 0.6739
##
## P value adjustment: tukey method for comparing a family of 6 estimatesBueno, es como mucho no? Esto compara TODO CONTRA TODO (el argumento “pairwise” hace que se hagan comparaciones de a pares… más sobre esto pronto). La verdad que no me interesa compara H,medias contra M,injustas…
Entonces, si sólo nos interesa hacer comparaciones entre sexos según tipo de oferta, al hacer las medias marginales, podemos especificar en vez de c("sexo","tipo_oferta"), especificar sólo "sexo", y agregar "by=tipo_oferta" como otro argumento a emmeans.
enojo.lm1.emm.sexo.tipoOferta_2 = emmeans(enojo.lm1.rg,"sexo",by="tipo_oferta")
contrast(enojo.lm1.emm.sexo.tipoOferta_2 ,"pairwise")## tipo_oferta = justas:
## contrast estimate SE df t.ratio p.value
## H - M 0.3089593 0.2079912 329 1.485 0.1384
##
## tipo_oferta = injustas:
## contrast estimate SE df t.ratio p.value
## H - M 0.3089593 0.2079912 329 1.485 0.1384
##
## tipo_oferta = medias:
## contrast estimate SE df t.ratio p.value
## H - M 0.3089593 0.2079912 329 1.485 0.1384Ahora sí! Según, cada tipo de oferta, tengo el contraste entre los sexos. Perfecto.
No?
Digamos que quiero el contraste inverso, dentro de cada sexo, comparar niveles de enojo entre las ofertas. Fácil, sólo hay que invertir los argumentos del emmeans:
enojo.lm1.emm.sexo.tipoOferta_2 = emmeans(enojo.lm1.rg,by="sexo","tipo_oferta")
contrast(enojo.lm1.emm.sexo.tipoOferta_2 ,"pairwise")## sexo = H:
## contrast estimate SE df t.ratio p.value
## justas - injustas -3.567568 0.2446002 329 -14.585 <.0001
## justas - medias -1.524024 0.2446002 329 -6.231 <.0001
## injustas - medias 2.043544 0.2446002 329 8.355 <.0001
##
## sexo = M:
## contrast estimate SE df t.ratio p.value
## justas - injustas -3.567568 0.2446002 329 -14.585 <.0001
## justas - medias -1.524024 0.2446002 329 -6.231 <.0001
## injustas - medias 2.043544 0.2446002 329 8.355 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesExcelente!!!
No?
Pero… todos los contrastes son iguales, no importa el sexo… Y en la otra, los contrastes entre sexos eran los mismos, sin importar el nivel del tipo de oferta..
WHAT’S GOING ON?! Qué está pasando?
WHAT’S GOING ON?!
A VER QUIÉN RESUELVE ESTE ENTUERTO!
Comparaciones múltiples
Bueno, digamos que tenemos una chorrera de contrastes, como los de hace un momento (nota: se ha corregido lo que andaba mal hasta hace un momento, y los contrastes ahora son diferentes…. claro que no se los muestro aquí :P):
contrast(enojo.lm1.emm.sexo.tipoOferta,"pairwise")## contrast estimate SE df t.ratio p.value
## H,justas - M,justas 0.43051643 0.3608115 327 1.193 0.8400
## H,justas - H,injustas -3.30416667 0.4080967 327 -8.097 <.0001
## H,justas - M,injustas -3.28544601 0.3608115 327 -9.106 <.0001
## H,justas - H,medias -1.55416667 0.4080967 327 -3.808 0.0023
## H,justas - M,medias -1.07652582 0.3608115 327 -2.984 0.0359
## M,justas - H,injustas -3.73468310 0.3608115 327 -10.351 <.0001
## M,justas - M,injustas -3.71596244 0.3063119 327 -12.131 <.0001
## M,justas - H,medias -1.98468310 0.3608115 327 -5.501 <.0001
## M,justas - M,medias -1.50704225 0.3063119 327 -4.920 <.0001
## H,injustas - M,injustas 0.01872066 0.3608115 327 0.052 1.0000
## H,injustas - H,medias 1.75000000 0.4080967 327 4.288 0.0003
## H,injustas - M,medias 2.22764085 0.3608115 327 6.174 <.0001
## M,injustas - H,medias 1.73127934 0.3608115 327 4.798 <.0001
## M,injustas - M,medias 2.20892019 0.3063119 327 7.211 <.0001
## H,medias - M,medias 0.47764085 0.3608115 327 1.324 0.7717
##
## P value adjustment: tukey method for comparing a family of 6 estimatesY bueno, esto es una chorrera de p-valores… Una cosa frecuente en esta situación, es hacer correcciones por comparaciones múltiples. Estas correcciones pueden ir desde atacar los p-valores con un martillo (Bonferroni), es decir, multiplicar los p-valores por el número de comparaciones realizadas, a estrategias un tanto menos violentas. Hay muchísimas formas de introducir esta corrección, aquí les dejo la forma en que yo suelo hacerlas.
pairs(enojo.lm1.emm.sexo.tipoOferta,adjust="bonferroni")## contrast estimate SE df t.ratio p.value
## H,justas - M,justas 0.43051643 0.3608115 327 1.193 1.0000
## H,justas - H,injustas -3.30416667 0.4080967 327 -8.097 <.0001
## H,justas - M,injustas -3.28544601 0.3608115 327 -9.106 <.0001
## H,justas - H,medias -1.55416667 0.4080967 327 -3.808 0.0025
## H,justas - M,medias -1.07652582 0.3608115 327 -2.984 0.0459
## M,justas - H,injustas -3.73468310 0.3608115 327 -10.351 <.0001
## M,justas - M,injustas -3.71596244 0.3063119 327 -12.131 <.0001
## M,justas - H,medias -1.98468310 0.3608115 327 -5.501 <.0001
## M,justas - M,medias -1.50704225 0.3063119 327 -4.920 <.0001
## H,injustas - M,injustas 0.01872066 0.3608115 327 0.052 1.0000
## H,injustas - H,medias 1.75000000 0.4080967 327 4.288 0.0004
## H,injustas - M,medias 2.22764085 0.3608115 327 6.174 <.0001
## M,injustas - H,medias 1.73127934 0.3608115 327 4.798 <.0001
## M,injustas - M,medias 2.20892019 0.3063119 327 7.211 <.0001
## H,medias - M,medias 0.47764085 0.3608115 327 1.324 1.0000
##
## P value adjustment: bonferroni method for 15 testsEn lugar de bonferroni puede ir mvt (prueba t multivariada) ,holm (Holm-Bonferroni) y fdr (False Discovery Rate, creo que es Benjamini-Hochberg, lean la documentación).
Contraste de contrastes (tumba de campeones)
Hay veces, que más que explorar las diferencias, queremos caracterizarlas, queremos poder saber si dos diferencias, son estadísticamente significativas.
Por ejemplo, la diferencia entre el enojo antes ofertas justas e injustas, es en los hombres de 3.3, y en las mujeres de 3.71. ¿Es esta diferencia significativa?
a=contrast(enojo.lm1.emm.sexo.tipoOferta,"pairwise")
a[c(2,7),]## contrast estimate SE df t.ratio p.value
## H,justas - H,injustas -3.304167 0.4080967 327 -8.097 <.0001
## M,justas - M,injustas -3.715962 0.3063119 327 -12.131 <.0001a[c(11,14),]## contrast estimate SE df t.ratio p.value
## H,injustas - H,medias 1.75000 0.4080967 327 4.288 <.0001
## M,injustas - M,medias 2.20892 0.3063119 327 7.211 <.0001Y la diferencia entre el enojo antes ofertas medias e injustas, es en los hombres de 1.75, y en las mujeres de 2.21. Pero, ¿es esta diferencia significativa?
Para responder a esta pregunta, utilizaremos la función contrast, pero con un nuevo argumento. Y para esto, tendremos que volver un minuto al listado de medias marginales (sí, acá ya están corregidos :P).
summary(enojo.lm1.emm.sexo.tipoOferta)## sexo tipo_oferta emmean SE df lower.CL upper.CL
## H justas 0.7333333 0.2885679 327 0.1656495 1.3010172
## M justas 0.3028169 0.2165952 327 -0.1232790 0.7289128
## H injustas 4.0375000 0.2885679 327 3.4698162 4.6051838
## M injustas 4.0187793 0.2165952 327 3.5926834 4.4448752
## H medias 2.2875000 0.2885679 327 1.7198162 2.8551838
## M medias 1.8098592 0.2165952 327 1.3837633 2.2359551
##
## Confidence level used: 0.95Primero que nada, ¿cuál es el contraste que queremos hacer? Digamos que sería la diferencia entre injustas y justas de los hombres, menos la misma diferencia para las mujeres:
[ (H, injustas) - (H, medias) ] - [ (M, injustas) - (M, medias) ]
o lo que es igual:
(H, injustas) - (H, medias) - (M, injustas) + (M, medias)
Ahora, identifiquemos los números de fila relevantes a nuestra comparación.
fila3 - fila5 - fila4 + fila6
Podemos ahora pasar como parámetro a la función contrast, los coeficientes de esta suma y resta (-1 y 1) para cada una de las filas de la tabla anterior. En nuestro caso, el coeficiente para la filas 1 y 2 es 0, y luego son 1, -1 -1 y 1 para las filas 3, 4, 5 y 6. Así queda entonces:
contrast(enojo.lm1.emm.sexo.tipoOferta,list(hombresVsMujeresEnInjustasMenosMedias=c(0,0,1,-1,-1,1)))## contrast estimate SE df t.ratio
## hombresVsMujeresEnInjustasMenosMedias -0.4589202 0.5102645 327 -0.899
## p.value
## 0.3691Este contraste tiene un estimado de -0.46 (recuerden que la resta entre condiciones daba 1.75 para hombres y 2.21 para mujeres). Eso medio que ya lo sabíamos: pero ahora tenmos una prueba t y un p-valor asociado y podemos decirlo: “no hay diferencias entre los sexos para ese contraste entre injustas y medias”.
¿Quieren probar más conrtrastes? Simplemente agreguen ítems a la lista: