ANOVA DE VARIABLES DEPENDIENTES DE MEDIAS REPETIDAS

1.- INTRODUCCIÓN

Cuando las variables a comparar son mediciones distintas; pero sobre los mismos sujetos, no se cumple la condición de independencia, por lo que se requiere un ANOVA específico que realice comparaciones considerando que los datos son pareados (de forma similar como se hace en los t-test pareados pero para comparar más de dos grupos).

Es la única condición para poder aplicar este tipo de análisis. La esfericidad implica que la varianza de las diferencias entre todos los pares de variables a comparar sea igual. Si el ANOVA se realiza para comparar la media de una variable cuantitativa entre tres niveles.

La función aov() permite realizar ANOVA de medidas repetidas; pero no comprueba la esfericidad ni permite incorporar las correcciones. La función Anova() del paquete car es mas completa, devuelve, junto con los resultados del anova, el test de esfericidad de Mauchly y los valores que se obtiene con las correcciones de Greenhouse-Geisser y la de Huynh- Feldt. La funciön Anova() requiere que los datos esten almacenados en una matriz en formato de “tabla ancha”.

2.- EJERCICIO

Supóngase un estudio en el que se quiere comprobar si el precio de la compra varía entre 4 cadenas de supermercado distintas. Para ellos se selecciona una serie de elementos de la compra cotidiana y se registra su valor en cada uno de los supermercados ¿Existen evidencias de que el precio medio de la compra es diferente dependiendo del supermercado?

elemento <- c("lettuce", "potatoes", "milk", "eggs", "bread", "cereal", "ground.beef",
              "tomato.soup", "laundry.detergent", "aspirin")
tienda_A <- c(1.755, 2.655, 2.235, 0.975, 2.370, 4.695, 3.135, 0.930, 8.235, 6.690)
tienda_B <- c(1.78, 1.98, 1.69, 0.99, 1.70, 3.15, 1.88, 0.65, 5.99, 4.84)
tienda_C <- c(1.29, 1.99, 1.79, 0.69, 1.89, 2.99, 2.09, 0.65, 5.99, 4.99)
tienda_D <- c(1.29, 1.99, 1.59, 1.09, 1.89, 3.09, 2.49, 0.69, 6.99, 5.15)

datos <- data.frame(elemento, tienda_A, tienda_B, tienda_C, tienda_D)
datos
##             elemento tienda_A tienda_B tienda_C tienda_D
## 1            lettuce    1.755     1.78     1.29     1.29
## 2           potatoes    2.655     1.98     1.99     1.99
## 3               milk    2.235     1.69     1.79     1.59
## 4               eggs    0.975     0.99     0.69     1.09
## 5              bread    2.370     1.70     1.89     1.89
## 6             cereal    4.695     3.15     2.99     3.09
## 7        ground.beef    3.135     1.88     2.09     2.49
## 8        tomato.soup    0.930     0.65     0.65     0.69
## 9  laundry.detergent    8.235     5.99     5.99     6.99
## 10           aspirin    6.690     4.84     4.99     5.15

Para poder visualizar los datos con ggplot2 y realizar la exploración inicial se pasa de formato de “tabla ancha” a “tabla larga” con una columna por variable.

require(tidyr)
## Loading required package: tidyr
## Warning: package 'tidyr' was built under R version 4.2.3
datos_tabla_larga <- gather(data = datos, key = "tienda", value = "precio", 2:5)
head(datos_tabla_larga, 5)
##   elemento   tienda precio
## 1  lettuce tienda_A  1.755
## 2 potatoes tienda_A  2.655
## 3     milk tienda_A  2.235
## 4     eggs tienda_A  0.975
## 5    bread tienda_A  2.370

Es recomendable calcular el precio total de cada uno de los grupos, así como una representación gráfica para hacerse una idea de cuáles podrían diferir significativamente.

with(data = datos_tabla_larga,expr = tapply(precio, tienda, mean))
## tienda_A tienda_B tienda_C tienda_D 
##   3.3675   2.4650   2.4360   2.6260
with(data = datos_tabla_larga,expr = tapply(precio, tienda, var))
## tienda_A tienda_B tienda_C tienda_D 
## 5.955412 2.915317 3.116271 3.951182
with(data = datos_tabla_larga,expr = tapply(precio, tienda, sd))
## tienda_A tienda_B tienda_C tienda_D 
## 2.440371 1.707430 1.765296 1.987758
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
ggplot(data=datos_tabla_larga, aes (x=tienda, y=precio, colour = tienda))+ geom_boxplot()+theme_bw()+theme(legend.position = "bottom")

Se aprecia que la media de las tiendas es similar en todos los casos. La dispersión de la tienda A es sustancialmente superior a la del resto. No obstante, la dispersión de B, C y D es similar. *Se detectan valores atípicos en A, B, C y D.

Podemos concluir que son datos dependientes.

2.2.- ANOVA DE DATOS PAREADOS

anova_pareado <- aov(precio ~ tienda + Error (elemento/tienda), data = datos_tabla_larga)
summary(anova_pareado)
## 
## Error: elemento
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  139.5    15.5               
## 
## Error: elemento:tienda
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## tienda     3  5.737  1.9124   13.03 1.9e-05 ***
## Residuals 27  3.964  0.1468                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tienda: El estadístico F asceinde a 13.03 con un p valor de 0.000 < 0.05. Tenemos evidencia empírica suficiente para indicar que la media de la cesta de la compra es diferente en cada tienda.

eta2 <- 5.737/(5.737+3.964)
eta2
## [1] 0.5913823

Al estar por encima del 30%, indica que el tamaño del efecto es grande (por eso se rechaza H0, y las medias son diferentes).

El análisis de varianza resulta significativo con un tamaño de efecto grande, pero al no haberse comprobado la condición de esfericidad no se pueden considerar válidos.

Con library ANOVA

Realizar un ANOVA de datos pareados con la función Anova() requiere tener los datos almacenados en un formato específico. Los valores de la variable continua se almacenan en una matriz en la que cada columna representa un nivel distinto (grupo) de la variable cualitativa.

datos <- as.matrix(datos[-1])
datos
##       tienda_A tienda_B tienda_C tienda_D
##  [1,]    1.755     1.78     1.29     1.29
##  [2,]    2.655     1.98     1.99     1.99
##  [3,]    2.235     1.69     1.79     1.59
##  [4,]    0.975     0.99     0.69     1.09
##  [5,]    2.370     1.70     1.89     1.89
##  [6,]    4.695     3.15     2.99     3.09
##  [7,]    3.135     1.88     2.09     2.49
##  [8,]    0.930     0.65     0.65     0.69
##  [9,]    8.235     5.99     5.99     6.99
## [10,]    6.690     4.84     4.99     5.15

El siguiente paso es crear un modelo lineal multivariable. Los coeficientes estimados coinciden con la media de cada grupo.

modelo_lm <- lm(datos~1)
summary(modelo_lm)
## Response tienda_A :
## 
## Call:
## lm(formula = tienda_A ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4375 -1.4925 -0.8550  0.9375  4.8675 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.3675     0.7717   4.364  0.00181 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.44 on 9 degrees of freedom
## 
## 
## Response tienda_B :
## 
## Call:
## lm(formula = tienda_B ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8150 -0.7725 -0.6350  0.3925  3.5250 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.4650     0.5399   4.565  0.00136 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.707 on 9 degrees of freedom
## 
## 
## Response tienda_C :
## 
## Call:
## lm(formula = tienda_C ~ 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.786 -1.021 -0.496  0.329  3.554 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.4360     0.5582   4.364  0.00181 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.765 on 9 degrees of freedom
## 
## 
## Response tienda_D :
## 
## Call:
## lm(formula = tienda_D ~ 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.936 -1.261 -0.686  0.314  4.364 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.6260     0.6286   4.178  0.00238 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.988 on 9 degrees of freedom
tienda <- factor(c("tienda_A", "tienda_B", "tienda_C", "tienda_D"))

Por último se emplea la función Anova() para realizar el análisis de varianza, el test de esfericidad y las correcciones.

library(car)
## Warning: package 'car' was built under R version 4.2.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.2.3
anova_pareado <- Anova(modelo_lm, idata = data.frame(tienda),
                       idesign = ~ tienda, type = "III")
summary(anova_pareado, multivariate = F)
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##              Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept) 296.725      1  139.479      9  19.146  0.001782 ** 
## tienda        5.737      3    3.964     27  13.025 1.898e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##        Test statistic   p-value
## tienda        0.12901 0.0078973
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##         GG eps Pr(>F[GG])   
## tienda 0.46824   0.001747 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           HF eps Pr(>F[HF])
## tienda 0.5287146 0.00103603

Ante la falta de esfericidad, se analiza la Greenhouse-Geisser.

El estadístico HF asciende a 0.52 con p valor de 0.001 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Por tanto, las medias de las cestas de la compra son diferentes en cada tienda.

##2.3.- Comparaciones múltiples

No es posible crear los intervalos TUKEY HSD para datos pareados. Para las comparaciones múltiples después de que el ANOVA haya resultado significativo hay que recurrir a otrsas correcciones.

pairwise.t.test(x = datos_tabla_larga$precio, g = datos_tabla_larga$tienda,
                p.adjust.method = "holm", paired = TRUE, alternative = "two.sided")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  datos_tabla_larga$precio and datos_tabla_larga$tienda 
## 
##          tienda_A tienda_B tienda_C
## tienda_B 0.022    -        -       
## tienda_C 0.014    0.695    -       
## tienda_D 0.014    0.491    0.331   
## 
## P value adjustment method: holm

A-B: El p valor es 0.022 < 0.05. Las medias son diferentes.

A-C: El p valor es 0.014 < 0.05. Las medias son diferentes.

A-D: El p valor es 0.014 < 0.05. Las medias son diferentes.

B-C: El p valor es 0.695 > 0.05. Las medias son similares.

B-D: El p valor es 0.491 > 0.05. Las medias son similares.

C-D: El p valor es 0.331 > 0.05. Las medias son similares.

#2.4.- CONCLUSIÓN

El análisis ANOVA (incluyendo las correcciones, dada la falta de esfericidad) encuentra diferencias significativas en el precio de los alimentos entre al menos 2 tiendas, siendo el tamaño del efecto grande. La posterior comparación dos a dos por t-student con corrección de significancia holm identifica como significativas las diferencias entre las tiendas A-B, A-C y A-D pero no entre B-C, B-D y C-D.

##2.5.- INFO DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 70 × 3
##    package    loadedversion source        
##    <chr>      <chr>         <chr>         
##  1 abind      1.4-5         CRAN (R 4.2.0)
##  2 bslib      0.4.2         CRAN (R 4.2.2)
##  3 cachem     1.0.6         CRAN (R 4.2.2)
##  4 callr      3.7.3         CRAN (R 4.2.3)
##  5 car        3.1-2         CRAN (R 4.2.3)
##  6 carData    3.0-5         CRAN (R 4.2.3)
##  7 cli        3.6.0         CRAN (R 4.2.2)
##  8 colorspace 2.1-0         CRAN (R 4.2.3)
##  9 crayon     1.5.2         CRAN (R 4.2.3)
## 10 devtools   2.4.5         CRAN (R 4.2.3)
## # ℹ 60 more rows