29.5: Análisis de varianza (Sección 28.6.1)

A menudo queremos comparar varios medios diferentes, para determinar si alguno de ellos es diferente de los demás. En este caso, veamos los datos de NHANES para determinar si el estado civil está relacionado con la calidad del sueño. Primero limpiamos los datos:

NHANES_sleep_marriage <- 
  NHANES_adult %>%
  dplyr::select(SleepHrsNight, MaritalStatus, Age) %>%
  drop_na()

En este caso vamos a tratar el conjunto completo de datos NHANES como nuestra muestra, con el objetivo de generalizar a toda la población estadounidense (de la cual el conjunto de datos NHANES es medio para ser una muestra representativa). Primero veamos la distribución de los diferentes valores de la variable maritalStatus:

NHANES_sleep_marriage %>%
  group_by(MaritalStatus) %>%
  summarize(n=n()) %>%
  kable()

MaritalStatus	n
divorciados	437
LivePartner	370
Casado	2434
Nunca Casado	889
Separados	134
Viudo	329

Hay números razonables de la mayoría de estas categorías, pero eliminemos la categoría Separados ya que tiene relativamente pocos miembros:

NHANES_sleep_marriage <-
  NHANES_sleep_marriage %>%
  dplyr::filter(MaritalStatus!="Separated")

Ahora usemos lm () para realizar un análisis de varianza. Dado que también sospechamos que la Edad está relacionada con la cantidad de sueño, también incluiremos Edad en el modelo.

lm_sleep_marriage <- lm(SleepHrsNight ~ MaritalStatus + Age,
                        data=NHANES_sleep_marriage)
summary(lm_sleep_marriage)

## 
## Call:
## lm(formula = SleepHrsNight ~ MaritalStatus + Age, data = NHANES_sleep_marriage)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.016 -0.880  0.107  1.082  5.282 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                6.51758    0.09802   66.49  < 2e-16 ***
## MaritalStatusLivePartner   0.14373    0.09869    1.46  0.14536    
## MaritalStatusMarried       0.23494    0.07094    3.31  0.00093 ***
## MaritalStatusNeverMarried  0.25172    0.08404    3.00  0.00276 ** 
## MaritalStatusWidowed       0.26304    0.10327    2.55  0.01090 *  
## Age                        0.00318    0.00141    2.25  0.02464 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.4 on 4453 degrees of freedom
## Multiple R-squared:  0.00458,    Adjusted R-squared:  0.00347 
## F-statistic:  4.1 on 5 and 4453 DF,  p-value: 0.00102

Esto nos dice que existe un efecto muy significativo del estado civil (basado en la prueba F), aunque representa una cantidad muy pequeña de varianza (menos del 1%).

También es útil observar con más detalle qué grupos difieren de cuáles otros, lo que podemos hacer examinando las medias marginales estimadas para cada grupo utilizando la función emmeans ().

# compute the differences between each of the means
leastsquare <- emmeans(lm_sleep_marriage, 
                      pairwise ~ MaritalStatus,
                      adjust="tukey")
 
# display the results by grouping using letters

CLD(leastsquare$emmeans, 
    alpha=.05,  
    Letters=letters)

##  MaritalStatus emmean    SE   df lower.CL upper.CL .group
##  Divorced         6.7 0.066 4453      6.5      6.8  a    
##  LivePartner      6.8 0.073 4453      6.7      7.0  ab   
##  Married          6.9 0.028 4453      6.8      7.0   b   
##  NeverMarried     6.9 0.050 4453      6.8      7.0   b   
##  Widowed          6.9 0.082 4453      6.8      7.1  ab   
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05

Las letras en la columna de grupo nos indican qué condiciones individuales difieren de cuáles otras; cualquier par de condiciones que no comparten un identificador de grupo (en este caso, las letras a y b) son significativamente diferentes entre sí. En este caso, vemos que las personas divorciadas duermen menos que las personas casadas o viudas; ninguna otra pareja difiere significativamente.

29.5.1 Análisis de varianza de medidas repetidas

El análisis estándar de varianza supone que las observaciones son independientes, lo que debería ser cierto para diferentes personas en el conjunto de datos de NHANES, pero puede no ser cierto si los datos se basan en medidas repetidas del mismo individuo. Por ejemplo, el conjunto de datos NHANES implica tres mediciones de la presión arterial para cada individuo. Si queremos probar si hay alguna diferencia entre esos, entonces necesitaríamos usar un análisis de varianza de medidas repetidas. Podemos hacer esto usando lmer () como hicimos anteriormente. Primero, necesitamos crear una versión “larga” del conjunto de datos.

NHANES_bp_all <- NHANES_adult %>%
  drop_na(BPSys1,BPSys2,BPSys3) %>%
  dplyr::select(BPSys1,BPSys2,BPSys3, ID) %>%
  gather(test, BPsys, -ID)

Entonces encajamos un modelo que incluye una intercepción separada para cada individuo.

repeated_lmer <-lmer(BPsys ~ test + (1|ID), data=NHANES_bp_all)
summary(repeated_lmer)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BPsys ~ test + (1 | ID)
##    Data: NHANES_bp_all
## 
## REML criterion at convergence: 89301
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.547 -0.513 -0.005  0.495  4.134 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 280.9    16.8    
##  Residual              16.8     4.1    
## Number of obs: 12810, groups:  ID, 4270
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  122.0037     0.2641 4605.7049   462.0   <2e-16 ***
## testBPSys2    -0.9283     0.0887 8538.0000   -10.5   <2e-16 ***
## testBPSys3    -1.6215     0.0887 8538.0000   -18.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) tsBPS2
## testBPSys2 -0.168       
## testBPSys3 -0.168  0.500

Esto nos muestra que la segunda y tercera prueba son significativamente diferentes de la primera prueba (la cual se asignó automáticamente como línea base por lmer ()). También podríamos querer saber si hay un efecto general de la prueba. Podemos determinar esto comparando el ajuste de nuestro modelo con el ajuste de un modelo que no incluya la variable de prueba, que encajaremos aquí. Luego comparamos los modelos usando la función anova (), que realiza una prueba F para comparar los dos modelos.

repeated_lmer_baseline <-lmer(BPsys ~ (1|ID), data=NHANES_bp_all)
anova(repeated_lmer,repeated_lmer_baseline)

## Data: NHANES_bp_all
## Models:
## repeated_lmer_baseline: BPsys ~ (1 | ID)
## repeated_lmer: BPsys ~ test + (1 | ID)
##                        Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## repeated_lmer_baseline  3 89630 89652 -44812    89624                        
## repeated_lmer           5 89304 89341 -44647    89294   330      2     <2e-16
##                           
## repeated_lmer_baseline    
## repeated_lmer          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Esto demuestra que la presión arterial difiere significativamente entre las tres pruebas.