Contrastes planeados

Un contraste planeado es un contraste que el experimentador a escogido previo a ver los resultados de un estudio. Volviendo a nuestro ejemplo anterior, usamos dos contrastes: uno comparando el efecto de la música envasada y el efecto de la música envivo, y otro comparando el efecto de ambos tratamientos con una condición sin música.

aov1 = aov(satisfaccion ~ condicion, 
           data = df)
summary(aov1)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## condicion    2     48  24.000   11.25 0.00104 **
## Residuals   15     32   2.133                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Requiere el paquete emmeans
#install.packages("emmeans")
require(emmeans)

## Loading required package: emmeans

## Warning: package 'emmeans' was built under R version 3.5.3

aov1.emm = emmeans(aov1, specs =~condicion)

# Agreguemos un contraste para la diferencia entre la 
# condición sin música y las otras dos condiciones 
# combinadas
Contrasts.2 = list(Musica_en_vivo_vivo_vs_Musica_Envasada  = c(-1,1,0),
                 Todos_vs_Sin_Musica     = c(1/2,1/2,-1))

contrast(aov1.emm, Contrasts.2)

##  contrast                               estimate    SE df t.ratio p.value
##  Musica_en_vivo_vivo_vs_Musica_Envasada        2 0.843 15  2.372  0.0315 
##  Todos_vs_Sin_Musica                          -3 0.730 15 -4.108  0.0009

Para determinar el valor p ajustado por el número de contrastes simplemente tenemos que dividir el nivel de \(\alpha\) especificado por el número de contrastes aplicados, y comparar el valor p ajustado con el valor p para cada contraste. En este caso, para mantener el error Tipo I en .05, dividimos .05/2=.025. En este caso, sólo el segundo contraste (Todos_vs_Sin música) es estadísticamente significativo. Este ajuste se llama ajuste de Bonferroni11 En el caso de contrastes no ortogonales también se puede hacer el mismo ajuste, sólo que es más conservador que cuando los contrastes son ortogonales (Maxwell, Delaney, and Kelley ⊕2018Maxwell, S. E., H. D. Delaney, and K. Kelley. 2018. Designing Experiments and Analyzing Data: A Model Comparison Perspective. 3rd ed. New York, NY: Routledge., 227)

El ajuste de Bonferroni también se puede solicitar directamente a través de la función contrast especificando la opción “bonferroni” en ajustes (adjust). En este caso, el ajuste se ve reflejado directamente en el valor p ajustado.22 Nótese que también se reportan los intervalos de confianza para los contrastes. Hasta ahora no he podido verificar cómo emmeans produce estos límites pero puede que no sean correctos si asume que C=3 en vez de C=2 como en nuestro caso

contrast(aov1.emm,Contrasts.2,adjust = "bonferroni")

##  contrast                               estimate    SE df t.ratio p.value
##  Musica_en_vivo_vivo_vs_Musica_Envasada        2 0.843 15  2.372  0.0630 
##  Todos_vs_Sin_Musica                          -3 0.730 15 -4.108  0.0019 
## 
## P value adjustment: bonferroni method for 2 tests

Una alternativa al ajuste de Bonferroni es el ajuste de Sidak que asegura que la tasa de error Tipo será igual o menor que .05 incluso cuando los contrastes son no ortogonales.

contrast(aov1.emm,Contrasts.2,adjust = "sidak")

##  contrast                               estimate    SE df t.ratio p.value
##  Musica_en_vivo_vivo_vs_Musica_Envasada        2 0.843 15  2.372  0.0620 
##  Todos_vs_Sin_Musica                          -3 0.730 15 -4.108  0.0019 
## 
## P value adjustment: sidak method for 2 tests

Como se puede observar, aunque la conclusión sustantiva es la misma, los valores p asociados al ajuste de Sidak con algo más pequeños que los obtenidos vía Bonferroni. Es decir, el ajuste de Sidak es más poderoso que el ajuste de Bonferroni.

Comparaciones post-hoc

Cuando el investigador tiene la intención de hacer comparaciones habiendo observado los resultados, es crítico controlar por la tasa de error Tipo I por cada una de las comparaciones posibles.Estos de llaman test post-hoc en general. En el caso de la comparación entre pares de medias,podemos usar el ajuste de Bonferroni para todos los pares de medias o \(a(a-1)/2\). En el caso de 3 grupos (a = 3) el valor ajustado para p sería

\[C=a(a-1)/2=3(3-1)/2=3\]

##  contrast                               estimate    SE df t.ratio p.value
##  Musica_en_vivo_vivo_vs_Musica_Envasada        2 0.843 15 2.372   0.0315 
##  Musica_en_vivo_vivo_vs_Sin_Musica             4 0.843 15 4.743   0.0003 
##  Musica_Envasada_vs_Sin_Musica                 2 0.843 15 2.372   0.0315

Comparando los valores p observados con los valores p ajustados (0.5/3 = 0.017) sólo la comparación entre los grupos extremos es estadísticamente significativa. Sin embargo, el ajuste de Bonferroni es menos poderoso el el ajuste de Tukey (llamado Tukey HSD). Podemos hacer las comparaciones entre medias simplemente invocando el ajuste de Tukey en el paquete emmeans

contrast(aov1.emm, "tukey")

##  contrast                         estimate    SE df t.ratio p.value
##  Música en vivo - Música envasada       -2 0.843 15 -2.372  0.0761 
##  Música en vivo - Sin música            -4 0.843 15 -4.743  0.0007 
##  Música envasada - Sin música           -2 0.843 15 -2.372  0.0761 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

O también mediante la función TukeyHSD() del paquete base pero aplicado al objeto aov1 no a aov1.emm que es del tipo emmGrid

TukeyHSD(aov1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = satisfaccion ~ condicion, data = df)
## 
## $condicion
##                                diff        lwr      upr     p adj
## Música envasada-Música en vivo    2 -0.1903792 4.190379 0.0760902
## Sin música-Música en vivo         4  1.8096208 6.190379 0.0007194
## Sin música-Música envasada        2 -0.1903792 4.190379 0.0760902

Una última alternativa que veremos por ahora es el ajuste de Dunnett que especial para comparar un grupo control vs varios grupos experimentales. Por defecto, la función contrast toma como grupo referencia el nivel 1 del factor, que en nuestro cso es el grupo con música en vivo

contrast(aov1.emm, "dunnett")

##  contrast                         estimate    SE df t.ratio p.value
##  Música envasada - Música en vivo        2 0.843 15 2.372   0.0588 
##  Sin música - Música en vivo             4 0.843 15 4.743   0.0005 
## 
## P value adjustment: dunnettx method for 2 tests

Para nuestro caso nos interesa usar como grupo control el grupo sin música por lo que el grupo de referencia en este caso corresponde al nivel 3 lo que especificamos usando ref=3 en la función anterior.

contrast(aov1.emm, "dunnett", ref=3)

##  contrast                     estimate    SE df t.ratio p.value
##  Música en vivo - Sin música        -4 0.843 15 -4.743  0.0005 
##  Música envasada - Sin música       -2 0.843 15 -2.372  0.0588 
## 
## P value adjustment: dunnettx method for 2 tests

Para nuestra “desgracia” el valor p de sólo es significativo para una de estas comparaciones.