Interpretando matrices de correlación en modelos multinivel con {`lme4`}

correlaciones_multinivel.R

# Interpretando matrices de correlación en modelos multinivel usando lme4
# Gonzalo García-Castro, gonzalo.garciadecastro@upf.edu

#### cargar paquetes #######################
library(dplyr)     # para manipular variables y usar pipes
library(lme4)      # para ajustar modelos mixtos
library(tibble)    # para pasar nombres de filas como una columna
library(tidyr)     # para pasar de una tabla en formato ancho a una en formato largo
library(janitor)   # para limpiar nombres de variables
library(ggplot2)   # para viualizar datos
library(patchwork) # para juntar gráficos

#### ajustar el modelo ######################
fit <- lmer(height ~ mass + (1 + mass + birth_year | gender) + (1 + mass + birth_year | hair_color), data = starwars)

# nuestro modelo incluye interceptos y slopes aleatorios pata "gender" y "hair_colour"
# por lo tanto, tendremos dos matriz de correlaciones: una para intercepts y slopes de "gender", y otra para "hair_colour"

# ATENCIÓN:
# en este caso, tenemos pocos datos y por lo tanto los coeficentes de correlación tendrán valores un poco extraños:
# por ejemplo, varianzas nulas(sd = 0) o correlaciones perfectas (corr = 1 o -1).
# esto nos lo dice la función lmer() con un aviso: "boundary (singular) fit: see ?isSingular" (del ue todo el mundo pasa;
# ver https://www.google.com/url?sa=i&url=https%3A%2F%2Fes-la.facebook.com%2FRmemes0%2Fphotos%2F&psig=AOvVaw3VxrvJeeJxNlNXcSPp6Ov2&ust=1589914271963000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCJj6ma2KvukCFQAAAAAdAAAAABAD)
# probablemente la estructura de efectos aleatorios sea demasiado compleja para los pocos datos que tenemos, y haya que simplificarla.
# no es un gran problema, porque la estimación de estos parámetros se realizar tras estimar todos los demás parámetros del modelo,
# pero hay que tener cuidado al interpretar estos coeficientes.
# sugerencia: Bates et al. (2015) https://arxiv.org/abs/1506.04967, donde se sugiere utilizar la función rePCA para identifiar grupos de parámetros que no explicar demasiada varianza y se pueden eliminar del modelo.

# pero hagamos como que todo está bien y sigamos

### matrices de var-cov y de correlación para el efecto aleatorio "gender"
var_corr <- VarCorr(fit) # esto nos da una matriz varianzas y covarianzas, y una matrix de correlaciones para cada efecto aleatorio

var_corr_gender <- VarCorr(fit)$gender # esto nos da lo anterior para el efecto aleatorio "gender"
varcov_gender <- attr(var_corr_gender, "stddev") # esto nos da la matriz de varianzas y covarianzas para "gender"
corr_gender <- attr(var_corr_gender, "correlation") # esto nos da la matriz de correlaciones para "gender"

# nos interesa interpretar la matriz de correlaciones de cada efecto aleatorio
# ¿Está asociado el tamaño de efecto de la masa sobre la altura a la altura media de cada género?
# Si es así, ¿cuál es la dirección de la asociación? Ejemplo:
#   - Si la asociación entre masa y altura es positiva, y afecta más a grupos más altos (e.j., género masculino), la correlación entre interceptos y slopes en el efecto aleatorio "género" será positiva.
#   - Si la asociación entre masa y altura es positiva, y afecta más a grupos más bajos (e.j., género masculino), la correlación entre interceptos y slopes en el efecto aleatorio "género" será negativa.

# visualización:
corr_gender_long <- corr_gender %>%
    as.data.frame() %>%
    clean_names() %>%
    rownames_to_column("term") %>%
    mutate(term = c("intercept", "mass", "birth_year")) %>%
    pivot_longer(cols = c(intercept, mass), names_to = "term2", values_to = "value")
    
plot_corr_gender <- ggplot(corr_gender_long, aes(x = term, y = term2, fill = value)) +
    geom_tile(colour = "white") +
    labs(fill = "Correlation") +
    theme_minimal() +
    theme(axis.title = element_blank(),
          legend.position = "top")

# ahora hacemos lo mismo para el efecto aleatorio "hair_color"
### matrices de var-cov y de correlación para el efecto aleatorio "hair_color"
var_corr_hair_color <- VarCorr(fit)$hair_color # esto nos da lo anterior para el efecto aleatorio "hair_color"
varcov_hair_color <- attr(var_corr_hair_color, "stddev") # esto nos da la matriz de varianzas y covarianzas para "hair_color"
corr_hair_color <- attr(var_corr_hair_color, "correlation") # esto nos da la matriz de correlaciones para "hair_color"

corr_hair_color_long <- corr_hair_color %>%
    as.data.frame() %>%
    clean_names() %>%
    rownames_to_column("term") %>%
    mutate(term = c("intercept", "mass", "birth_year")) %>%
    pivot_longer(cols = c(intercept, mass), names_to = "term2", values_to = "value")

# ahora visualizamos las dos matrices juntas
plot_corr_hair_color <- ggplot(corr_hair_color_long, aes(x = term, y = term2, fill = value)) +
    geom_tile(colour = "white") +
    theme_minimal() +
    theme(axis.title = element_blank(),
          legend.position = "none")

    
# juntamos los dos plots ¡Suerte con la interpretación! ;)
guide_area() / (plot_corr_gender + plot_corr_hair_color) +
    plot_layout(guides = "collect", heights = c(0.2, 0.8)) +
    plot_annotation(tag_levels = "A", title = "Correlation between intercepts and slopes for gender (A) and hair colour (B)")