library(tidyverse)
library(ggpubr)
library(rstatix)
library(broom)
library(emmeans)Ejemplo de Análisis de Covarianza: Control de la Presión Sanguínea Sistólica mediante Edad
Introducción
Se presenta un Análisis de Covarianza para verificar si la Presión Sanguínea Sistólica (SPB) varía según el Sexo (SEX), controlando por la covariable Edad (AGE).
Los datos están disponibles en el libro de Kleinbaum, Kupper, Nizam & Rosenberg (2014) y contiene las siguientes variables:
OBS: Número de Observación.
SEX: Sexo (Male (Masculino) y Female (Femenino)).
SBP: Presión Sanguínea Sistólica.
AGE: Edad en años.
Carga de los paquetes
Preparación de los datos
Lectura de los datos
datos <- read.table("Kleinbaum_Ancova.txt", header = TRUE)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(datos$SBP ~ datos$SEX, xlab = "Sexo", ylab = "Presión Sanguínea Sistólica")
Se puede observar que la Presión Sanguínea Sistólica es ligeramente superior para los hombres.
Prueba t para comparación de medias
Mediante la prueba t:
t.test(datos$SBP ~ datos$SEX, alternative = "less")
Welch Two Sample t-test
data: datos$SBP by datos$SEX
t = -3.6622, df = 58.47, p-value = 0.0002698
alternative hypothesis: true difference in means between group Female and group Male is less than 0
95 percent confidence interval:
-Inf -8.310877
sample estimates:
mean in group Female mean in group Male
139.8621 155.1500 Con lo cual podemos concluir que el promedio de la Presión Sanguínea Sistólica de los hombres es significativamente mayor que el de las mujeres (p-valor = 0,0002698).
Análisis de Covarianza
Con la finalidad de mejorar el ajuste, se efectúa un Análisis de Covarianza incorporando la covariable Edad (AGE).
Chequeo de los supuestos
Linealidad
ggscatter(
datos, x = "AGE", y = "SBP",
color = "SEX", add = "reg.line"
)+
stat_regline_equation(
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~"), color = SEX)
)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 pareciera existir una relación lineal entre la la Edad y la Preión Sanguínea Sistólica para cada Sexo.
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:
datos %>% anova_test(SBP ~ SEX*AGE)ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 ges
1 SEX 1 65 38.221 4.69e-08 * 0.370000
2 AGE 1 65 175.958 3.68e-20 * 0.730000
3 SEX:AGE 1 65 0.007 9.34e-01 0.000106En donde se observa que la interacción no es significativa (p-valor = 0,934), 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(SBP ~ AGE + SEX, data = datos)
# 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
SBP AGE SEX .resid .cooksd .std.resid
<int> <int> <chr> <dbl> <dbl> <dbl>
1 158 41 Male 8.51 0.00884 0.972
2 185 60 Male 17.3 0.0495 1.99
3 152 41 Male 2.51 0.000771 0.287Normalidad 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.956 0.0157Y al ser el p-valor (0,0105) menor que 0,05, concluimos que pareciera violarse el supuesto de normalidad de los residuos, por lo que se suguiere usar una alternativa no paramétrica.
Homogeneidad de las varianzas (homocedasticidad)
La homogeneidad de las varianza puede verificarse mediante el test de Levene:
metricas.modelo %>% levene_test(.resid ~ SEX)Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.# A tibble: 1 × 4
df1 df2 statistic p
<int> <int> <dbl> <dbl>
1 1 67 0.693 0.408Y al ser el p-valor (0,408) mayor 0,05 podemos concluir que no se viola el supuesto de homogeneidad de las varianzas para los dos sexos (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] SBP AGE SEX .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 pudiera haber problemas con 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 <- datos %>% anova_test(SBP ~ AGE + SEX)
get_anova_table(res.aov)ANOVA Table (type II tests)
Effect DFn DFd F p p<.05 ges
1 AGE 1 66 178.646 1.9e-20 * 0.73
2 SEX 1 66 38.805 3.7e-08 * 0.37Y en la misma se puede observar que luego del ajuste para la covariable (Presión Sanguínea Sistólica), existe diferencia significativa entre los sexos (p-valor = 3,7 x 10-8).
En este caso no es necesario realizar pruebas a posteriori, pero creamos el objeto para la obtención de los ajustes:
# Prueba a posteriori
pwc <- datos %>%
emmeans_test(
SBP ~ SEX, covariate = AGE,
p.adjust.method = "bonferroni"
)
pwc# A tibble: 1 × 9
term .y. group1 group2 df statistic p p.adj p.adj.signif
* <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 AGE*SEX SBP Female Male 66 -6.23 0.0000000370 3.70e-8 **** 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: 2 × 8
AGE SEX emmean se df conf.low conf.high method
<dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 46.1 Female 141. 1.65 66 138. 144. Emmeans test
2 46.1 Male 154. 1.40 66 152. 157. Emmeans testCon lo que podemos concluir que la media la Presión Sanguínea Sistólica es significativamente más alta para el los hombre (154,4 +/- 1,4) que para las mujeres (140,9 +/- 1,7).
Referencias
Kleinbaum, D., Kupper, L., Nizam, A., & Rosenberg, E. (2014). Applied Regression Analysis and Other Multivariable Methods. %ta. edición. Boston, MA, USA : Cengage Learning.