class: center, middle, inverse, title-slide # Introducción a los modelos mixtos con R y tidymodels ### Alejandra Tapia ### R-Ladies Talca ### 19/03/2021 --- class: left, middle ## Antes de empezar 🧘🏻 #### > Este taller es un punto de partida para el uso de modelos mixtos con R y {tidymodels} #### > No se utilizan fórmulas matemáticas, como un intento de hacerlo amigable e intuitivo #### > Si surgen dudas, no te preocupes, es normal #### > El código y la presentación quedarán disponibles #### > Habrá espacio para preguntas --- class: center, middle ## ¿Qué es el modelamiento con efectos mixtos y por qué es importante?  --- class: center, middle ## `> Datos de las diferentes áreas del conocimiento pueden presentar heterogeneidad` --- class: center, middle ## `> Como consecuencia de estar agrupados, ser provenientes de estudios longitudinales o presentar algún otro tipo de medidas repetidas` --- class: center, middle ##`> Por ejemplo, datos de los capítulos de R-Ladies por continente` .center[<img src=imgs/R-Ladies_Chapters.png width="60%">] --- class: center, middle ##`> Datos de las integrantes de R-Ladies por capítulos` .center[<img src=imgs/R-Ladies_Woman.jpeg width="80%">] --- class: center, middle ## `> El tipo de medidas repetidas induce una estructura de correlación, que si no es considerada puede llevar a estimaciones sesgadas` --- class: center, middle ## `> Afectando inferencias y predicciones y, por tanto, la toma de decisiones basadas en datos (data-driven decision-making)`  --- class: left .pull-left[ ## Modelos lineales `>` Errores aleatorios independientes con varianza constante `>` Errores aleatorios siguen una distribución normal ] .pull-right[ ## Modelos lineales mixtos `>` Incorpora efectos aleatorios para acomodar la correlación entre las observaciones `>` Condicionado a los efectos aleatorios, los errores aleatorios son independientes con varianza constante `>` Errores aleatorios con distribución normal `>` Efectos aleatorios con distribución normal `>` Efectos aleatorios y errores aleatorios son independientes ] --- class: left ## Tidymodels <img src="imgs/hex_tidymodels.png" alt="Sharingan" width="6%" align="center"/> ### `>` .left[<img src=imgs/hex_broom.png width="15%">] ### `> {broom.mixed}` https://cran.r-project.org/web/packages/broom.mixed/index.html --- class: left ## Modelos mixtos 📦 ### `> {lme4}` https://cran.r-project.org/web/packages/lme4/index.html ### `> {lmerTest}` https://cran.r-project.org/web/packages/lmerTest/index.html --- class: center, middle ## Datos de dragones 🐉 ### .left[`> Datos ficticios de Gabriela Hajduk y taller inspirado en Athanasia Mowinckel`] --- class: center, middle ## Datos de dragones 🐉 ### .left[`> Decidimos entrenar dragones y, por tanto, salimos a diferentes sitios de las montañas y recopilamos datos sobre la inteligencia de los dragones`] --- class: center, middle ## Datos de dragones 🐉 ### .left[`> Inteligencia (testScore), longitud corporal (bodyLength), sitios (site) y cadenas montañosas (mountainRange)`] --- class: left, middle ## Cargar paquetes, leer y mirar los datos ```r library(tidyverse) library(readr) library(dplyr) ``` ```r url<-"https://raw.githubusercontent.com/alejandraandrea/slides-xaringan-mixed-models/master/dragons.tsv" download.file(url,"dragons.tsv") dragones<-read_tsv("dragons.tsv") glimpse(dragones) ``` ``` ## Rows: 480 ## Columns: 4 ## $ testScore <dbl> 16.147309, 33.886183, 6.038333, 18.838821, 33.862328, 4… ## $ bodyLength <dbl> 165.5485, 167.5593, 165.8830, 167.6855, 169.9597, 168.6… ## $ mountainRange <chr> "Bavarian", "Bavarian", "Bavarian", "Bavarian", "Bavari… ## $ site <chr> "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", … ``` --- class: center, middle ## Limpieza y transformación de datos ✅ --- class: left, middle ## Resumen preliminar de los datos ```r library(skimr) skim(dragones) %>% tibble::as_tibble() ``` ``` ## # A tibble: 4 x 17 ## skim_type skim_variable n_missing complete_rate character.min character.max ## <chr> <chr> <int> <dbl> <int> <int> ## 1 character mountainRange 0 1 6 8 ## 2 character site 0 1 1 1 ## 3 numeric testScore 0 1 NA NA ## 4 numeric bodyLength 0 1 NA NA ## # … with 11 more variables: character.empty <int>, character.n_unique <int>, ## # character.whitespace <int>, numeric.mean <dbl>, numeric.sd <dbl>, ## # numeric.p0 <dbl>, numeric.p25 <dbl>, numeric.p50 <dbl>, numeric.p75 <dbl>, ## # numeric.p100 <dbl>, numeric.hist <chr> ``` --- class: center, middle ## `Supongamos que queremos saber cómo afecta la longitud del cuerpo (bodyLength) de los dragones a sus puntuaciones en las pruebas (testScore)` --- class: left, middle ## Análisis exploratorio de datos 🕵 .pull-left[ ```r dragones %>% ggplot(aes(x=bodyLength, y=testScore)) + geom_jitter(alpha=.2)+ theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## `Si hay más covariables` 👀 `>` Se deben estandarizar `>` Garantiza que las estimaciones estén en la misma escala, facilitando la comparación de los tamaños de los efectos `>` Se puede usar la función scale() de {base}: .center[ ```r dragones$bodyLength2 <- scale(dragones$bodyLength) ``` ] --- class: center, middle ## `Retomando: ¿la puntuación de la prueba (testScore) se ve afectada por la longitud del cuerpo (bodyLength)?` --- class: center, middle ## `A partir del análisis exploratorio de datos se propone un modelo lineal` --- class: Left, middle ## Usaremos {stats} 📦 `>` lm(): Ajusta modelos lineales `>` ¡existen muchas funciones más! --- class: Left, middle ## Usaremos {broom} <img src="imgs/hex_broom.png" alt="Sharingan" width="6%" align="center"/> `>` tidy(): Construye un tibble que resume los hallazgos estadísticos del modelo `>` augment(): Agrega columnas a los datos originales que se modelaron `>` glance(): Construye un resumen de otras informaciones relacionadas al modelo --- class: left, middle ## Cargar paquetes .left[ ```r library(tidymodels) library(broom) ``` ] --- class: left, middle .left[ ## Ajustar un modelo lineal ```r ajuste_lm <- lm(testScore ~ bodyLength, data=dragones) ``` ] --- class: left, middle ## Obtener información sobre el modelo .left[ ```r tidy(ajuste_lm) %>% tibble::as_tibble() ``` ``` ## # A tibble: 2 x 5 ## term estimate std.error statistic p.value ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 (Intercept) -61.3 12.1 -5.08 5.38e- 7 ## 2 bodyLength 0.555 0.0597 9.29 5.59e-19 ``` ```r glance(ajuste_lm) %>% tibble::as_tibble() ``` ``` ## # A tibble: 1 x 12 ## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 0.153 0.151 21.2 86.2 5.59e-19 1 -2146. 4298. 4311. ## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int> ``` ] --- class: left, middle ## Graficar el ajuste del modelo .pull-left[ ```r info_ajuste_lm<- augment_columns(ajuste_lm, dragones) info_ajuste_lm %>% ggplot(aes(x=bodyLength,y=testScore))+ geom_jitter(alpha=.2)+ geom_line(aes(x=bodyLength,y=.fitted))+ theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ❗ ### Linealidad .pull-left[ ```r plot(ajuste_lm, which=1) ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ❗ ### Normalidad .pull-left[ ```r plot(ajuste_lm, which=2) ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ❗ ### Homocedasticidad .pull-left[ ```r plot(ajuste_lm, which=3) ``` ] .pull-right[ <!-- --> ] --- class: center, middle ## `¿Qué hacemos si el modelo no es adecuado?`  --- class: center, middle ## `¡Dejémos que los datos hablen!` 🗣 --- class: left, middle ## Análisis exploratorio de datos 🕵 🐉 .pull-left[ ```r dragones %>% ggplot(aes(x=bodyLength,y=testScore)) + geom_jitter(alpha=2) + theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Análisis exploratorio de datos 🕵 🐉 ⛰ .pull-left[ ```r dragones %>% ggplot(aes(x=bodyLength,y=testScore, colour=mountainRange)) + geom_jitter(alpha=2) + theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Análisis exploratorio de datos 🕵 🐉 ⛰ .pull-left[ ```r dragones %>% ggplot(aes(bodyLength,testScore, colour = mountainRange))+ geom_jitter(alpha=2) + facet_wrap(~ mountainRange) + theme_bw()+ theme(strip.background =element_rect(fill="white")) ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Análisis exploratorio de datos 🕵 🐉 ⛰ .pull-left[ ```r dragones %>% ggplot(aes(x=mountainRange, y=testScore, colour=mountainRange))+ geom_boxplot(alpha=.5)+ coord_flip()+ theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: center, middle ## `¡No podemos ignorar las cadenas montañosas!` ⛰ --- class: center, middle ## `¿Cómo las incorporamos al modelo?` ⛰ --- class: center, middle ## `¡Con los modelos lineales mixtos!`  --- class: center, middle ##`A partir del análisis exploratorio de datos se propone un modelo lineal mixto` --- class: Left, middle ## `Usaremos {lme4}` 📦 `>` lmer(): Ajusta un modelo lineal mixto `>` ¡existen muchas funciones más! --- class: Left, middle ## `Usaremos {broom.mixed}` <img src="imgs/hex_broom.png" alt="Sharingan" width="6%" align="center"/> `>` tidy(): Construye un tibble que resume los hallazgos estadísticos del modelo `>` augment(): Agrega columnas a los datos originales que se modelaron `>` glance(): Construye un resumen de otras informaciones relacionadas al modelo --- class: left, middle ## Cargar paquetes .left[ ```r library(broom.mixed) library(lme4) ``` ] --- class: left, middle ## Ajustar un modelo lineal mixto .left[ ### `Incorporaremos las cadenas montañosas con la expresión (1|mountainRange)` ```r ajuste_lmer <- lmer(testScore ~ bodyLength + (1|mountainRange), data = dragones) ``` ] --- class: left, middle ## Obtener información sobre el modelo .left[ ```r tidy(ajuste_lmer) %>% tibble::as_tibble() ``` ``` ## # A tibble: 4 x 6 ## effect group term estimate std.error statistic ## <chr> <chr> <chr> <dbl> <dbl> <dbl> ## 1 fixed <NA> (Intercept) 43.7 17.1 2.55 ## 2 fixed <NA> bodyLength 0.0332 0.0786 0.422 ## 3 ran_pars mountainRange sd__(Intercept) 18.4 NA NA ## 4 ran_pars Residual sd__Observation 15.0 NA NA ``` ```r glance(ajuste_lmer) %>% tibble::as_tibble() ``` ``` ## # A tibble: 1 x 6 ## sigma logLik AIC BIC REMLcrit df.residual ## <dbl> <dbl> <dbl> <dbl> <dbl> <int> ## 1 15.0 -1996. 3999. 4016. 3991. 476 ``` ] --- class: center, middle ## `¿No hay p-valor?`  --- class: left, middle ## Adicionando p-valor ✊🏻 .left[ ```r library(lmerTest) ajuste_lmer <- lmer(testScore ~ bodyLength + (1|mountainRange), data = dragones) ``` ```r tidy(ajuste_lmer) %>% tibble::as_tibble() ``` ``` ## # A tibble: 4 x 8 ## effect group term estimate std.error statistic df p.value ## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 fixed <NA> (Intercept) 43.7 17.1 2.55 177. 0.0116 ## 2 fixed <NA> bodyLength 0.0332 0.0786 0.422 473. 0.673 ## 3 ran_pars mountainRan… sd__(Interce… 18.4 NA NA NA NA ## 4 ran_pars Residual sd__Observat… 15.0 NA NA NA NA ``` ] --- class: left, middle ## Graficar el ajuste del modelo .pull-left[ ```r info_ajuste_lmer<- augment_columns(ajuste_lmer,dragones) info_ajuste_lmer %>% ggplot(aes(x=bodyLength,y=testScore, colour=mountainRange))+ geom_jitter(alpha=2)+ facet_wrap(~ mountainRange)+ geom_line(aes(x=bodyLength,y=.fitted), colour="black")+ theme_bw() ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ### Linealidad y homocedasticidad .pull-left[ ```r y.fit <- fitted(ajuste_lmer) res.fit <- residuals(ajuste_lmer) plot(y.fit, res.fit) abline(h=0, lty=2,col="red") ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ### Normalidad del error .pull-left[ ```r y.fit <- fitted(ajuste_lmer) res.fit <- residuals(ajuste_lmer) qqnorm(res.fit) qqline(res.fit) ``` ] .pull-right[ <!-- --> ] --- class: left, middle ## Verificar los supuestos ### Normalidad del efecto aleatorio .pull-left[ ```r pred.fit<-ranef(ajuste_lmer)[[1]][[1]] qqnorm(pred.fit) qqline(pred.fit) ``` ] .pull-right[ <!-- --> ] --- class: center, middle ##`¿Y si ahora incorporamos los sitios?` 👀 .left[ `Podemos incorporar los sitios con la expresión (1|mountainRange/site) o (1|mountainRange) + (1|mountainRange:site)` ] --- class: left ##`Referencias` 📖 ### > Gabriela Hajduk https://gkhajduk.github.io/2017-03-09-mixed-models/ ### > Athanasia Mowinckel https://athanasiamo.github.io/LME_introduction_workshop/ ### > Ben Bolker https://github.com/bbolker ### > Douglas Bates https://github.com/dmbates --- class: left ##`Referencias`📖 ### > Max Kuhn, Hadley Wickham and RStudio https://www.tidymodels.org/ ### > Max Kuhn and Julia Silge https://www.tmwr.org/ ### > David Robinson y + https://github.com/tidymodels/broom ### > Ben Bolker y + https://github.com/bbolker/broom.mixed --- class: center, middle ##`¡Comienza a usar modelos mixtos y {tidymodels} con tus propios datos!`👩🏻‍💻 --- class: center, middle .center[<img src=imgs/R-Ladies_Talca_hex.png width="30%">] ##### Presentación creada con el paquete [**xaringan**](https://github.com/yihui/xaringan) de [**Yihui Xie**](https://github.com/yihui) y el tema [**rladies**](https://github.com/rbind/apreshill/blob/master/static/slides/rladies-demo-slides.Rmd) de [**Alison Hill**](https://github.com/apreshill) <img src="imgs/logo_twitter.png" alt="Sharingan" width="6%" align="center"/> `@aleants @RLadiesTalca`