Debido a que ANCOVA es un modelo híbrido entre regresión lineal y ANOVA, la función de R para hacerlo es lm. En clase introductoria de ANCOVA vimos las maneras de especificar el modelo, los significados de cada tipo de especificación y la representación visual. En este tutorial analizaremos la base de datos stress del paquete datarium.

La base stress contiene datos de un experimento clínico sobre los efectos de un tratamiento farmaceutico, ejercicio físico y la edad sobre los niveles de estrés. El análisis comenzará entonces con una revisión rápida de la base de datos y la medición de las características estadísticas de las variables.

estres <- datarium::stress
str(estres)
## Classes 'tbl_df', 'tbl' and 'data.frame':    60 obs. of  5 variables:
##  $ id       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ score    : num  95.6 82.2 97.2 96.4 81.4 83.6 89.4 83.8 83.3 85.7 ...
##   ..- attr(*, "label")= chr "Cholesterol concentration (in mmol/L)"
##   ..- attr(*, "format.spss")= chr "F8.2"
##   ..- attr(*, "display_width")= int 9
##  $ treatment: Factor w/ 2 levels "yes","no": 1 1 1 1 1 1 1 1 1 1 ...
##  $ exercise : Factor w/ 3 levels "low","moderate",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ age      : num  59 65 70 66 61 65 57 61 58 55 ...
##   ..- attr(*, "label")= chr "Body weight (in kg)"
##   ..- attr(*, "format.spss")= chr "F8.2"

Con la función str (structure) podemos dar un vistazo al nombre y tipo de variables contenidas en cada columna. La base estrés contiene cinco columnas, cuyos nombres aparecen después del operador $ junto con el tipo. La columna score contiene es la variable de respuesta y consiste de la medición de los niveles de estrés, mientras que las demás contienen las variables independientes categóricas (treatment y exercise) y contínuas (age).

Debido a la presencia de tres variables independientes podemos intuir que se trata de un experimento de tres vías, y debido a que no hay mediciones repetidas o bloques experimentales se tratará de un análisis con efectos fijos únicamente. En esta ocasión nos saltaremos la evaluación previa de la distribución estadística de age, y comenzaremos por verificar aspectos pertinentes a un ANCOVA:

  1. Potenciales diferencias entre los interceptos (\(\beta_0\)) de cada grupo (treatment y exercise).
  2. Potenciales diferencies entre las pendientes (\(\beta_1\)) de cada grupo en relación al factor contínuo (age).

0.1 Análisis exploratorio

Antes de continuar con el análisis visual de las asociaciones entre las diferentes variables, discutamos brevemente el primer punto. Sabemos que el intercepto es el valor de \(y\) cuando \(x=0\), y debido a que age tiene valores de ~45 - 75, es poco probable que logremos verificar esto visualmente. Sin embargo, sabemos que si las líneas de regresión para treatment:yes y treatment:no son paralelas habrá un efecto de treatment sobre el intercepto únicamente. Y, si las líneas no son paralelas pero se cruzan en un punto donde $ x $ también hay un efecto de treatment (o exercise en su defecto) sobre el intercepto.

library(ggplot2)

ggplot(estres) + geom_point(aes(x = age, y = score, colour = treatment)) + 
      geom_smooth(aes(x = age, y = score, colour = treatment, fill = treatment), alpha = 0.3, method = "lm") + 
      theme_bw()
Efectos de edad y tratamiento.

Figure 1: Efectos de edad y tratamiento.

Las líneas de regresión no son paralelas y se cruzan en un punto donde \(x \neq 0\), por lo tanto podemos inferir que treatment sí afectará al intercepto. Y con respecto a exercise:

ggplot(estres) + geom_point(aes(x = age, y = score, colour = exercise)) + 
      geom_smooth(aes(x = age, y = score, colour = exercise, fill = exercise), alpha = 0.3, method = "lm") + 
      theme_bw()
Efectos de edad y ejercicio

Figure 2: Efectos de edad y ejercicio

Observamos efectos diferentes, para comenzar las líneas de regresión de exercise:moderate y exercise:high tienen pendientes muy similares (están prácticamente paralelas), y la línea de exercise:low cruza con la línea de exercise:moderate en un punto donde \(x \neq 0\). Para resumir, en el caso de exercise podemos también inferir que hay efectos del ejercicio sobre el intercepto (\(\beta_0\)).

En cuanto a las interacciones, queda más claro que treatment sí interactúa con age, mientras que tal efecto para exercise es menos claro. Queda entonces por explorar si hay una posible interacción treatment * exercise, hagamos para ello las siguientes gráficas de cajas con ggplot2:

ggplot(estres) + geom_boxplot(aes(x = treatment, y = score, colour = exercise, fill = exercise), alpha = 0.3) + 
      facet_wrap(~exercise) +
      theme_bw()
Interacción tratamiento-ejercicio.

Figure 3: Interacción tratamiento-ejercicio.

Las gráficas de cajas sugieren que los niveles de estrés son más bajos cuando hay tratamiento farmacológico que cuando no lo hay y que además los niveles de estrés son menores cuano el tratamiento se combina con ejercicios de mayor intensidad, lo cual sugiere que sí hay interacción entre ambos factores. Con base en esto propongo comenzar en análisis con el modelo más complicado posible:

age ~ treatment * exercise * age

y posteriormente mediante el análisis de los efectos estimados y las sumas de cuadrados de cada término podremos decidir si se elimina dicho término de un modelo subsecuente.

0.2 Ajustando los modelos

m1 <- lm(score ~ treatment * exercise * age, estres)
summary(m1)
## 
## Call:
## lm(formula = score ~ treatment * exercise * age, data = estres)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.2940 -3.4969  0.3342  2.2295 10.2988 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                       59.18087   22.87155   2.588   0.0127 *
## treatmentno                        8.41860   33.82489   0.249   0.8045  
## exercisemoderate                  -1.05570   46.02625  -0.023   0.9818  
## exercisehigh                     -16.40905   50.25374  -0.327   0.7454  
## age                                0.46482    0.36973   1.257   0.2148  
## treatmentno:exercisemoderate     -18.45149   55.34236  -0.333   0.7403  
## treatmentno:exercisehigh          12.85565   63.49342   0.202   0.8404  
## treatmentno:age                   -0.11070    0.54501  -0.203   0.8399  
## exercisemoderate:age               0.01664    0.76488   0.022   0.9827  
## exercisehigh:age                   0.04428    0.86727   0.051   0.9595  
## treatmentno:exercisemoderate:age   0.30217    0.91128   0.332   0.7416  
## treatmentno:exercisehigh:age      -0.08848    1.08428  -0.082   0.9353  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.204 on 48 degrees of freedom
## Multiple R-squared:  0.6657, Adjusted R-squared:  0.5891 
## F-statistic: 8.689 on 11 and 48 DF,  p-value: 3.358e-08

El modelo es demasiado complicado y ningún término con excepción del intercepto es significativo, por lo que debemos comenzar a tomar decisiones. Para ello, nos ayudaremos de la tabla ANOVA del modelo y veremos qué términos explicaron mayor variabilidad:

anova(m1)
## Analysis of Variance Table
## 
## Response: score
##                        Df  Sum Sq Mean Sq F value   Pr(>F)    
## treatment               1  351.38  351.38 12.9757 0.000747 ***
## exercise                2 1776.27  888.13 32.7965 1.05e-09 ***
## age                     1  222.74  222.74  8.2252 0.006121 ** 
## treatment:exercise      2  220.94  110.47  4.0793 0.023109 *  
## treatment:age           1    0.38    0.38  0.0142 0.905701    
## exercise:age            2   12.75    6.38  0.2355 0.791115    
## treatment:exercise:age  2    3.96    1.98  0.0730 0.929667    
## Residuals              48 1299.85   27.08                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como bien intuimos al principio, hay poca evidencia de la interacción exercise:age, y por ende de la itneracción treatment:exercise:age, por lo que podemos proponer el modelo simplificado sin esas interacciones. Originalmente, sospechábamos de la posible interaccción treatment:age, así que la dejaremos para ver qué pasa con su efecto una vez que eliminemos las otras interacciones en el modelo:

m2 <- lm(score ~ treatment * (exercise + age), estres)
summary(m2)
## 
## Call:
## lm(formula = score ~ treatment * (exercise + age), data = estres)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.2844 -2.9981  0.2526  2.4911 10.6950 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   58.57851   17.92087   3.269  0.00192 ** 
## treatmentno                   -1.18118   22.11609  -0.053  0.95761    
## exercisemoderate              -0.04339    2.33083  -0.019  0.98522    
## exercisehigh                 -13.83948    2.62917  -5.264 2.72e-06 ***
## age                            0.47458    0.28930   1.640  0.10696    
## treatmentno:exercisemoderate   0.22627    3.24307   0.070  0.94465    
## treatmentno:exercisehigh       8.41676    3.56832   2.359  0.02213 *  
## treatmentno:age                0.04382    0.35580   0.123  0.90245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.032 on 52 degrees of freedom
## Multiple R-squared:  0.6614, Adjusted R-squared:  0.6158 
## F-statistic: 14.51 on 7 and 52 DF,  p-value: 2.691e-10

Aquí comenzamos a ver que sí hay coeficientes significativos, lo cual es buena noticia. Veamos entonces si hay más términos que debemos eliminar del modelo:

anova(m2)
## Analysis of Variance Table
## 
## Response: score
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## treatment           1  351.38  351.38 13.8786 0.0004814 ***
## exercise            2 1776.27  888.13 35.0786 2.271e-10 ***
## age                 1  222.74  222.74  8.7976 0.0045475 ** 
## treatment:exercise  2  220.94  110.47  4.3632 0.0177124 *  
## treatment:age       1    0.38    0.38  0.0152 0.9024520    
## Residuals          52 1316.56   25.32                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El único término que, aparentemente deberíamos eliminar es la interacción treatment:age, así que tendremos un modelo con interacción y efectos aditivos entre treatment y exerceise. Ahora, veamos directamente la tabla ANOVA:

m3 <- lm(score ~ treatment * exercise + age, estres)
anova(m3)
## Analysis of Variance Table
## 
## Response: score
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## treatment           1  351.38  351.38 14.1414 0.0004249 ***
## exercise            2 1776.27  888.13 35.7428 1.488e-10 ***
## age                 1  222.74  222.74  8.9641 0.0041770 ** 
## treatment:exercise  2  220.94  110.47  4.4458 0.0164087 *  
## Residuals          53 1316.94   24.85                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Y podemos confirmar nuestra selección de modelo comparado los tres:

anova(m1, m2, m3)
## Analysis of Variance Table
## 
## Model 1: score ~ treatment * exercise * age
## Model 2: score ~ treatment * (exercise + age)
## Model 3: score ~ treatment * exercise + age
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     48 1299.8                           
## 2     52 1316.6 -4  -16.7084 0.1542 0.9602
## 3     53 1316.9 -1   -0.3841 0.0142 0.9057

Vemos que la suma de cuadrados de los residualtes es mayor en el tercer modelo, es que elegimos como el más adecuado. Esto es debido a que un modelo más complicado siempre va a tender a explicar mayor variabilidad que otro modelo más sencillo. Necesitamos entonces una manera de tomar en cuenta la complejidad del modelo en la comparación, y una manera estándar de hacerlo en estadística es con el Criterio de Información de Akaike, (AIC , por sus siglas en inglés), cuya fórmula es:

\[ \mathrm{AIC} = 2 n \ln{\mathcal{L}}\] donde \(n\) es el número de parámetros o coeficientes estimados por el modelo (medida de complejidad), y \(\mathcal{L}\) es la verosimilitud. En modelos lineales (regresión, ANOVA y ANCOVA), la verosimilitud es la suma de cuadrados del modelo, de modo que un AIC más pequeño es mejor (hay excepciones, pues AIC no sustituye al análisis posterior). Para calcular el AIC de un modelo en R, hacemos:

AIC(m1, m2, m3)
##    df      AIC
## m1 13 380.8121
## m2  9 373.5784
## m3  8 371.5959

Y podemos ver que m3, con base en este criterio, es mejor que m1. Este proceso de selección del modelo lo terminaremos con el diagnóstico de los residuales.

0.3 Diagnóstico posterior

El diangóstico de modelos ANCOVA es igual que para modelos de regresión lineal y ANOVA:

par(mfrow = c(2,2))
plot(m3)
Gráficas de diagnóstico de `m3`.

Figure 4: Gráficas de diagnóstico de m3.

Basándonos en las gráficas de la primera fila, podemos suponer que los residuales están bien contenidos, con la pequeña execpción de los residuales para el score = 95. Y para ver si esta desviación podría resultar en problemas para probar las hipótesis veamos la distribución de los residuales:

shapiro.test(residuals(m3))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m3)
## W = 0.98226, p-value = 0.5313

Y vemos que sí es posible confiar, con base en los supuestos estadísticos del modelo ANCOVA en las pruebas de hipótesis pertinentes a m3:

  1. El efecto del tratamiento depende de que se practique ejercicio (interacción)
  2. La edad tenderá a incrementar los niveles de estrés, independientemente del tratamiento y ejercicio (efecto aditivo)

El script ANCOVA-selección.R tiene ejemplos de código para combinar los tipos de modelos vistos en esta unidad.

Regresar al índice del curso