Capítulo 9 Después del ANOVA

En este capítulo vamos a revisar que podemos hacer después de un ANOVA.

9.1 ¿Qué hacemos después de un ANOVA?

Revisa este video (12’) y trata de responder:

¿Que hago después de una prueba omnibus?
¿De qué depende el tipo de análisis que se hace después del ANOVA?
¿Qué significa que desglocemos la varianza explicada en los distintos contrastes en el análisis de contrastes planificados?
¿En un experimento dónde analizo el efecto de dos dosis de droga comparada con un placebo, que contrastes planificados se podrían hacer y que indican el intercepto y los betas de este modelo?
¿Por que es importante que el contraste que hagamos sea un contraste ortogonal?

9.2 Comparaciones planificadas

Revisa este video (18’) y trata de responder:

¿Para que sirven los contrastes que evaluan una tendencia líneal o cuadrática?

9.3 Cálculos con R

Ahora vamos a hacer el análisis. Primero seteamos nuestro directorio de trabajo y cargamos las librerías que necesitemos (si no las tienes instaladas debes instalarlas).

setwd(Sys.getenv("BOOKDOWN_PROJECT_PATH"))

library(ggplot2)
library(Hmisc)
library(Rmisc)
library(effsize)
library(pastecs)
library(reshape2)
library(car)
library(effsize)

Primero, vamos a revisar el código para realizar los contrastes planificados. Recuerda que se utiliza un dummy coding en un modelo de regresión para que los valores b representen las diferencias entre los promedios de los contrastes que se quieren evaluar.

Importamos el set de datos. Aprovechamos de agregar unas etiquetas que nos sirvan a visualizar las condiciones experimentales. Y le damos una mirada.

contrastData <- read.csv("data/inteligencia_c_Contrast_Planif.csv", header = TRUE)
contrastData$Dosis <- factor(contrastData$Dosis, 
                             levels = c(1:3), labels = c("Placebo", "Baja", "Alta"))
contrastData

##      Dosis Inteligencia Compars1 Compars2
## 1  Placebo            3       -2        0
## 2  Placebo            2       -2        0
## 3  Placebo            1       -2        0
## 4  Placebo            1       -2        0
## 5  Placebo            4       -2        0
## 6     Baja            5        1       -1
## 7     Baja            2        1       -1
## 8     Baja            4        1       -1
## 9     Baja            2        1       -1
## 10    Baja            3        1       -1
## 11    Alta            7        1        1
## 12    Alta            4        1        1
## 13    Alta            5        1        1
## 14    Alta            3        1        1
## 15    Alta            6        1        1

Usemos las columnas que corresponden a las variables dummy. Compars1 refleja la diferencia para los niveles de inteligencia entre el placebo y cualquier dosis de droga (en realidad el promedio de ambas). Compars2 refleja la diferencia para los niveles de inteligencia entre la dosis baja y la dosis alta.

m_contrastedAOV <- lm(Inteligencia ~ Compars1 + Compars2, data = contrastData)
summary(m_contrastedAOV)

## 
## Call:
## lm(formula = Inteligencia ~ Compars1 + Compars2, data = contrastData)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.0   -1.2   -0.2    0.9    2.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
## Compars1      0.6333     0.2560   2.474   0.0293 *  
## Compars2      0.9000     0.4435   2.029   0.0652 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.402 on 12 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3704 
## F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

¿Qué nos dice este modelo? El F es el mismo. Lo que cambia son los coeficientes de regresión. Primero. El intercepto es el gran promedio (3.467). Verifiquemos.

# promedio a través de todas las condiciones
mean(contrastData$Inteligencia)

## [1] 3.466667

# promedio a través de todas las condiciones experimentales
varX <- contrastData[contrastData$Dosis != "Placebo",]
mean(varX$Inteligencia)

## [1] 4.1

El contraste 1 (o b1) refleja la diferencia en los niveles de inteligencia entre la aplicación de la droga y el placebo. El coeficiente de regresión (0.63) para el contraste 1 es 1/3 de esta diferencia (1.9/3 = 0.633). El FWE se controla haciendo que el coeficiente de regresión sea igual a la diferencia real dividida por el número de grupos en el contraste (en este caso 3).

El contraste 2 (o b2) refleja la diferencia en los niveles de inteligencia entre la aplicación de una dosis alta droga y una dosis baja de droga. El coeficiente de regresión (0.9) para el contraste 1 es 1/2 de esta diferencia (1.8/2 = 0.9). El FWE se controla haciendo que el coeficiente de regresión sea igual a la diferencia real dividida por el número de grupos en el contraste (en este caso 3).

Uno podría pensar que siempre hay que crear las dos columnas del dummy coding en los datos. En realidad eso se puede hacer en R.

Importemos los datos originales.

inteligData2 <- read.csv("data/inteligencia.csv", header = TRUE)
inteligData2$Dosis <- factor(inteligData2$Dosis, levels = c(1:3),
                     labels = c("Placebo", "Baja", "Alta"))

Si revisas los valores de dosis encontrarás que tienen los nombres correspondientes a las condiciones experimentales, y al final veras que tiene 3 niveles: Placebo, Baja y Alta

inteligData2$Dosis

##  [1] Placebo Placebo Placebo Placebo Placebo Baja    Baja    Baja    Baja    Baja    Alta    Alta    Alta    Alta    Alta   
## Levels: Placebo Baja Alta

Por defecto R toma el primer grupo base y crea dos variables dummy que se comparan con ese grupo base. Puedes mirar estos contrastes con la función contrasts.

contrasts(inteligData2$Dosis)

##         Baja Alta
## Placebo    0    0
## Baja       1    0
## Alta       0    1

Pero, tu puedes crear manualmente contrastes específicos. Por ejemplo, con la función contr.treatment puedes crear 2 variables dummy (a partir de 3 condiciones exprimentales) dónde decidas que el grupo base es el segundo.

contr.treatment(3, base = 2)

##   1 3
## 1 1 0
## 2 0 0
## 3 0 1

Luego para que estos contraste surgan efecto lo inyectas como contrastes

contrasts(inteligData2$Dosis) <- contr.treatment(3, base = 2)

Luego verificas como quedaron los contrastes.

inteligData2$Dosis

##  [1] Placebo Placebo Placebo Placebo Placebo Baja    Baja    Baja    Baja    Baja    Alta    Alta    Alta    Alta    Alta   
## attr(,"contrasts")
##         1 3
## Placebo 1 0
## Baja    0 0
## Alta    0 1
## Levels: Placebo Baja Alta

Y finalmente realizas el ANOVA nuevamente.

m_AOV3 <- aov(Inteligencia ~ Dosis, data = inteligData2)
summary.lm(m_AOV3)

## 
## Call:
## aov(formula = Inteligencia ~ Dosis, data = inteligData2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.0   -1.2   -0.2    0.9    2.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.2000     0.6272   5.102 0.000261 ***
## Dosis1       -1.0000     0.8869  -1.127 0.281584    
## Dosis3        1.8000     0.8869   2.029 0.065192 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.402 on 12 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3704 
## F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

Fíjate los betas son diferentes.

Ahora, en R hay ciertos contrastes que ya están preparados. Uno de ellos es el contraste Helmert. Trata de entender que contrastes se hacen en este caso. Igual que antes creas el contraste lo inyectas en la base de datos y verificas.

contr.helmert(3)

##   [,1] [,2]
## 1   -1   -1
## 2    1   -1
## 3    0    2

contrasts(inteligData2$Dosis) <- contr.helmert(3)
inteligData2$Dosis

##  [1] Placebo Placebo Placebo Placebo Placebo Baja    Baja    Baja    Baja    Baja    Alta    Alta    Alta    Alta    Alta   
## attr(,"contrasts")
##         [,1] [,2]
## Placebo   -1   -1
## Baja       1   -1
## Alta       0    2
## Levels: Placebo Baja Alta

Luego calculas el ANOVA nuevamente.

m_AOV4 <- aov(Inteligencia ~ Dosis, data = inteligData2)
summary.lm(m_AOV4)

## 
## Call:
## aov(formula = Inteligencia ~ Dosis, data = inteligData2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.0   -1.2   -0.2    0.9    2.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
## Dosis1        0.5000     0.4435   1.127   0.2816    
## Dosis2        0.7667     0.2560   2.994   0.0112 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.402 on 12 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3704 
## F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

Más importante que esto es que puedas crear tu propios contrastes. Para ello, te recomiendo darles nombres que te permitan entender esos contrastes. Por ejemplo, siguiendo nuestro ejemplo anterior podríamos tener dos hipótesis específicas. Primero, comparamos el efecto del tratamiento de la droga versus el grupo placebo. Para ello concatenamos -2, 1 y 1, que corresponden respectivamente a Placebo, Low y High.

C_vs_Treatment <- c(-2,1,1)

Segundo, comparamos el efecto del dosis alta versus la dsis baja. Para ello concatenamos 0, -1 y 1, que corresponden respectivamente a Placebo, Low y High.

Low_vs_High <- c(0,-1,1)

Luego ponemos los dos contrastes en la base de datos y calculamos un nuevo ANOVA.

contrasts(inteligData2$Dosis) <- cbind(C_vs_Treatment, Low_vs_High)
m_AOV4 <- aov(Inteligencia ~ Dosis, data = inteligData2)
summary.lm(m_AOV4)

## 
## Call:
## aov(formula = Inteligencia ~ Dosis, data = inteligData2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.0   -1.2   -0.2    0.9    2.0 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.4667     0.3621   9.574 5.72e-07 ***
## DosisC_vs_Treatment   0.6333     0.2560   2.474   0.0293 *  
## DosisLow_vs_High      0.9000     0.4435   2.029   0.0652 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.402 on 12 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3704 
## F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

Fíjate que obtenemos los mismos resultados que cuando utilizamos la base de datos dónde orginalmente se había definido las variables dummy.

Por último, algo que también podriamos hacer es un análisis de tendencia líneal. Si pensamos en los efectos de las drogas (vean el gráfico arriba) esperaríamos que a medida que aumenta la dosis aumentará su efecto de manera lineal. Para ello podemos crear un contraste que de cuenta de este efecto lineal. Para ello usamos la función contr.poly y lo inyectamos en la base de datos.

contr.poly(3)

##                 .L         .Q
## [1,] -7.071068e-01  0.4082483
## [2,] -7.850462e-17 -0.8164966
## [3,]  7.071068e-01  0.4082483

contrasts(inteligData2$Dosis) <- contr.poly(3)
contrasts(inteligData2$Dosis)

##                    .L         .Q
## Placebo -7.071068e-01  0.4082483
## Baja    -7.850462e-17 -0.8164966
## Alta     7.071068e-01  0.4082483

Finalmente, calculamos el ANOVA nuevamente

trend_M <- aov(Inteligencia ~ Dosis, data = inteligData2)
summary.lm(trend_M)

## 
## Call:
## aov(formula = Inteligencia ~ Dosis, data = inteligData2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.0   -1.2   -0.2    0.9    2.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
## Dosis.L       1.9799     0.6272   3.157  0.00827 ** 
## Dosis.Q       0.3266     0.6272   0.521  0.61201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.402 on 12 degrees of freedom
## Multiple R-squared:  0.4604, Adjusted R-squared:  0.3704 
## F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469

La fila “L” indica un efecto lineal y la fila “Q” indica un efecto cuadrático. ¿Qué ves?

9.4 Pruebas post-hoc

Revisa este video (15’) y trata de responder:

¿Qué tipos de pruebas post-hoc existen?
¿Cuál es el riesgo que corremos si usamos una prueba post-hoc muy exigente (o conservadora)?
¿Cuál es el riesgo que corremos si usamos una prueba post-hoc poco exigente (o conservadora)?

9.5 Cálculos con R

Ahora vamos a revisar el código para realizar las pruebas post-hoc.

Usemos el set de datos originales. Aprovechamos de agregar unas etiquetas que nos sirvan a visualizar las condiciones experimentales. Y le damos una mirada.

dat1 <- read.csv("data/inteligencia.csv", header = TRUE)
dat1$Dosis <- factor(dat1$Dosis, 
                            levels = c(1:3), labels = c("Placebo", "Baja", "Alta"))

En las pruebas post-hoc comparamos todo con todo. Para ello usamos la función pairwise.t.test. Pero, además podemos aplicar una corrección (en el parámetro p.adjust.method) para lidiar con el error de tipo 1.

Veamos que pasa cuando no agregamos ninguna corrección.

pairwise.t.test(dat1$Inteligencia, dat1$Dosis, p.adjust.method = "none")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dat1$Inteligencia and dat1$Dosis 
## 
##      Placebo Baja  
## Baja 0.2816  -     
## Alta 0.0083  0.0652
## 
## P value adjustment method: none

Veamos que pasa cuando usamos una corrección de Bonferroni.

pairwise.t.test(dat1$Inteligencia, dat1$Dosis, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dat1$Inteligencia and dat1$Dosis 
## 
##      Placebo Baja 
## Baja 0.845   -    
## Alta 0.025   0.196
## 
## P value adjustment method: bonferroni

Veamos que pasa cuando usamos una corrección de Benjamini-Hochberg.

pairwise.t.test(dat1$Inteligencia, dat1$Dosis, p.adjust.method = "BH")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dat1$Inteligencia and dat1$Dosis 
## 
##      Placebo Baja 
## Baja 0.282   -    
## Alta 0.025   0.098
## 
## P value adjustment method: BH