ANOVA

ANOVA como un caso especial de regresión

Gonzalo J. Muñoz

2019-04-03

No es inusual escuchar que ANOVA es un caso especial del análisis de regresión. Aunque esto es fundamentalmente cierto existen razones fundamentadas para continuar usando ANOVA como modelo de análisis (e.g., Maxwell, Delaney, & Kelley (2018). Por ello, me remitiré a ejemplificar el caso de ANOVA de una vía mediante regresión, pero no volveremos sobre este problema más adelante.

El análisis de regresión permite estudiar el efecto de una o más variables sobre una variable dependiente continua. Aunque típicamente los predictores de un modelo de regresión son variables continuas también es posible especificar un modelo donde los predictores sean variables categóricas o discretas. Luego, se dice que ANOVA es un caso especial de regresión ya que todos los análisis que se pueden hacer con ANOVA se pueden hacer con regresión, pero no viceversa (e.g., ANOVA no puede usar predictores continuos).

# Agreguemos una variables (dummy code) para representar la pertenencia de cada 
# sujeto (s1...s18) a las distintas condiciones experimentales.
df$X1 <- c(rep(0,6),rep(-1,6),rep(0,6))
df$X2 <- c(rep(0,6),rep(0,6),rep(1,6))

Nivel de satisfacción con el restaurant para los distintos tipos de música.

Como se puede observar, agregamos 2 predictores (X1 y X2) a la data df. ¿Pero si tenemos 3 condiciones experimentales por qué se agregan sólo dos predictores?.11 En general, al incorporar variables dummy basta con agregar (a - 1) variables dummy parastimar el efecto de un factor La categoría que no recibe ningún código–que en este caso corresponde a la condición Música en vivo–queda implicada como categoría de referencia (esto quedará claro a continuación). Si hacemos una regresión lineal22 En este link hay una guía rápida para correr modelos de regresión https://www.statmethods.net/stats/regression.html usando X1 y X2 como predictores tenemos que:

# Ejemplo Regresión Lineal Múltiple
fit.reg <- lm(satisfaccion ~ X1 + X2, data =df)
summary(fit.reg)

## 
## Call:
## lm(formula = satisfaccion ~ X1 + X2, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##     -2     -1      0      1      2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.0000     0.5963   3.354 0.004349 ** 
## X1           -2.0000     0.8433  -2.372 0.031518 *  
## X2            4.0000     0.8433   4.743 0.000261 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.461 on 15 degrees of freedom
## Multiple R-squared:    0.6,  Adjusted R-squared:  0.5467 
## F-statistic: 11.25 on 2 and 15 DF,  p-value: 0.001036

En la tabla de coeficientes (Coefficients) el valor para el intercepto corresponde a la condición Música en vivo. De hecho, 2 es la media para esta condición. Luego, los coeficientes para X1 y X2 representan los que hay que “agregarle” o “quitarle” al intercepto para obtener la media para el grupo Música envasada y para el grupo Sin música.

Específicamente, usando estos contrastes, se puede calcular la media para cada grupo de la siguiente manera: \[\hat{Y}_{ij}=b_0+b_1X_1+b_2X_2\] \[\hat{Y}_{ij}=2+(-2)X_1+(4)X_2\] Para el grupo con Música en vivo el valor predicho para cada individuo (o la media del grupo) sería: \[\hat{Y}_{mev}=2+(-2)0+(4)0=2\] Para el grupo con Música envasada el valor predicho para cada individuo (o la media del grupo) sería: \[\hat{Y}_{mev}=2+(-2)(-1)+(4)0=4\] Y para el grupo con Sin música el valor predicho para cada individuo (o la media del grupo) sería: \[\hat{Y}_{mev}=2+(-2)(0)+(4)1=6\]

require(psych)

## Loading required package: psych

describeBy(df[,"satisfaccion"],condicion)

## 
##  Descriptive statistics by group 
## group: Música en vivo
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 6    2 1.41      2       2 1.48   0   4     4    0    -1.58 0.58
## -------------------------------------------------------- 
## group: Música envasada
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 6    4 1.41      4       4 1.48   2   6     4    0    -1.58 0.58
## -------------------------------------------------------- 
## group: Sin música
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 6    6 1.55      6       6 1.48   4   8     4    0    -1.96 0.63

Es interesante comparar este resultado con el output de ANOVA que vimos anteriormente.

fit.anova<-aov(satisfaccion~condicion,data=df)
summary(fit.anova)

##             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

Podemos volver a realizar el mismo análisis pero con un esquema de códigos diferentes. Por ejemplo, tomemos ahora el grupo Sin música como categoría de referencia modificando levemente el esquema de codificación:

# Agreguemos una variables (dummy code) para representar la pertenencia de cada 
# sujeto (s1...s18) a las distintas condiciones experimentales ()
df$X1 <- c(rep(0,6),rep(-1,6),rep(0,6))
df$X2 <- c(rep(1,6),rep(0,6),rep(0,6))

Nivel de satisfacción con el restaurant para los distintos tipos de música.

Y volvemos a realizar el análisis de regresión tal como antes:

# Ejemplo Regresión Lineal Múltiple
fit.reg <- lm(satisfaccion ~ X1 + X2, data =df)
summary(fit.reg)

## 
## Call:
## lm(formula = satisfaccion ~ X1 + X2, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##     -2     -1      0      1      2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0000     0.5963  10.062  4.6e-08 ***
## X1            2.0000     0.8433   2.372 0.031518 *  
## X2           -4.0000     0.8433  -4.743 0.000261 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.461 on 15 degrees of freedom
## Multiple R-squared:    0.6,  Adjusted R-squared:  0.5467 
## F-statistic: 11.25 on 2 and 15 DF,  p-value: 0.001036

Por último, es importante considerar que para evaluar este modelo el efecto de las variables dummy (X1 y X2) deben ser evaluado incorporando ambas variables concurrentemente.

Llamémosle fit1 va al modelo restringido y fit2 al modelo saturado. En R podemos comparar dos modelos usando la función anova como se muestra en el ejemplo a continuación.

# Nótese que en el modelo fit1 uso "~ 1" para indicar que este parámetro es una # constante, que usualmente es la media poblacional y que corresponde a la gran media. Sólo estimamos en este modelo la varianza residual

fit1 <- lm(satisfaccion ~ 1, data =df)
summary(fit1)

## 
## Call:
## lm(formula = satisfaccion ~ 1, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -4.00  -1.75   0.00   1.00   4.00 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.0000     0.5113   7.823 4.94e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.169 on 17 degrees of freedom

# A continuación agregamos nuestro factor como un set
fit2 <- lm(satisfaccion ~ X1 + X2, data =df)
summary(fit2)

## 
## Call:
## lm(formula = satisfaccion ~ X1 + X2, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##     -2     -1      0      1      2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0000     0.5963  10.062  4.6e-08 ***
## X1            2.0000     0.8433   2.372 0.031518 *  
## X2           -4.0000     0.8433  -4.743 0.000261 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.461 on 15 degrees of freedom
## Multiple R-squared:    0.6,  Adjusted R-squared:  0.5467 
## F-statistic: 11.25 on 2 and 15 DF,  p-value: 0.001036

# Y comparamos el ajuste de ambos modelos
anova(fit1, fit2)

## Analysis of Variance Table
## 
## Model 1: satisfaccion ~ 1
## Model 2: satisfaccion ~ X1 + X2
##   Res.Df RSS Df Sum of Sq     F   Pr(>F)   
## 1     17  80                               
## 2     15  32  2        48 11.25 0.001036 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El valor de RSS (residual sum of squares) en el modelo restringido es 80, y en el modelo saturado es 32. Los grados de libertad del modelo restringido es N - 1 y en el modelo saturado N - a