Un poco de EDA y modelos (WIP)

---

1. Introducción

Hemos visto en tutoriales anteriores como cargar, manejar, arreglar y visualizar datos à la tidyverse. Una vez somos capaces de manejar y arreglar nuestros datos con R podemos pasar a intentar efectuar un análisis “realista” con un conjunto de datos.

En este tutorial y siguientes, vamos a utilizar un conjunto de datos sobre nacimientos de bebes en España para mostrar algunas técnicas de EDA y modelización. No se mostrarán los detalles técnicos de las técnicas, sino solo su aplicación práctica y, a veces, la interpretación de los resultados.

Antes de proceder a la estimación de modelos estadísticos formales o contrastar hipótesis se suele realizar una etapa conocida como análisis exploratorio. El análisis exploratorio (EDA) es una parte importante de todo análisis de datos. Suele tener lugar antes de la especificación y estimación de modelos estadísticos formales. Su objetivo último es “comprender” los datos y las relaciones existentes entre las variables.

Una de las máximas de un data scientist es, como señalan aquí, “Know Your Data”. Además nos explican que:

Data exploration is the art of looking at your data, rapidly generating hypotheses, quickly testing them, then repeating again and again and again. The goal of data exploration is to generate many promising leads that you can later explore in more depth.

With exploratory data analysis, you’ll combine visualisation and transformation with your curiosity and scepticism to ask and answer interesting questions about data

Un buen ejemplo de qué es un análisis exploratorio son los “screencast” de David Robinson: cada semana graba un vídeo con un análisis rápido de un conjunto de datos. Puedes verlos, en su canal de youtube ,aquí. En este repo de Github tienes la transcripción del código, en ficheros .Rmd que acaba utilizando en los screencasts. Como hay mucha gente que sigue sus vídeos hay, como explican en este post, varias personas que mantienen un listado diseccionando las tareas que hace David en cada screencast, puedes encontrarlos aquí y aquí.

Por último, un libro, bueno un bookdown, sobre EDA: Exploratory Data Analysis with R de Roger Peng.

¿Qué queremos saber?

Generalmente un estudio cuantitativo comienza por una pregunta(s) que guía el análisis. En nuestro caso, el objetivo del tutorial es meramente mostrar algunas técnicas/funciones útiles en la fase de exploración de datos (EDA) de un estudio, aún así, nuestro “estudio” estará guiado por varias preguntas relacionadas con los bebes; por ejemplo: 1) vamos a intentar ver/contestar si los bebitos nacen con más peso que las bebitas 2) intentaremos contestar a la pregunta(s) de si determinadas características de las madres/padres, como la edad o el nivel educativo, afectan al peso de los recién nacidos.

El siguiente paso de todo análisis empírico es la búsqueda de un conjunto de datos que permita (hasta cierto grado de confianza) dar una respuesta a las preguntas planteadas. ¿De donde sacamos esos datos? En nuestro caso, utilizaremos la estadística de nacimientos del INE. Los detalles en la sección siguiente.

---

2. Los datos

Origen de los datos

Los datos provienen del INE, concretamente de la Estadística de nacimientos. Como puedes ver en la web del INE hay una serie de tablas que muestran resultados para las principales variables; por ejemplo, hay una tabla llamada “Nacimientos por edad de la madre, mes y sexo”. Estas tablas están muy bien para mostrar los principales resultados de la Estadística de nacimientos, pero nosotros queremos analizar “de verdad” los datos, así que tendremos que descargarnos los microdatos de la encuesta.

En los microdatos de la Estadística de Nacimientos hay datos sobre: (i) Nacimientos, (ii) Muertes fetales tardías y (iii) Partos. Vamos a usar los ficheros de datos sobre PARTOS.

Hay un fichero para cada año desde 1975, pero como han habido diversos cambios en la metodología de la estadística, solo vamos a trabajar los ficheros de los años 2007 a 2018 que comparten diccionario.

Me descargue (a mano) los ficheros de microdatos de nacimientos para los años 2007 a 2018. Los datos están en formato texto, y cada registro es una cadena larga de caracteres. Los datos van acompañados de un diccionario que sirve para poder separar las cadenas de caracteres en los valores de cada variable. Os ahorro los detalles del proceso.

Los datos eran de partos, pero los preparé para que cada fila pertenezca a los datos de un solo bebe: cada bebe tiene su propia fila.

Cargando el fichero de datos y dicc

Mi ordenador tuvo problemas para mover el fichero con todos los datos (2007 a 2018), y, como además, el objetivo del tutorial es tan sólo mostrar algunas técnicas y funciones útiles para EDA, por lo que solamente se utilizarán los datos referentes a 2017. Adicionalmente, se restringió la muestra a partos con un solo bebe que habían tenido lugar en un centro sanitario.

df <- rio::import(here::here("datos", "df_bebes_2017.rds"))             #- 374.933 x 13
df_dicc <- rio::import(here::here("datos", "df_bebes_2017_dicc.xlsx"))  #- 13 x 18

Con el anterior chunk hemos cargado el fichero de datos df, que cuenta con 374.933 registros de bebes y 13 columnas o variables para cada bebe.

---

3. Cosas de EDA

Hemos visto en tutoriales anteriores como cargar, arreglar y manejar datos con R à la tidyverse. En este tutorial vamos a ver algunos paquetes y funciones que conviene conocer porque nos pueden ayudar a acelerar el análisis exploratorio inicial de nuestros datos.

En este repo de Github, Mateusz Staniak mantiene un listado de paquetes y recursos relacionados con la EDA.

En este post y en este artículo nos hablan del paquete autoEDA que trata de automatizar el proceso inicial exploratorio de datos (EDA). Además el post nos muestra los 11 paquetes de R, relacionados con la EDA, más descargados de CRAN. El objetivo de estos 11/12 paquetes es similar, tratar de automatizar/facilitar el análisis preliminar o exploratorio de los datos (EDA). Evidentemente automatizar totalmente la EDA es muy-muy complicado, por no decir imposible, pero al menos estos paquetes (y otros muchos) contienen funciones que nos pueden ayudar en las primeras etapas de nuestros análisis de datos.

Por ejemplo, en este post, Laura Ellis, nos explica como hacer EDA en R mostrándonos algunas de sus funciones/trucos favoritos como la función skimr::skim(), visdat::vis_dat() o DataExplorer::create_report(). El post tuvo una segunda parte (aquí), donde Laura nos explica el uso de algunas funciones útiles para EDA del paquete inspectdf como por ejemplo: inspectdf::inspect_types(), inspectdf::inspect_na, etc…

En este otro post, Sharla Gelfand nos habla sobre estrategias para afrontar la EDA cuando se trabaja con un conjunto de datos nuevo. Sharla utiliza, entre otros, los paquetes visdat, skimr y assertr.

DE entre los muchos paquetes y funciones relacionados con EDA, conviene conocer, en mi opinión, estos¹:

3.1 Manipular los nombres de las variables

Con los datos de bebes no vamos a cambiar los nombres, ya que son los nombres que les dio el INE y que figuran en el diccionario oficial, pero muchas veces hay que cambiar/arreglar los nombres de las variables, así que vamos a hacerlo con nuestro df, aunque al final volveremos a utilizar los nombres originales.

Podemos acceder a los nombres de las columnas de un data.frame con la función base::names(). Veámoslo:

nombres_originales <- names(df)  #- almacenamos los nombres originales en el vector nombres_originales

Veamos cuales son los nombres de las variables de df:

names(df)
#>  [1] "MESPAR"     "MUNPAR"     "MUNREM"     "NORMA.f"    "CESAREA.f" 
#>  [6] "SEXO1"      "SEMANAS"    "PESON1"     "EDADM"      "EDADP"     
#> [11] "ESTUDIOM.f" "ESTUDIOP.f" "PAISNXM"

Si quisieramos, por ejemplo, cambiar el nombre de la segunda variable, podemos hacerlo con:

names(df)[2] <- "nuevo_nombre_variable_2"
names(df)
#>  [1] "MESPAR"                  "nuevo_nombre_variable_2"
#>  [3] "MUNREM"                  "NORMA.f"                
#>  [5] "CESAREA.f"               "SEXO1"                  
#>  [7] "SEMANAS"                 "PESON1"                 
#>  [9] "EDADM"                   "EDADP"                  
#> [11] "ESTUDIOM.f"              "ESTUDIOP.f"             
#> [13] "PAISNXM"

Si tuviésemos un data.frame con nombres realmente extraños, incluso con nombres “no sintácticos”, y necesitamos arreglarlos de forma rápida, podemos hacerlo con la función janitor::clean_names(). Arregla los nombres de las variables de forma automática, además tiene algunas opciones para hacerlo como más te guste.

janitor::clean_names(df) %>% names()
#>  [1] "mespar"                  "nuevo_nombre_variable_2"
#>  [3] "munrem"                  "norma_f"                
#>  [5] "cesarea_f"               "sexo1"                  
#>  [7] "semanas"                 "peson1"                 
#>  [9] "edadm"                   "edadp"                  
#> [11] "estudiom_f"              "estudiop_f"             
#> [13] "paisnxm"

Volvamos a dejar los nombres originales de df:

names(df) <- nombres_originales

---

3.2 Estructura y tipo de variables

Siempre-siempre hay que saber de que tipo son nuestras variables. Esto se puede hacer de muchas formas. Os muestro dos:

a) con str()

str(df)            #- str() muestra la estructura interna de un objeto R
#> tibble [374,933 × 13] (S3: tbl_df/tbl/data.frame)
#>  $ MESPAR    : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ MUNPAR    : chr [1:374933] "28079" "26036" "14021" "47186" ...
#>  $ MUNREM    : chr [1:374933] NA "26018" "14021" "47186" ...
#>  $ NORMA.f   : Factor w/ 2 levels "Parto con complicaciones",..: 1 2 1 1 1 1 2 1 1 1 ...
#>  $ CESAREA.f : Factor w/ 3 levels "Con cesárea",..: 1 2 1 1 1 1 2 1 1 1 ...
#>  $ SEXO1     : Factor w/ 2 levels "bebita","bebito": 2 1 2 2 1 1 2 2 2 2 ...
#>  $ SEMANAS   : num [1:374933] 26 26 26 25 26 26 23 27 NA 28 ...
#>  $ PESON1    : num [1:374933] 420 500 540 560 592 600 600 600 635 650 ...
#>  $ EDADM     : num [1:374933] 39 30 31 36 42 41 36 32 34 40 ...
#>  $ EDADP     : num [1:374933] NA 32 31 36 42 44 39 30 26 47 ...
#>  $ ESTUDIOM.f: Factor w/ 3 levels "Primarios","Medios",..: NA 1 3 1 3 1 NA 1 1 NA ...
#>  $ ESTUDIOP.f: Factor w/ 3 levels "Primarios","Medios",..: NA 2 3 1 3 1 NA 1 1 NA ...
#>  $ PAISNXM   : chr [1:374933] "122" "228" "108" "407" ...

b) con inspectdf::inspect_types()

inspectdf::inspect_types(df)   #- muestra de que tipo son las variables

type	cnt	pcnt	col_name
factor	6	46.15385	MESPAR , NORMA.f , CESAREA.f , SEXO1 , ESTUDIOM.f, ESTUDIOP.f
numeric	4	30.76923	SEMANAS, PESON1 , EDADM , EDADP
character	3	23.07692	MUNPAR , MUNREM , PAISNXM

c) 2 funciones útiles del paquete pjpv2020.01

La función pjpv2020.01::pjp_f_estadisticos_basicos(). Sí, el nombre de la función es un poco largo, pero es que es un paquete para uso personal. Yo estoy acostumbrado a usarla. Da un resumen rápido de las variables del data.frame:

zz <- pjpv2020.01::pjp_f_estadisticos_basicos(df) 
gt::gt(zz)

variable	q_zeros	p_zeros	q_na	p_na	type	unique	min	max	mean	sd	NN	NN_ok
MESPAR	0	0	0	0.00	factor	12	NA	NA	NA	NA	374933	374933
MUNPAR	0	0	7821	2.09	character	613	NA	NA	NA	NA	374933	367112
MUNREM	0	0	68464	18.26	character	750	NA	NA	NA	NA	374933	306469
NORMA.f	0	0	0	0.00	factor	2	NA	NA	NA	NA	374933	374933
CESAREA.f	0	0	0	0.00	factor	2	NA	NA	NA	NA	374933	374933
SEXO1	0	0	0	0.00	factor	2	NA	NA	NA	NA	374933	374933
SEMANAS	0	0	42842	11.43	numeric	27	20	46	39.12	1.76	374933	332091
PESON1	0	0	17290	4.61	numeric	2889	300	6350	3248.34	506.24	374933	357643
EDADM	0	0	0	0.00	numeric	46	12	58	32.47	5.60	374933	374933
EDADP	0	0	11665	3.11	numeric	64	14	77	35.44	6.10	374933	363268
ESTUDIOM.f	0	0	137206	36.59	factor	3	NA	NA	NA	NA	374933	237727
ESTUDIOP.f	0	0	134748	35.94	factor	3	NA	NA	NA	NA	374933	240185
PAISNXM	0	0	2971	0.79	character	64	NA	NA	NA	NA	374933	371962

Muchas veces también me es muy útil la función pjpv2020.01::pjp_f_unique_values() porque devuelve un df con los valores únicos de cada variable:

zz <- pjpv2020.01::pjp_f_unique_values(df, truncate = TRUE, nn_truncate = 47) 
gt::gt(zz)

variables	nn_unique	unique_values
MESPAR	12	1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 10 - 11 - 1
MUNPAR	614	NA - 01059 - 02003 - 02009 - 02025 - 02037 - 02
MUNREM	751	NA - 01002 - 01036 - 01059 - 02003 - 02009 - 02
NORMA.f	2	Parto con complicaciones - Parto normal
CESAREA.f	2	Con cesárea - Sin cesárea
SEXO1	2	bebita - bebito
SEMANAS	28	NA - 20 - 21 - 22 - 23 - 24 - 25 - 26 - 27 - 28
PESON1	2890	NA - 300 - 410 - 420 - 430 - 440 - 450 - 465 -
EDADM	46	12 - 13 - 14 - 15 - 16 - 17 - 18 - 19 - 20 - 21
EDADP	65	NA - 14 - 15 - 16 - 17 - 18 - 19 - 20 - 21 - 22
ESTUDIOM.f	4	NA - Primarios - Medios - Universidad
ESTUDIOP.f	4	NA - Primarios - Medios - Universidad
PAISNXM	65	NA - 102 - 103 - 104 - 107 - 108 - 110 - 112 -

• La función dplyr::distinct() es muy útil para ver los con valores únicos de una o varias variables. Por ejemplo:

zz <- df %>% distinct(NORMA.f, CESAREA.f)
zz %>% gt::gt()

NORMA.f	CESAREA.f
Parto con complicaciones	Con cesárea
Parto normal	Sin cesárea
Parto normal	Con cesárea
Parto con complicaciones	Sin cesárea

---

3.3 Estadísticos básicos de las variables

a) con pjpv2020.01

Ya he dicho que yo suelo utilizar estas 2 funciones:

zz <- pjpv2020.01::pjp_f_estadisticos_basicos(df)   #- estadísticos básicos del df
zz <- pjpv2020.01::pjp_f_unique_values(df)          #- valores únicos de cada variable de df

b) con summarytools

La función summarytools::dfSummary() da un informe/resumen de las variables de un data.frame muy útil. Es muy útil para hacerse una idea rápida de las propiedades de las variables, pero tiene la pega de que el output es muy voluminosos para mostrarlo aquí, así que si quieres verlo tendrás que ejecutar el anterior chunk de forma interactiva.

zz <- summarytools::dfSummary(df) #- genera un fichero con un resumen útil y agradable de ver 
summarytools::view(zz)  #- para visualizar el informe

c) con skimr

Otra alternativa es usar `skimr::skim(df). Ocurre lo mismo: como su output es voluminoso, no lo voy a mostrar aquí.

skimr::skim(df)

---

3.4 NA’s

Siempre hay que chequear la existencia de NA’s en nuestro df. NA representa una característica que existe pero no sabemos su valor En R, los valores no disponibles se representan con NA. Se puede chequear si un valor es NA con la función is.na(). En datos generados por otros paquetes estadísticos, los valores ausentes pueden estar codificados con diferentes valores como por ejemplo “-999”, “N/A”, un espacio en blanco, etc …

Tres paquetes que nos pueden ayudar en esta tarea son DataExplorer, visdat y naniar.

DataExplorer::plot_missing(df)

naniar::gg_miss_var(df, show_pct = TRUE)        

El paquete visdat es útil para ver los NA’s pero el siguiente chunk no nos va a funcionar para las 374.933 filas que tiene nuestro df.

visdat::vis_dat(df) #- no funciona , hay demasiadas filas

Pero veamos como funcionaría para, por ejemplo, las primeras 1.000 filas de df:

df %>% slice(1:1000) %>% visdat::vis_dat()

• Algo parecido pero con visdat::vis_miss()

df %>% slice(1:1000) %>% visdat::vis_miss(cluster = TRUE) 

• Otra forma de visualizar la co-ocurrencia de NA’s en las variables de df:

naniar::gg_miss_upset(df)

trabajando con los NA's

• Por ejemplo, podemos querer seleccionar las variables que tienen NA’s:

zz_con_NAs  <- df %>% select(where(anyNA))

• A lo mejor, nos puede interesar ver la ocurrencia de NA’s para los grupos/categorías de una variable categórica de df, por ejemplo los estudios de la madre(ESTUDIOM.f)

naniar::gg_miss_var(df, facet = ESTUDIOM.f, show_pct = TRUE) #- faceted por la variable ESTUDIOM.f

¿Que hacemos con los NA’s?

Pues no está claro, depende … pero supongamos que queremos dejar nuestro df SIN NA’s;

zz <- df %>% tidyr::drop_na()              #- con tidyr

#- otras formas de hacer lo mismo ---
zz <- df[complete.cases(df), ]             #- con base-R 
zz <- df %>% filter(complete.cases(.))     #- con dplyr y base-R

Si quitáramos todos los registros/bebes que tienen algún NA en alguna de las 13 variables, nos quedaríamos con 152.286 registros.

Se puede refinar la eliminación de NA's. Puedes ser que igual Igual nos interesase eliminar sólo las filas que tengan un NA en una determinada variable, por ejemplo, si quisieramos quitar las filas que tuviesen NA en la variable PESO1 y/o en la variable SEMANAS, lo haríamos así:

zz <- df %>% tidyr::drop_na(c(PESON1, SEMANAS))   #- quito filas con NA en PESON1 o en SEMANAS

• Imagina que hubiese una fila o una columna con todos sus valores NA. Lógicamente habría que eliminarlas. Se puede hacer fácil con `janitor::remove_empty()

zz <- df %>% janitor::remove_empty(c("rows", "cols"))    #- quita variables y filas vacías

---

br>

3.6 Variables numéricas

Vamos a ver algunas funciones útiles para hacer un análisis inicial de las variables numéricas.

Si quisieramos crear un df sólo con las variables numéricas:

df_numeric <- df %>% select_if(is.numeric)

a) Estadísticos descriptivos

Por ejemplo, summarytools::descr() nos ayuda a calcular rápidamente estadísticos descriptivos de las variables numéricas.

summarytools::descr(df, style = "rmarkdown")

summarytools::descr(df)
#> Descriptive Statistics  
#> df  
#> N: 374933  
#> 
#>                         EDADM       EDADP      PESON1     SEMANAS
#> ----------------- ----------- ----------- ----------- -----------
#>              Mean       32.47       35.44     3248.34       39.12
#>           Std.Dev        5.60        6.10      506.24        1.76
#>               Min       12.00       14.00      300.00       20.00
#>                Q1       29.00       32.00     2970.00       38.00
#>            Median       33.00       35.00     3260.00       39.00
#>                Q3       36.00       39.00     3570.00       40.00
#>               Max       58.00       77.00     6350.00       46.00
#>               MAD        5.93        5.93      444.78        1.48
#>               IQR        7.00        7.00      600.00        2.00
#>                CV        0.17        0.17        0.16        0.04
#>          Skewness       -0.41        0.18       -0.56       -2.30
#>       SE.Skewness        0.00        0.00        0.00        0.00
#>          Kurtosis        0.03        0.97        2.12       11.20
#>           N.Valid   374933.00   363268.00   357643.00   332091.00
#>         Pct.Valid      100.00       96.89       95.39       88.57

b) Histogramas o f. de densidad

En general, los métodos estadísticos requieren que las variables sigan una determinada distribución. Esto puede chequearse formalmente o utilizar histogramas y/o funciones de densidad estimadas.

Podemos usar el paquete DataExplorer para obtener histogramas o gráficos de densidad para las variables numéricas.

• Histogramas:

DataExplorer::plot_histogram(df, ncol = 2)

• Funciones de densidad estimadas:

DataExplorer::plot_density(df, ncol = 2)

c) Matriz de correlaciones

Otro aspecto importante de un análisis exploratorio consiste en analizar la existencia de relaciones entre las variables del conjunto de datos. Esto puede hacerse de diversas maneras, dos de las más sencillas y usadas son los gráficos de dispersión y las matrices de correlación.

Podemos obtener la matriz de correlaciones de diversas formas:

• Primero con la función cor() de R-base:

stats::cor(df_numeric, use = "pairwise.complete.obs") %>%  #- devuelve una matriz, no un df
      round(. , 2)
#>         SEMANAS PESON1 EDADM EDADP
#> SEMANAS    1.00   0.55 -0.02 -0.01
#> PESON1     0.55   1.00  0.01  0.03
#> EDADM     -0.02   0.01  1.00  0.65
#> EDADP     -0.01   0.03  0.65  1.00

#df_numeric %>% GGally::ggcorr(label = TRUE)
df_numeric %>% GGally::ggpairs()

• Ahora con la función corrr::correlate()

df %>% select_if(is.numeric) %>% 
       corrr::correlate() %>% 
       pjpv2020.01::pjp_f_decimales(nn = 2) %>%  
       gt::gt()

term	SEMANAS	PESON1	EDADM	EDADP
SEMANAS	NA	0.55	-0.02	-0.01
PESON1	0.55	NA	0.01	0.03
EDADM	-0.02	0.01	NA	0.65
EDADP	-0.01	0.03	0.65	NA

• Otra forma más de visualizar las correlaciones, con inspectdf::inspect_cor()

df %>% inspectdf::inspect_cor()

col_1	col_2	corr	p_value	lower	upper	pcnt_nna
EDADP	EDADM	0.6466554	0.0000000	0.6447594	0.6485435	96.88878
PESON1	SEMANAS	0.5538819	0.0000000	0.5514923	0.5562625	86.54826
EDADP	PESON1	0.0286436	0.0000000	0.0253196	0.0319670	92.59441
EDADM	SEMANAS	-0.0224419	0.0000000	-0.0258410	-0.0190422	88.57343
EDADP	SEMANAS	-0.0094270	0.0000001	-0.0128786	-0.0059752	85.98070
EDADM	PESON1	0.0062598	0.0001814	0.0029825	0.0095370	95.38851

• Además, después, con inspectdf::show_plot() se pueden visualizar las correlaciones en un gráfico:

df %>% inspectdf::inspect_cor() %>% inspectdf::show_plot()

• El paquete correlation facilita el cálculo de diferentes estadísticos de correlación, por defecto calcula el coeficiente de correlación de Pearson:

correlation::correlation(df_numeric)    #- remotes::install_github("easystats/correlation")

Parameter1	Parameter2	r	CI_low	CI_high	t	df	p	Method	n_Obs
SEMANAS	PESON1	0.5538819	0.5514923	0.5562625	378.955289	324496	0.0000000	Pearson	324498
SEMANAS	EDADM	-0.0224419	-0.0258410	-0.0190422	-12.935864	332089	0.0000000	Pearson	332091
SEMANAS	EDADP	-0.0094270	-0.0128786	-0.0059752	-5.352658	322368	0.0000002	Pearson	322370
PESON1	EDADM	0.0062598	0.0029825	0.0095370	3.743631	357641	0.0001814	Pearson	357643
PESON1	EDADP	0.0286436	0.0253196	0.0319670	16.883950	347165	0.0000000	Pearson	347167
EDADM	EDADP	0.6466554	0.6447594	0.6485435	510.957420	363266	0.0000000	Pearson	363268

d) Boxplots frente a una v. categórica

Para ver diferencias en la distribución de las variables numéricas frente a una variable categórica es muy útil DataExplorer::plot_boxplot(). Por ejemplo frente a la variable ESTUDIOM.f:

DataExplorer::plot_boxplot(df, by = "ESTUDIOM.f")

O, por ejemplo, frente a la variable CESAREA.f

my_vv <- names(df)[3]
my_vv <- "CESAREA.f"
DataExplorer::plot_boxplot(df, by = my_vv)

• También se pueden hacer boxplots de una variables numérica frente a 2 categóricas:

df %>% explore::explore(EDADM, ESTUDIOM.f, target = CESAREA.f)

• Con explore::explore_all() se muestran gráficos de todas las variables del df frente a una variable target u objetivo:

df %>% select(NORMA.f, EDADM, PESON1,  ESTUDIOM.f, CESAREA.f) %>% explore::explore_all(target = CESAREA.f)

Veamos la relación entre “Parto normal” y CESAREA.f:

df %>% janitor::tabyl(NORMA.f, CESAREA.f) %>% janitor::adorn_percentages()

NORMA.f	Con cesárea	Sin cesárea	Otros_xx
Parto con complicaciones	0.7292581	0.2707419	0
Parto normal	0.1796008	0.8203992	0

e) Scatterplots entre variables numéricas

DataExplorer::plot_scatterplot() nos permite hacer rápidamente un scatterplot de todas las variables del df frente a una variable, por ejemplo el peso del bebe: PESON1

my_vv <- "PESON1"
DataExplorer::plot_scatterplot(df, by = my_vv, sampled_rows = 500L)

---

3.7 Variables categóricas

Si necesitásemos un df con todas las variables categóricas, ¿cómo lo obtendrías?

df %>% select_if(is.numeric) %>% names()      #- old fashion
#> [1] "SEMANAS" "PESON1"  "EDADM"   "EDADP"
df %>% select(where(is.numeric)) %>% names()  
#> [1] "SEMANAS" "PESON1"  "EDADM"   "EDADP"

¿Y si quisieras seleccionar las variables no-numéricas?

df %>% select_if(!(is.numeric(.))) #- NO funciona
df %>% select_if(~ !is.numeric(.)) %>% names()     #- anonymous function but select_if
df %>% select(where(~ !is.numeric(.))) %>% names() #- anonymous functions & where
df %>% select(where(purrr::negate(is.numeric)))  %>% names()   #- purrr::negate()!!!!

Para visualizar rápidamente los valores de las variables categóricas, tenemos inspectdf::inspect_cat(), que nos devuelve un gráfico con la distribución de todas las variables categóricas.

inspectdf::inspect_cat(df) %>% inspectdf::show_plot(high_cardinality = 1)

• También podemos hacerlo con DataExplorer::plot_bar():

DataExplorer::plot_bar(df)

---

3.8 Paquetes con shiny

a) burro

El paquete burro:

burro attempts to make EDA accessible to a larger audience by exposing datasets as a simple Shiny App

#- devtools::install_github("laderast/burro")
df_small <- df %>% slice(1:1000)
burro::explore_data(df_small, outcome_var = colnames(df))

b) explore

El paquete explore.

Instead of learning a lot of R syntax before you can explore data, the explore package enables you to have instant success. You can start with just one function - explore() - and learn other R syntax later step by step

• Podemos crear un shiny para explorar nuestro df completo:

explore::explore(df)

• o solo explorar una o varias variables.

df %>% explore::explore(CESAREA.f)
df %>% explore::explore(EDADM, ESTUDIOM.f, target = CESAREA.f)

c) ExPanDaR

El paqueteExPanDaR también permite explorar los datos a través de un shiny. Explican su funcionamiento en este post y este otro.

library(ExPanDaR) #- remotes::install_github("joachim-gassen/ExPanDaR")

ExPanD(mtcars)

---

4. Tablas

Muchas veces hay que presentar resultados básicos como una tabla de casos o de frecuencias. para ello tenemos muchas posibilidades:

4.1 Tablas con summarytools

El paquete summarytools.

summarytools is an R package providing tools to neatly and quickly summarize data. It can also make R a little easier to learn and to use, especially for data cleaning and preliminary analysis.

• La función freq() provee tablas con conteos y frecuencias

summarytools::freq(df$SEXO1, style = "rmarkdown")

Frequencies

df$SEXO1

Type: Factor

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
bebita	181708	48.46	48.46	48.46	48.46
bebito	193225	51.54	100.00	51.54	100.00
<NA>	0			0.00	100.00
Total	374933	100.00	100.00	100.00	100.00

• tabulación cruzada entre dos variables categóricas:

summarytools::ctable(df$SEXO1, df$CESAREA.f)

Cross-Tabulation, Row Proportions
SEXO1 * CESAREA.f
Data Frame: df

	CESAREA.f	Con cesárea	Sin cesárea	Otros_xx	Total
SEXO1
bebita		43247 (23.8%)	138461 (76.2%)	0 (0.0%)	181708 (100.0%)
bebito		49763 (25.8%)	143462 (74.2%)	0 (0.0%)	193225 (100.0%)
Total		93010 (24.8%)	281923 (75.2%)	0 (0.0%)	374933 (100.0%)

• Incluso se pueden hacer test chi-cuadrado

summarytools::ctable(df$ESTUDIOM.f, df$NORMA.f, chisq = TRUE)

Cross-Tabulation, Row Proportions
ESTUDIOM.f * NORMA.f
Data Frame: df

	NORMA.f	Parto con complicaciones	Parto normal	Total
ESTUDIOM.f
Primarios		11898 (12.6%)	82201 (87.4%)	94099 (100.0%)
Medios		8394 (12.8%)	56953 (87.2%)	65347 (100.0%)
Universidad		8553 (10.9%)	69728 (89.1%)	78281 (100.0%)
		17860 (13.0%)	119346 (87.0%)	137206 (100.0%)
Total		46705 (12.5%)	328228 (87.5%)	374933 (100.0%)

Chi.squared	df	p.value
161.1	2	0

---

4.2 Tablas con janitor

La verdad es que janitor es un paquete fantástico!!!

janitor has simple functions for examining and cleaning dirty data. It was built with beginning and intermediate R users in mind and is optimized for user-friendliness. Advanced R users can already do everything covered here, but with janitor they can do it faster and save their thinking for the fun stuff.

• hacer (y dar formato) a tablas de 1, 2 o 3 variables:

df %>% janitor::tabyl(CESAREA.f) %>% gt::gt()

CESAREA.f	n	percent
Con cesárea	93010	0.248071
Sin cesárea	281923	0.751929
Otros_xx	0	0.000000

df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f) %>% gt::gt()

CESAREA.f	Primarios	Medios	Universidad	NA_
Con cesárea	22486	17446	20401	32677
Sin cesárea	71613	47901	57880	104529
Otros_xx	0	0	0	0

df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f, SEXO1)  #- xq no podemos usar gt::gt()
#> $bebita
#>    CESAREA.f Primarios Medios Universidad   NA_
#>  Con cesárea     10497   8059        9489 15202
#>  Sin cesárea     35456  23486       28150 51369
#>     Otros_xx         0      0           0     0
#> 
#> $bebito
#>    CESAREA.f Primarios Medios Universidad   NA_
#>  Con cesárea     11989   9387       10912 17475
#>  Sin cesárea     36157  24415       29730 53160
#>     Otros_xx         0      0           0     0

---

4.3 Tablas con R-base

Las tablas con R-base tampoco están mal. El problema que tienen es que los resultados no se almacenan en data.frames, si no en matrices.

table(df$CESAREA.f, df$ESTUDIOM.f, useNA = "always") 
#>              
#>               Primarios Medios Universidad   <NA>
#>   Con cesárea     22486  17446       20401  32677
#>   Sin cesárea     71613  47901       57880 104529
#>   Otros_xx            0      0           0      0
#>   <NA>                0      0           0      0

• Como se almacenan en matrices, para poder graficarlas con gt primero las hemos de convertir a data.frame y casi seguro que habrá que hacer uso de pivot_wider()

my_tabla <- table(df$CESAREA.f, df$ESTUDIOM.f) 
my_tabla %>% as.data.frame() %>% 
            pivot_wider(names_from = 2, values_from = 3) %>% 
            gt::gt()

Var1	Primarios	Medios	Universidad
Con cesárea	22486	17446	20401
Sin cesárea	71613	47901	57880
Otros_xx	0	0	0

• mosaicplot(), un gráfico que me gusta del sistema de R-base

zz <- prop.table(my_tabla, 1)
zz <- zz[1:2,]   #- tengo que hacer esto xq hay NA's en la tabla
mosaicplot(t(zz), color = TRUE, main = "% de Cesáreas para niveles educativos de la madre")

---

5. Contrastes

Evidentemente R permite implementar múltiples técnicas y modelos estadísticos, así como hacer múltiples y variados contrastes de hipótesis. En esta pagina web tenéis un buena introducción a estos tema. Veamos algún ejemplo:

• t-test: https://statistics.berkeley.edu/computing/r-t-tests

t.test(df$PESON1, mu = 3250)
#> 
#>  One Sample t-test
#> 
#> data:  df$PESON1
#> t = -1.9636, df = 357642, p-value = 0.04958
#> alternative hypothesis: true mean is not equal to 3250
#> 95 percent confidence interval:
#>  3246.679 3249.997
#> sample estimates:
#> mean of x 
#>  3248.338
t.test(df$PESON1 ~ df$CESAREA.f)
#> 
#>  Welch Two Sample t-test
#> 
#> data:  df$PESON1 by df$CESAREA.f
#> t = -21.261, df = 125952, p-value < 0.00000000000000022
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>  -51.55217 -42.84971
#> sample estimates:
#> mean in group Con cesárea mean in group Sin cesárea 
#>                  3212.853                  3260.053

• correlación entre dos variables cuantitativas

library(Hmisc)
Hmisc::rcorr(as.matrix(df_numeric)) 

• un ejemplo para ver si hay diferencias en el peso en función de los estudios de la madre (!!!!)

library(purrr)
library(broom)
df %>% group_by(ESTUDIOM.f) %>% 
  summarise(t_test = list(t.test(PESON1))) %>% 
  mutate(tidied = map(t_test, tidy)) %>% 
  tidyr::unnest(tidied) %>% 
  ggplot(aes(estimate, ESTUDIOM.f)) + geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) + 
  labs(title = "Peso del bebe para diferentes niveles educativos de la madre")

• chi-squared test: https://data-flair.training/blogs/chi-square-test-in-r/

chisq.test(df$CESAREA.f, df$SEXO1)
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  df$CESAREA.f and df$SEXO1
#> X-squared = 191.51, df = 1, p-value < 0.00000000000000022
chisq.test(df$CESAREA.f, df$ESTUDIOM.f) 
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  df$CESAREA.f and df$ESTUDIOM.f
#> X-squared = 188.48, df = 2, p-value < 0.00000000000000022

• Aquí puedes encontrar un curso sobre como implementar los contrastes estadísticos más habituales con R.

En este post, Jonas Kristoffer Lindeløv, nos presenta un cuadro con los contrastes de hipótesis más frecuentes y cómo efectuar esos contrastes en el contexto de los modelos de regresión con R.

---

6. Modelos

En este apartado vamos a ver (un poco²) como estimar modelos lineales y no lineales con R. Para una introducción a la estimación de modelos en R puedes ir aquí.

Además, tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo de Machine Learning.

Antes restrinjamos aún un poco más el df con los datos sobre bebes:

df_m <- df %>% select(PESON1, SEMANAS, SEXO1, CESAREA.f, EDADM, EDADP, ESTUDIOM.f, ESTUDIOP.f)
df_m <- df_m %>% drop_na()

6.1 Modelos lineales

La función para estimar modelos lineales es lm(), así que vamos a utilizarla para estimar nuestro primer modelo con R. Estimaremos un modelo lineal con variable a explicar el peso del bebe (PESON1) en función de todas las demás variables en df_m.

mod_1 <- lm(PESON1 ~ . , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ ., data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2784.6  -270.8   -12.2   258.8  3063.5 
#> 
#> Coefficients:
#>                         Estimate Std. Error  t value             Pr(>|t|)    
#> (Intercept)           -3088.6457    22.1993 -139.133 < 0.0000000000000002 ***
#> SEMANAS                 158.5187     0.5426  292.158 < 0.0000000000000002 ***
#> SEXO1bebito             134.4620     1.9022   70.686 < 0.0000000000000002 ***
#> CESAREA.fSin cesárea    -13.5279     2.2083   -6.126       0.000000000903 ***
#> EDADM                    -0.7358     0.2414   -3.048             0.002306 ** 
#> EDADP                     3.0648     0.2055   14.913 < 0.0000000000000002 ***
#> ESTUDIOM.fMedios         -9.7726     2.5635   -3.812             0.000138 ***
#> ESTUDIOM.fUniversidad   -10.3297     2.8377   -3.640             0.000273 ***
#> ESTUDIOP.fMedios         11.9535     2.4401    4.899       0.000000965227 ***
#> ESTUDIOP.fUniversidad    25.1406     2.9862    8.419 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 420.2 on 195579 degrees of freedom
#> Multiple R-squared:  0.3167, Adjusted R-squared:  0.3166 
#> F-statistic: 1.007e+04 on 9 and 195579 DF,  p-value: < 0.00000000000000022

Generalmente las variables categóricas se introducen en los modelos mediante variables dummies. Crear dummies en R es sencillo, solo tienes que tener los datos como factor o como texto y R creará las dummies por ti cuando introduzcas la variable en lm(). Eso sí, es más fácil elegir la categoría de referencia si la variable es un factor.

Para entender como fija los regresores para las variables categóricas, mira esto:

levels(df$SEXO1)
#> [1] "bebita" "bebito"
levels(df$ESTUDIOM.f)
#> [1] "Primarios"   "Medios"      "Universidad"

Si necesitas cambiar la categoría de referencia siempre puedes usar forcast::fct_relevel()

zz <- forcats::fct_relevel(df$SEXO1, "bebito")
levels(zz)
#> [1] "bebito" "bebita"

Veamos que hay en el objeto mod_1

str(mod_1)
#listviewer::jsonedit(mod_1, mode = "view") ## Interactive option

Especificación del modelo

Lógicamente, a veces querremos seleccionar las variables explicativas:

mod_1 <- lm(PESON1 ~ SEMANAS, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1 + CESAREA.f + EDADM , data = df_m)
mod_1 <- lm(log(PESON1) ~ log(SEMANAS), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = log(PESON1) ~ log(SEMANAS), data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.73125 -0.08115  0.00386  0.08791  1.52234 
#> 
#> Coefficients:
#>               Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)  -0.245233   0.024450  -10.03 <0.0000000000000002 ***
#> log(SEMANAS)  2.269906   0.006671  340.27 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.142 on 195587 degrees of freedom
#> Multiple R-squared:  0.3718, Adjusted R-squared:  0.3718 
#> F-statistic: 1.158e+05 on 1 and 195587 DF,  p-value: < 0.00000000000000022

Si queremos introducir interacciones entre los regresores, podemos usar el operador :, aunque casi mejor hacerlo directamente con I():

mod_1 <- lm(PESON1 ~ SEMANAS + SEMANAS:EDADM, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*EDADM), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * EDADM), data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2780.80  -275.75   -13.51   262.58  3118.43 
#> 
#> Coefficients:
#>                       Estimate  Std. Error t value            Pr(>|t|)    
#> (Intercept)        -2898.86012    21.38108  -135.6 <0.0000000000000002 ***
#> SEMANAS              155.79835     0.56164   277.4 <0.0000000000000002 ***
#> I(SEMANAS * EDADM)     0.04972     0.00448    11.1 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 425.9 on 195586 degrees of freedom
#> Multiple R-squared:  0.298,  Adjusted R-squared:  0.2979 
#> F-statistic: 4.15e+04 on 2 and 195586 DF,  p-value: < 0.00000000000000022

Si queremos introducir las variables originales y también las interacciones entre ellas, podemos hacerlo directamente o o utilizar el operador *:

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM + SEMANAS:EDADM, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM + I(SEMANAS*EDADM), data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*EDADM, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*ESTUDIOM.f, data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS * ESTUDIOM.f, data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2792.55  -276.51   -14.41   263.49  3113.84 
#> 
#> Coefficients:
#>                                 Estimate Std. Error t value            Pr(>|t|)
#> (Intercept)                   -2862.7148    33.5901 -85.225 <0.0000000000000002
#> SEMANAS                         156.3815     0.8593 181.996 <0.0000000000000002
#> ESTUDIOM.fMedios               -104.1508    52.0312  -2.002              0.0453
#> ESTUDIOM.fUniversidad             6.3027    51.2415   0.123              0.9021
#> SEMANAS:ESTUDIOM.fMedios          2.6794     1.3301   2.014              0.0440
#> SEMANAS:ESTUDIOM.fUniversidad     0.1272     1.3093   0.097              0.9226
#>                                  
#> (Intercept)                   ***
#> SEMANAS                       ***
#> ESTUDIOM.fMedios              *  
#> ESTUDIOM.fUniversidad            
#> SEMANAS:ESTUDIOM.fMedios      *  
#> SEMANAS:ESTUDIOM.fUniversidad    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 426 on 195583 degrees of freedom
#> Multiple R-squared:  0.2976, Adjusted R-squared:  0.2976 
#> F-statistic: 1.658e+04 on 5 and 195583 DF,  p-value: < 0.00000000000000022

Recuerda que si queremos introducir algún regresor que sea la multiplicación de dos variables, tendremos que hacerlo con I()

mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*SEMANAS), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * SEMANAS), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2794.4  -274.4   -14.4   261.9  3609.5 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)          -5978.7888   151.5234  -39.46 <0.0000000000000002 ***
#> SEMANAS                324.5984     8.1524   39.82 <0.0000000000000002 ***
#> I(SEMANAS * SEMANAS)    -2.2567     0.1097  -20.57 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 425.6 on 195586 degrees of freedom
#> Multiple R-squared:  0.299,  Adjusted R-squared:  0.299 
#> F-statistic: 4.172e+04 on 2 and 195586 DF,  p-value: < 0.00000000000000022

También puede sernos de utilidad la función poly()

mod_1 <- lm(PESON1 ~ poly(SEMANAS, degree = 3), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ poly(SEMANAS, degree = 3), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2804.1  -270.2   -10.2   260.8  3316.1 
#> 
#> Coefficients:
#>                               Estimate  Std. Error t value            Pr(>|t|)
#> (Intercept)                  3254.6209      0.9586 3395.14 <0.0000000000000002
#> poly(SEMANAS, degree = 3)1 122613.9643    423.9502  289.22 <0.0000000000000002
#> poly(SEMANAS, degree = 3)2  -8755.4731    423.9502  -20.65 <0.0000000000000002
#> poly(SEMANAS, degree = 3)3 -16393.2469    423.9502  -38.67 <0.0000000000000002
#>                               
#> (Intercept)                ***
#> poly(SEMANAS, degree = 3)1 ***
#> poly(SEMANAS, degree = 3)2 ***
#> poly(SEMANAS, degree = 3)3 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 424 on 195585 degrees of freedom
#> Multiple R-squared:  0.3043, Adjusted R-squared:  0.3043 
#> F-statistic: 2.852e+04 on 3 and 195585 DF,  p-value: < 0.00000000000000022

Resultados de estimación

Una vez que sabemos la sintaxis de stats::lm(), vamos a ver como podemos acceder a la información que devuelve lm(). Ya hemos visto que los resultados se almacenan en una lista, así que podemos utilizar los operadores habituales de las listas para acceder a la información. Por ejemplo:

zz_betas     <- mod_1[[1]]

zz_residuals <- mod_1[[2]]
zz_residuals <- mod_1[["residuals"]]
zz_residuals <- mod_1$residuals

zz_betas_1   <- mod_1[[1]]      #- [[ ]] doble corchete
zz_betas_2 <- mod_1[1]          #- []    corchete simple

zz_betas_1a <- mod_1[["coefficients"]]   #- doble corchete 
zz_betas_1c <- mod_1$coefficients        #- $

zz_betas_2a <- mod_1["coefficients"]     #- single corchete

Afortunadamente ya hay construidas algunas funciones para manipular/ver los resultados de estimación. Por ejemplo:

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM + SEXO1 , data = df_m)

summary(mod_1)                   #- tabla resumen
summary(mod_1)$coefficients      #- tabla resumen con los coeficientes
coefficients(mod_1)              #- coeficientes estimados
confint(mod_1, level = 0.95)     #- Intervalos de confianza para los coeficientes
#fitted(mod_1)                   #- predicciones (y-sombrero, y-hat)
#residuals(mod_1)                #- vector de residuos
#model.matrix(mod_1)             #- extract the model matrix
anova(mod_1)                     #- ANOVA
vcov(mod_1)                      #- matriz de var-cov para los coeficientes 
#influence(mod_1)                #- regression diagnostics 
# diagnostic plots
layout(matrix(c(1, 2, 3, 4), 2, 2))   #- optional 4 graphs/page
plot(mod_1)                      #-
library(ggfortify)
autoplot(mod_1, which = 1:6, ncol = 2, colour = "steelblue")

Más utilidades

Hay muchas más funciones para valorar la idoneidad de un modelo. Por ejemplo aquí.

El paquete GGally permite hacer muchos análisis, por ejemplo:

• gráfico de los coeficientes estimados

GGally::ggcoef(mod_1)

• y muchas más cosas

df_jugar <- df_m %>% select(EDADM, SEXO1, SEMANAS)
GGally::ggpairs(df_jugar)

Predicciones

Para hacer predicciones para observaciones fuera de la muestra has de usar la función predict(). Has de proporcionar las observaciones a predecir como un df.

nuevas_observaciones <- df_m %>% slice(c(3, 44, 444))

predict(mod_1, newdata = nuevas_observaciones)  #- predicciones puntuales
predict(mod_1, newdata = nuevas_observaciones, type = 'response', se.fit = TRUE)  #- tb errores estándar  predictions
predict(mod_1, newdata = nuevas_observaciones, interval = "confidence")  #- intervalo (para el valor esperado)
predict(mod_1, newdata = nuevas_observaciones, interval = "prediction")  #- intervalo (para valores individuales)

Errores estándar robustos

Lo habitual es que al estimar un modelo, tengamos situaciones de heterocedasticidad, clustering etc … Podemos ajustar los errores estándar a estas situaciones.

El paquete de referencia para estos temas solía ser sandwich, aunque quizás ya haya ocupado su lugar el paquete estimatr.

Por ejemplo se pueden obtener errores estándar robustos con estimatr::lm_robust()³.

#- install.packages("emmeans")
mod_1 <- lm(PESON1 ~ SEMANAS, data = df_m)
mod_1_ee <- estimatr::lm_robust(PESON1 ~ SEMANAS, data = df_m)

Paquete broom

Un opción interesante es utilizar el paquete broom. Este paquete tiene 3 funciones útiles:

• tidy()

• augment()

• glance()

zz <- broom::tidy(mod_1, conf.int = TRUE)
zz %>% pjpv2020.01::pjp_f_decimales(nn = 2)  %>% gt::gt()

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	3254.62	0.96	3395.14	0	3252.74	3256.50
poly(SEMANAS, degree = 3)1	122613.96	423.95	289.22	0	121783.03	123444.90
poly(SEMANAS, degree = 3)2	-8755.47	423.95	-20.65	0	-9586.41	-7924.54
poly(SEMANAS, degree = 3)3	-16393.25	423.95	-38.67	0	-17224.18	-15562.31

mod_1 %>% broom::glance() %>% select(adj.r.squared, p.value)

adj.r.squared	p.value
0.3043378	0

broom::glance(mod_1)

r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
0.3043485	0.3043378	423.9502	28522.9	0	3	-1460765	2921540	2921591	35153224020	195585	195589

zz <- broom::augment(mod_1)

• un ejemplo sencillo en el que se ve la utilidad de broom

mod_1 %>% broom::tidy() %>% filter(p.value < 0.05)

term	estimate	std.error	statistic	p.value
(Intercept)	3254.621	0.9586114	3395.14100	0
poly(SEMANAS, degree = 3)1	122613.964	423.9501651	289.21787	0
poly(SEMANAS, degree = 3)2	-8755.473	423.9501651	-20.65213	0
poly(SEMANAS, degree = 3)3	-16393.247	423.9501651	-38.66786	0

• Un ejemplo de uso de broom, pero antes vamos a recordar alguna cosa de ggplot2:

ggplot(data = df_m, mapping = aes(x = EDADM, y = PESON1,  color = ESTUDIOM.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM, y = PESON1,  color = SEXO1)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM, y = PESON1,  color = SEXO1)) +
      geom_point(alpha = 0.1) +  geom_smooth()

Ahora sí viene el ejemplo en que se usa el paquete broom

mod_1 <- lm(PESON1 ~ EDADM + SEXO1 , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ EDADM + SEXO1, data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2908.94  -278.67    16.88   314.39  3038.84 
#> 
#> Coefficients:
#>              Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept) 3173.6776     6.9225 458.456 < 0.0000000000000002 ***
#> EDADM          0.5554     0.2072   2.681              0.00735 ** 
#> SEXO1bebito  121.9336     2.2832  53.405 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 504.6 on 195586 degrees of freedom
#> Multiple R-squared:  0.01441,    Adjusted R-squared:  0.0144 
#> F-statistic:  1430 on 2 and 195586 DF,  p-value: < 0.00000000000000022

td_mod_1 <- mod_1 %>% broom::augment(data = df_m) 

td_mod_1 %>% ggplot(mapping = aes(x = EDADM, y = PESON1, color = SEXO1)) +
                geom_point(alpha = 0.1) +
                geom_line(aes(y = .fitted, group = SEXO1))

td_mod_1 %>% ggplot(mapping = aes(x = EDADM, y = .fitted,  color = SEXO1)) +
                geom_line(aes(group = SEXO1)) 

(!!!!) Por ejemplo también permite fácilmente estimar modelos por grupos:

library(broom)
df_m %>% group_by(SEXO1) %>% do(tidy(lm(PESON1 ~ SEMANAS, .)))

SEXO1	term	estimate	std.error	statistic	p.value
bebita	(Intercept)	-2744.6382	30.4932289	-90.00812	0
bebita	SEMANAS	151.6970	0.7784646	194.86689	0
bebito	(Intercept)	-3078.8921	29.2617194	-105.21911	0
bebito	SEMANAS	163.6906	0.7484992	218.69176	0

O hacer gráficos de los intervalos de los coeficientes:

td_mod_1 <- tidy(mod_1, conf.int = TRUE)
ggplot(td_mod_1, aes(estimate, term, color = term)) +
    geom_point() +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
    geom_vline(xintercept = 0)

o este:

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line()

o este:

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line() +
  facet_wrap(vars(ESTUDIOM.f))

Comparación de modelos

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM,  data = df_m)
mod_2 <- lm(PESON1 ~ SEMANAS + EDADM + I(SEMANAS*EDADM), data = df_m)

anova(mod_1, mod_2)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
195586	35476827424	NA	NA	NA	NA
195585	35475031137	1	1796287	9.903494	0.0016499

AIC(mod_1, mod_2)

	df	AIC
mod_1	4	2923330
mod_2	5	2923323

lmtest::lrtest(mod_1, mod_2)    

#Df	LogLik	Df	Chisq	Pr(>Chisq)
4	-1461661	NA	NA	NA
5	-1461656	1	9.903446	0.0016497

(!!!) Vamos a ordenar los modelos en función de su AIC. Para ello vamos a crear una lista con modelos

modelos <- list(mod_1 <- lm(PESON1 ~ SEMANAS + EDADM,  data = df_m),
                mod_2 <- lm(PESON1 ~ SEMANAS + EDADM + I(SEMANAS*EDADM), data = df_m)  )

modelos_ordered_AIC <- purrr::map_df(modelos, broom::glance, .id = "model") %>% arrange(AIC)

modelos_ordered_AIC %>% gt::gt()

model	r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
2	0.2979802	0.2979694	425.8863	27672.75	0	3	-1461656	2923323	2923373	35475031137	195585	195589
1	0.2979446	0.2979374	425.8959	41502.28	0	2	-1461661	2923330	2923371	35476827424	195586	195589

• una tabla con broom (antes creamos una función (!!!!)):

my_kable <- function(df){ gt::gt(mutate_if(df, is.numeric, round, 2)) }

tidy(mod_1) %>% my_kable

term	estimate	std.error	statistic	p.value
(Intercept)	-2961.01	22.27	-132.97	0
SEMANAS	157.41	0.55	288.09	0
EDADM	1.92	0.17	10.96	0

El paquete modelr

Con el paquete modelr podemos fácilmente comparar las predicciones de varios modelos:

zz <- df_m %>% modelr::gather_predictions(mod_1, mod_2)

Para ver mejor las predicciones de los 2 modelos habrá que pasar zz a formato ancho:

zz1 <- pivot_wider(zz, names_from = model, values_from = pred)

---

6.2 Modelos GLM

Por ejemplo un Logit:

mod_logit <- glm(CESAREA.f ~ PESON1 + SEMANAS + EDADM + ESTUDIOM.f, family = binomial(link = "logit"), data = df_m)

summary(mod_logit)
#> 
#> Call:
#> glm(formula = CESAREA.f ~ PESON1 + SEMANAS + EDADM + ESTUDIOM.f, 
#>     family = binomial(link = "logit"), data = df_m)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.2698  -1.1990   0.7021   0.7882   2.0867  
#> 
#> Coefficients:
#>                          Estimate  Std. Error z value             Pr(>|z|)    
#> (Intercept)           -2.73466194  0.12159914 -22.489 < 0.0000000000000002 ***
#> PESON1                -0.00009625  0.00001227  -7.845  0.00000000000000432 ***
#> SEMANAS                0.14435995  0.00345507  41.782 < 0.0000000000000002 ***
#> EDADM                 -0.04652829  0.00104477 -44.534 < 0.0000000000000002 ***
#> ESTUDIOM.fMedios      -0.03820180  0.01340958  -2.849              0.00439 ** 
#> ESTUDIOM.fUniversidad  0.08638628  0.01325995   6.515  0.00000000007277406 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 221737  on 195588  degrees of freedom
#> Residual deviance: 217286  on 195583  degrees of freedom
#> AIC: 217298
#> 
#> Number of Fisher Scoring iterations: 4

En los modelos lineales, calcular efectos marginales es sencillo, pero el Logit es un modelo no lineal. Para calcular efectos marginales en modelos no lineales podemos usar el paquete margins. Por ejemplo, calculemos el Efecto marginal medio o average marginal effect (AME) en mod_logit.

mod_1_AME <- mod_logit %>% margins::margins() %>% summary() 
mod_1_AME

factor	AME	SE	z	p	lower	upper
EDADM	-0.0086132	0.0001907	-45.166365	0.0000000	-0.0089869	-0.0082394
ESTUDIOM.fMedios	-0.0072001	0.0025291	-2.846869	0.0044152	-0.0121571	-0.0022431
ESTUDIOM.fUniversidad	0.0158127	0.0024263	6.517178	0.0000000	0.0110572	0.0205682
PESON1	-0.0000178	0.0000023	-7.849686	0.0000000	-0.0000223	-0.0000134
SEMANAS	0.0267234	0.0006307	42.371560	0.0000000	0.0254873	0.0279595

Si queremos calcular efectos marginales para unos valores concretos de los regresores

mod_logit %>% margins::margins(at = list(SEMANAS = c(25, 35), ESTUDIOM.f = c("Primarios", "Medios", "Universidad")),  variables = "EDADM" ) 

Si prefieres visualizarlo, utiliza margins::cplot():

margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM", what = "effect", drop = TRUE)
margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM")

El pkg ggeffects

El paquete ggeffects proporciona otra forma de calcular efectos marginales en modelos de regresión:

Results of regression models are typically presented as tables that are easy to understand. For more complex models that include interaction or quadratic / spline terms, tables with numbers are less helpful and difficult to interpret. In such cases, marginal effects are far easier to understand. In particular, the visualization of marginal effects allows to intuitively get the idea of how predictors and outcome are associated, even for complex models.

No lo he probado aún, pero se puede usar en una gran variedad de modelos. Estimated Marginal Means and Marginal Effects from Regression Models

6.3 Tablas para modelos

Con sjPLot y friends …

m1 <- lm(PESON1 ~ SEMANAS + SEXO1 + EDADM + EDADP.2, data = df)
m2 <- lm(PESON1 ~ SEMANAS + SEXO1 + EDADM + EDADP.2 + ESTUDIOM.f + ESTUDIOP.f, data = df)
sjPlot::tab_model(m1)
sjPlot::plot_model(m1, sort.est = TRUE)
sjPlot::tab_model(m1, m2)

Con stargazer

mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1 , data = df_m)
stargazer::stargazer(mod_1, type = "html")

	Dependent variable:
	PESON1
SEMANAS	158.119***
	(0.540)
SEXO1bebito	134.645***
	(1.904)
Constant	-2,995.949***
	(21.163)
Observations	195,589
R2	0.315
Adjusted R2	0.315
Residual Std. Error	420.682 (df = 195586)
F Statistic	44,976.540*** (df = 2; 195586)
Note:	p<0.1; p<0.05; p<0.01

Con modelsummary

#remotes::install_github('vincentarelbundock/modelsummary')
library(modelsummary)

mys_modelitos <- list()
mys_modelitos[["PESON1:  OLS 1"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1,            df_m)
mys_modelitos[["PESON1:  OLS 2"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1 + EDADM , df_m)
mys_modelitos[["CESAREA: Logit 1"]]  <- glm(CESAREA.f ~ SEMANAS + SEXO1 , data = df_m, family = binomial(link = "logit"))

mm <- msummary(mys_modelitos, title = "Resultados de estimación")
mm

Resultados de estimación

	PESON1: OLS 1	PESON1: OLS 2	CESAREA: Logit 1
(Intercept)	-2995.949	-3063.898	-3.971
	(21.163)	(22.037)	(0.110)
SEMANAS	158.119	158.277	0.131
	(0.540)	(0.540)	(0.003)
SEXO1bebito	134.645	134.620	-0.105
	(1.904)	(1.903)	(0.010)
EDADM		1.903
		(0.173)
Num.Obs.	195589	195589	195589
R2	0.315	0.315
R2 Adj.	0.315	0.315
AIC	2918512.0	2918392.7	219454.5
BIC	2918552.7	2918443.6	219485.1
Log.Lik.	-1459251.990	-1459191.341	-109724.269
F	44976.535	30043.250

Con reports

• informes con el paquete reports

library(report) #- devtools::install_github("neuropsychology/report")
my_model <- lm(PESON1 ~ SEMANAS + SEXO1, df_m)
rr <- report(my_model, target = 1)
rr

report::as.report(rr)

Recuerda también que se puede obtener la ecuación (en latex) de un modelo con:

library(equatiomatic)  #- remotes::install_github("datalorax/equatiomatic")
extract_eq(mod_1)
extract_eq(mod_1, use_coefs = TRUE)

6.4 Otros modelos/técnicas

En realidad sólo voy a insistir en que con R se pueden implementar una gran variedad de modelos y técnicas estadísticas. Puedes ver algunos ejemplos en este libro. Hay materiales excelentes, tanto libros, como tutoriales, posts, etc … que puedes encontrar fácilmente en internet.

My #rstats learning path:

1. Install R
2. Install RStudio
3. Google “How do I [THING I WANT TO DO] in R?”

Repeat step 3 ad infinitum.

— Jesse Mostipak (@kierisi) August 18, 2017

---

6.5 Machine Learning

Simplemente señalar que es un área de estudio enorme y en constante evolución. En R hay varios entornos/enfoques para hacer ML. Creo que el que más futuro tiene es tidymodels. Un post introductorio sobre tidymodels. Un ejemplo de uso de tidymodels. Otro ejemplo.

Simplemente algunas referencias:

• Machine Learning for Everyone. Un post introductorio.

• Hands-On Machine Learning with R. Un buen bookdown

• 101 Machine Learning Algorithms for Data Science with Cheat Sheets. Cheatsheet de algoritmos ML.

• 101 Data Science Interview Questions, Answers, and Key Concepts. Post con cosas que “deberías” saber si quieres un trabajo como Data Scientist.

• Machine learning essentials. Un tema de un curso sobre ML.

• Introduction to Machine Learning with R. Otro curso sobre ML.

• Machine Learning. Bookdown centrado en explicar las principales ideas/conceptos de ML.

• Machine Learning con R y caretpor Joaquín Amat Rodrigo. A pesar de que el paquete caret ha sido sustituido por tidymodels, continua siendo una buena entrada a ML con R y en castellano, por Joaquín Amat Rodrigo.

• A Gentle Introduction to tidymodels.

• Choosing the right estimator. Guía de scitit learn para selección de modelo.

• Una introducción visual al machine learning - I. Post introductorio pero muy bonito visualmente sobre ML. Aquí la segunda parte.

• R: MLR, Decision Trees and Random Forest to Predict MPG for 2019 Vehicles. Después lo implementan en TensorFlow

---

1. Seguro que se me olvidan muchos; además el ecosistema R está en constante evolución, así que seguro que surgen mejoras↩︎

2. La verdad es que con la cantidad de materiales fantásticos que hay sobre R, no hay necesidad de explicar o escribir sobre todos los temas↩︎

3. Por defecto, el paquete usa Eicker-Huber-White robust standard errors, habitualmente conocidos como errores estándar “HC2”. Se pueden especificar otros métodos con la opción se_type. Por ejemplo se puede utilizar el método usado por defecto en Stata. Aquí puedes encontrar porque los errores estándar de Stata difieren de los usados en R y Phyton↩︎

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Web creada por Pedro J. Pérez. © 2020. If you find any bug please report them to pedro.j.perez@uv.es.