Ejemplo de Análisis de Covarianza: Control de la ansiedad mediante ejercicio físico

Author

Rafael E. Borges

Introducción

Se presenta un análisis de covarianza para estudiar el efecto del ejercicio físico (bajo (grp1), medio (grp2) y alto (grp3)) sobre la ansiedad, usando como covariable la ansiedad medida antes del programa de entrenamiento físico de seis meses.

El análisis original de los datos se presenta en el enlace: https://www.datanovia.com/en/lessons/ancova-in-r/

Carga de los paquetes

library(tidyverse)
library(ggpubr)
library(rstatix)
library(broom)
library(emmeans)

Preparación de los datos

Los datos (anxiety) están disponibles en el paquete datarium.

# Carga y preparación de los datos
data("anxiety", package = "datarium")
anxiety <- anxiety %>%
  select(id, group, t1, t3) %>%
  rename(pretest = t1, posttest = t3)
anxiety[14, "posttest"] <- 19
# Inspección de los datos mostrando una fila aleatoria por grupo:
set.seed(123)
anxiety %>% sample_n_by(group, size = 1)

# A tibble: 3 × 4
  id    group pretest posttest
  <fct> <fct>   <dbl>    <dbl>
1 15    grp1     19.8     19.4
2 30    grp2     19.3     17.7
3 33    grp3     15.5     11

Ajuste del modelo sin corregir el efecto de la covariable

Gráfico de cajas múltipes de la ansiedad según grupo

En el gráfico de cajas multiples:

boxplot(anxiety$posttest ~ anxiety$group)

Se puede observar que el nivel de ansiedad posterior al programa de entrenamiento físico pareciera disminuir a medida que aumenta la actividad física, pero no pareciera existir diferencia entre los grupos 1 y 2.

Ajuste del modelo sin la covariable

Al ajustar el modelo si la covariable:

anova(lm(anxiety$posttest ~ anxiety$group))

Analysis of Variance Table

Response: anxiety$posttest
              Df  Sum Sq Mean Sq F value    Pr(>F)    
anxiety$group  2  72.134  36.067  13.996 2.199e-05 ***
Residuals     42 108.229   2.577                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se puede observar que existe diferencia significativa en las medias del nivel de ansiedad según los grupos (p-valor = 2,199 x 10^-5).

Prueba a posteriori

En la prueba a posteriori:

# Prueba a posteriori (con el paquete emmmeans)
pwc <- anxiety %>% 
  emmeans_test(
    posttest ~ group,
    p.adjust.method = "bonferroni"
    )
pwc

# A tibble: 3 × 9
  term  .y.      group1 group2    df statistic          p     p.adj p.adj.signif
* <chr> <chr>    <chr>  <chr>  <dbl>     <dbl>      <dbl>     <dbl> <chr>       
1 group posttest grp1   grp2      42      1.87 0.0691     0.207     ns          
2 group posttest grp1   grp3      42      5.22 0.00000519 0.0000156 ****        
3 group posttest grp2   grp3      42      3.36 0.00169    0.00507   **

Se puede observar que no existe diferencia significativa del nivel de ansiedad entre los grupos 1 y 2.

Análisis de Covarianza

Se efectúa un Análisis de Covarianza incorporando la covariable pretest (nivel de ansiedad antes del programa de entrenamiento).

Chequeo de los supuestos

Linealidad

ggscatter(
  anxiety, x = "pretest", y = "posttest",
  color = "group", add = "reg.line"
  )+
  stat_regline_equation(
    aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = group)
    )

Warning: The dot-dot notation (`..eq.label..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(eq.label)` instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.

En el gráfico anterior se puede observar que existe una relación lineal entre la ansiedad antes y después del programa de entrenamiento físico para cada grupo.

Homogeneidad de las pendientes de las rectas de regresión

La homogeneidad de las rectas de regresión puede chequearse verificando la no significación del término de interacción entre la covariable y la variable de grupo en la tabla Anova (tipo II) del modelo saturado:

anxiety %>% anova_test(posttest ~ group*pretest)

ANOVA Table (type II tests)

         Effect DFn DFd       F        p p<.05   ges
1         group   2  39 209.314 1.40e-21     * 0.915
2       pretest   1  39 572.828 6.36e-25     * 0.936
3 group:pretest   2  39   0.127 8.81e-01       0.006

En donde se observa que la interacción no es significativa (p-valor = 0,88), con lo que se comprueba la homogeneidad de las pendientes de las rectas de regresión.

Ajuste del modelo

Para verificar los supuestos de normalidad de los residuos, homogeneidad de las varianzas y la no presencia de datos atípicos, lo primero que hay que hacer es ajustar un modelo lineal usando la función lm y posteriormente generar las métrica mediante la función augment(model) del paquete broom:

# Ajuste del modelo (la covarianza va primero)
modelo <- lm(posttest ~ pretest + group, data = anxiety)
# Inspeccón de las métricas de diagnóstico del modelo
metricas.modelo <- augment(modelo) %>%
  select(-.hat, -.sigma, -.fitted) # Remove details
# Visualización de los primeros 3 datos:
head(metricas.modelo, 3)

# A tibble: 3 × 6
  posttest pretest group .resid .cooksd .std.resid
     <dbl>   <dbl> <fct>  <dbl>   <dbl>      <dbl>
1     14.1    14.1 grp1   0.550  0.101       1.46 
2     14.3    14.5 grp1   0.338  0.0310      0.885
3     14.9    15.7 grp1  -0.295  0.0133     -0.750

Normalidad de los residuos

La normalidad de los residuos puede verificarse a través de la prueba de Shapiro y Wilks:

# Chequeo de la normalidad dw los residuos usando el test de Shapiro y Wilks
shapiro_test(metricas.modelo$.resid)

# A tibble: 1 × 3
  variable               statistic p.value
  <chr>                      <dbl>   <dbl>
1 metricas.modelo$.resid     0.975   0.444

Y al ser el p-valor (0,444) mayor que 0,05, concluimos que no se viola el supuesto de normalidad de los residuos.

Homogeneidad de las varianzas (homocedasticidad)

La homogeneidad de las varianza puede verificarse mediante el test de Levene:

metricas.modelo %>% levene_test(.resid ~ group)

# A tibble: 1 × 4
    df1   df2 statistic     p
  <int> <int>     <dbl> <dbl>
1     2    42      2.27 0.116

Y al ser el p-valor (0,116) mayor 0,05 podemos concluir que no se viola el supuesto de homogeneidad de las varianzas para todos los grupos (homocedásticidad).

Detección de valores atípicos

La dectección de datos atípicos se puede hacer verificando si hay individios con residuos estandarizados con valor absoluto mayor a 3:

metricas.modelo %>% 
  filter(abs(.std.resid) > 3) %>%
  as.data.frame()

[1] posttest   pretest    group      .resid     .cooksd    .std.resid
<0 rows> (or 0-length row.names)

Y no se identifican individuos on residuos estandarizados con valor absoluto mayor a 3.

Verificación gráfica de los supuestos

La verificación de los dos supuestoa anteriores puede verse mediante el análisis de los residuos:

plot(modelo)

En el primer gráfico podemos observa que no se viola el supuesto de homogeneidad de las varianzas.

En el segundo gráfico se puede observar que no se viola el supuesto de normalidad de los residuos.

Y en el cuarto gráfico se puede observar que no existen datos atípicos.

Ajuste del modelo de Análisis de Covarianza:

La tabla Anova del Análisis de Covarianza puede ser obtenida eliminando primero el efecto de la covariable, lo cual se obtiene mediante:

res.aov <- anxiety %>% anova_test(posttest ~ pretest + group)
get_anova_table(res.aov)

ANOVA Table (type II tests)

   Effect DFn DFd       F        p p<.05   ges
1 pretest   1  41 598.321 4.48e-26     * 0.936
2   group   2  41 218.629 1.35e-22     * 0.914

Y en la misma se puede observar que luego del ajuste para la covariable (nivel de ansiedad antes del programa de entrenamiento físico), existe diferencia significativa entre los grupos (p-valor = 1,35 x 10^-22).

Prueba a posteriori

Y en la prueba a posteriori:

# Prueba a posteriori
pwc <- anxiety %>% 
  emmeans_test(
    posttest ~ group, covariate = pretest,
    p.adjust.method = "bonferroni"
    )
pwc

# A tibble: 3 × 9
  term        .y.   group1 group2    df statistic        p    p.adj p.adj.signif
* <chr>       <chr> <chr>  <chr>  <dbl>     <dbl>    <dbl>    <dbl> <chr>       
1 pretest*gr… post… grp1   grp2      41      4.24 1.26e- 4 3.77e- 4 ***         
2 pretest*gr… post… grp1   grp3      41     19.9  1.19e-22 3.58e-22 ****        
3 pretest*gr… post… grp2   grp3      41     15.5  9.21e-19 2.76e-18 ****

Se observa diferencia significativa en la media del nivel de ansiedad luego del programa de entrenamiento, entre todos los pares de grupos.

Visualización del ajuste del modelo

Las medias ajustadas mediante el modelo de Análisis de Covarianza pueden obtenerse mediante las medias marginales estimadas:

get_emmeans(pwc)

# A tibble: 3 × 8
  pretest group emmean    se    df conf.low conf.high method      
    <dbl> <fct>  <dbl> <dbl> <dbl>    <dbl>     <dbl> <chr>       
1    16.9 grp1    16.4 0.106    41     16.2      16.7 Emmeans test
2    16.9 grp2    15.8 0.107    41     15.6      16.0 Emmeans test
3    16.9 grp3    13.5 0.106    41     13.2      13.7 Emmeans test

Con lo que podemos concluir que la media del nivel de ansiedad es significativamente más alto para el grupo 1 (16,4 +/- 0,11) comparado con el grupo 2 (15,8 +/- 0,11) y el grupo 3 (13,5 +/- 0,11), con lo que se verifica la diferencia en la media del nivel de ansiedad posterior al entrenamiento para los tres grupos de ejercicio física del programa de entrenamiento (bajo, medio o alto), aspecto que no podía ser visualizado sin eliminar el efecto de la covariable.