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 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 métodos, 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, David graba un vídeo en el que efectúa 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, en ficheros .Rmd, del código que acaba utilizando en los screencasts. Hay bastante gente que sigue sus vídeos y , como explican en este post, varias personas mantienen un listado en el que diseccionan las tareas que hace David en cada screencast, puedes encontrarlos aquí.

Por último, un buen 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 funciones y técnicas útiles en la fase de exploración de datos (EDA); aún así, durante el tutorial, 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, y 2) intentaremos contestar a la pregunta 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 2019 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.

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 tanto, solamente utilizaremos datos referentes a 2016 y 2019 para bebes nacidos en la Comunidad Valenciana.

Cargando el fichero de datos y dicc

Ya dijimos que los datos provienen del INE, concretamente de la Estadística de nacimientos. He restringido la muestra a los partos en la Comunitat Valenciana para los años 2016 y 2019.

#- Cargo datos q voy a usar y el diccionario. 
mys_archivos <- fs::dir_ls(here::here("datos", "usar"))         #- AQUI
dicc_df_all <- readxl::read_excel(mys_archivos[[2]])
df_all <-  readr::read_rds(mys_archivos[[1]])

Los datos tienen: 107 variables y 79830 registros de bebes. Tenemos datos para los años: 2016 y 2019. Hay muchas variables, iremos seleccionado las variables que nos interesan para cada visulización

Hay muchas variables, concretamente variables; así que en algún momento tendremos que centrarnos en algunas de ellas.

Diccionario

Lo primero es conocer tus variables, para ello tiene que haber un diccionario y lo hay, pero en nuestro caso el nombre de las variables (casí) es suficiente.

names(df_all)
#>   [1] "prov_inscrip"                       "prov_inscrip.n"                    
#>   [3] "muni_inscrip"                       "muni_inscrip.n"                    
#>   [5] "year_parto"                         "fecha_parto"                       
#>   [7] "prov_parto"                         "prov_parto.n"                      
#>   [9] "muni_parto"                         "muni_parto.n"                      
#>  [11] "lugar_parto.f"                      "parto_asistido.f"                  
#>  [13] "nn_nacidos_parto_con_o_sin_vida"    "parto_normal.f"                    
#>  [15] "parto_cesarea.f"                    "parto_37_semanas.f"                
#>  [17] "parto_nn_semanas"                   "fecha_nac_madre"                   
#>  [19] "nacionalidad_esp_madre"             "pais_nacionalidad_madre.n"         
#>  [21] "nacionalidad_madre_cuando.f"        "prov_nac_madre"                    
#>  [23] "prov_nac_madre.n"                   "muni_nac_madre"                    
#>  [25] "muni_nac_madre.n"                   "pais_naci_madre.n"                 
#>  [27] "prov_res_madre"                     "prov_res_madre.n"                  
#>  [29] "muni_res_madre"                     "muni_res_madre.n"                  
#>  [31] "pais_resi_madre.n"                  "estudios_madre.f"                  
#>  [33] "ocupacion_madre.f"                  "estado_civil_madre.f"              
#>  [35] "primer_matrimonio_madre.f"          "fecha_actual_matri_madre"          
#>  [37] "nn_anyos_casada"                    "intervalo_parto_matrimonio_madre.1"
#>  [39] "intervalo_parto_pareja_madre.2"     "intervalo_parto_matrimonio_madre.2"
#>  [41] "pareja_hecho_madre.f"               "primera_pareja_estable_madre.f"    
#>  [43] "fecha_actual_pareja_madre"          "nn_anyos_pareja_estable"           
#>  [45] "intervalo_parto_pareja_madre.1"     "nn_hijos_tot"                      
#>  [47] "nn_hijos_vivos_anteriores"          "fecha_parto_anterior"              
#>  [49] "intervalo_parto_anterior.2"         "intervalo_intergenesico"           
#>  [51] "intervalo_parto_anterior.1"         "fecha_nac_padre"                   
#>  [53] "nacionalidad_esp_padre"             "pais_nacionalidad_padre.n"         
#>  [55] "nacionalidad_padre_cuando.f"        "prov_nac_padre"                    
#>  [57] "muni_nac_padre"                     "pais_naci_padre.n"                 
#>  [59] "mismo_domicilio_padres.f"           "prov_res_padre"                    
#>  [61] "muni_res_padre"                     "estudios_padre.f"                  
#>  [63] "ocupacion_padre.f"                  "tamanyo_muni_inscri.f"             
#>  [65] "edad_madre"                         "edad_madre.1"                      
#>  [67] "edad_madre.2"                       "edad_madre_al_matrimonio"          
#>  [69] "edad_madre_inicio_rel_estable"      "edad_padre"                        
#>  [71] "edad_padre.1"                       "edad_padre.2"                      
#>  [73] "parto_multiple_generos.f"           "parto_multiple_generos.ff"         
#>  [75] "MULTSEX_duplicado"                  "parto_multiple_muertes.f"          
#>  [77] "parto_multiple_muertes.ff"          "MULTMFT_duplicado"                 
#>  [79] "sexo_bebe"                          "peso_bebe"                         
#>  [81] "bebe_muerto_cuando.f"               "bebe_muerto.f"                     
#>  [83] "nn_hijos_vivos_totales"             "muni_parto_tamanyo"                
#>  [85] "rel_actividad_madre.f"              "rel_actividad_padre.f"             
#>  [87] "MUNI.pob_t"                         "MUNI.pob_m"                        
#>  [89] "PROI.pob_t"                         "PROI.pob_m"                        
#>  [91] "ine_ccaa"                           "ine_ccaa.n"                        
#>  [93] "CCAAI.pob_t"                        "CCAAI.pob_m"                       
#>  [95] "partos_muni_T"                      "partos_prov_T"                     
#>  [97] "partos_ccaa_T"                      "estudios_madre.ff"                 
#>  [99] "dif_edad_padre_madre"               "dif_edad_padre_madre.1"            
#> [101] "dif_edad_padre_madre.2"             "estudios_padre.ff"                 
#> [103] "id_parto"                           "id_bebe_parto"                     
#> [105] "nn_bebes_muni_inscrip_T"            "nn_bebes_prov_inscrip_T"           
#> [107] "nn_bebes_ccaa_inscrip_T"

2 tablas

Hago 2 tablas para poder buscar las variables de forma rápida y ver sus principales características, aunque lo mejor es usar dicc_df_all ya sea en el Global o el diccionario original en local:

dicc_show_1 <- dicc_df_all %>% 
  select(variable, type, nn_unique, unique_values, p_na, p_zeros) %>% 
  mutate(unique_values = stringr::str_sub(unique_values, 1, 45))

1. Tabla con DT

DT::datatable(dicc_show_1, filter = 'top', extensions = "Scroller", 
              options = list(autoWidth = TRUE,deferRender = TRUE, 
                             scroller = TRUE, scrollY = 450 ))

1. Tabla con reactable

# Tabla con el paquete [`reactable`](https://glin.github.io/reactable/index.html)
dicc_show_2 <- dicc_df_all %>% 
  select(variable, type, nn_unique, unique_values, p_na, p_zeros) 
   
library(reactable)
reactable::reactable(dicc_show_2, pagination = FALSE, height = 450, filterable = TRUE, searchable = TRUE, highlight = TRUE, columns = list( variable = colDef(
     # sticky = "left",
      # Add a right border style to visually distinguish the sticky column
      style = list(borderRight = "1px solid #eee"),
      headerStyle = list(borderRight = "1px solid #eee")
    )),
  defaultColDef = colDef(minWidth = 70)
)

Datos simplificados

Para que todo vaya más fluido vamos a utilizar un fichero de datos reducido:

Al final, las variables que usaremos son:

#- Al final usaremos estas variables
df <- df_2
names(df)
#>  [1] "fecha_parto"               "prov_inscrip.n"           
#>  [3] "tamanyo_muni_inscri.f"     "parto_cesarea.f"          
#>  [5] "parto_normal.f"            "parto_nn_semanas"         
#>  [7] "sexo_bebe.f"               "peso_bebe"                
#>  [9] "pais_nacionalidad_madre.n" "edad_madre.1"             
#> [11] "estudios_madre.ff"         "pais_nacionalidad_padre.n"
#> [13] "edad_padre.1"              "estudios_padre.ff"        
#> [15] "nn_hijos_tot"

Vamos a tener 15 variables y 79830 bebes

Hagamos un diccionario para la base de datos reducida que al final usaremos:

#- creo dicc. de los daros q voy a USAR 
df_aa <- pjpv2020.01::pjp_f_estadisticos_basicos(df)
df_bb <- pjpv2020.01::pjp_f_unique_values(df, truncate = TRUE, nn_truncate = 50)
dicc_df <- bind_cols(df_aa, df_bb)
dicc_df <- dicc_df %>% select(variable, type, nn_unique, unique_values, q_na, p_na, p_zeros, min, max, mean, sd, NN, NN_ok)

DT::datatable(dicc_df, filter = 'top', extensions = "Scroller", options = list(
    autoWidth = TRUE, 
    deferRender = TRUE,
    scroller = TRUE,
    scrollY = 500
))

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.

Un post más reciente, de 2021, nos habla de 3 paquetes que hacen informes exploratorios completos de los datos: DataExplorer, SmartEDA y dlookr.

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

3.0 Informes completos

Hay varios paquetes que permiten realizar informes completos de un dataset de forma completamente automatizada. Por ejemplo:

Con el paquete DataExplorer

# informe sin especificar variable objetivo
DataExplorer::create_report(df)

# informe con  variable objetivo
create_report(df, y = "peso_bebe")
create_report(df, y = "sexo_bebe.f")

Con el paquete dlookr

dlookr::diagnose_report(df) 

# informe con variable objetivo
dlookr::eda_report(df, target = peso_bebe, 
                   output_format = "html", output_file = "EDA_bebes.html")

# ejemplo con NA's
dlookr::transformation_report(df, target = peso_bebe)

# ejemplo con  outliers
transformation_report(df, target = peso_bebe)

Con el paquete SmartEDA

SmartEDA::ExpReport(df, op_file = 'smarteda.html')

# informe con variable objetivo
SmartEDA::ExpReport(df, op_file = 'smarteda.html', Target = "peso_bebe")

Con el paquete summarytools

summarytools::dfSummary(df)

3.1 Nombres de las variables

Podemos cambiar y ver el nombre de las variables de varias maneras:

Con dplyr::rename()

df_copy <- df %>% dplyr::rename(peso_bebe_nuevo = peso_bebe) 
names(df_copy)[1] <- paste0(names(df_copy)[1], "_new") #- cambio el primer nombre

Con base::names()

names(df_copy)[1] <- paste0(names(df_copy)[1], "_new") #- cambio el primer nombre
names(df_copy)[2] <- "nuevo_nombre"    #- cambio el nombre de la segunda variable

la función janitor::clean_names() arregla automáticamente los nombres de las variables, lo que pasa es que es nuestro ejemplo los nombres ya están suficientemente arreglados, pero si tuviésemos un data.frame con nombres extraños, muy largos o incluso con nombres “no sintácticos”, y necesitaramos arreglarlos de forma rápida, podemos hacerlo con la función janitor::clean_names(), que arregla los nombres de las variables de forma automática, y además tiene algunas opciones para hacerlo como más te guste.

janitor::clean_names(df)

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 [79,830 × 15] (S3: tbl_df/tbl/data.frame)
#>  $ fecha_parto              : POSIXct[1:79830], format: "2016-01-01" "2016-01-01" ...
#>  $ prov_inscrip.n           : chr [1:79830] "Valencia/València" "Castellón/Castelló" "Valencia/València" "Valencia/València" ...
#>  $ tamanyo_muni_inscri.f    : Factor w/ 6 levels "[ 0 - 10.000]",..: 3 1 3 2 6 6 1 4 6 2 ...
#>  $ parto_cesarea.f          : Factor w/ 2 levels "Sin cesárea",..: 1 1 1 2 1 2 2 1 2 1 ...
#>  $ parto_normal.f           : Factor w/ 2 levels "Parto normal",..: 1 1 1 1 1 2 2 1 2 2 ...
#>  $ parto_nn_semanas         : num [1:79830] 23 26 25 25 27 27 25 26 30 27 ...
#>  $ sexo_bebe.f              : Factor w/ 2 levels "bebita","bebito": 1 1 1 1 2 1 2 2 1 2 ...
#>  $ peso_bebe                : num [1:79830] 560 565 685 750 750 800 850 960 1000 1030 ...
#>  $ pais_nacionalidad_madre.n: chr [1:79830] "España" "España" "España" "Bulgaria" ...
#>  $ edad_madre.1             : num [1:79830] 41.2 27 26.5 25.5 35.7 ...
#>  $ estudios_madre.ff        : Factor w/ 3 levels "Primarios","Medios",..: 2 1 2 1 2 3 2 NA 2 2 ...
#>  $ pais_nacionalidad_padre.n: chr [1:79830] "España" "España" "España" "Bulgaria" ...
#>  $ edad_padre.1             : num [1:79830] 47.4 34.8 30.1 27 37.9 ...
#>  $ estudios_padre.ff        : Factor w/ 3 levels "Primarios","Medios",..: 2 1 2 1 2 3 2 NA 2 1 ...
#>  $ nn_hijos_tot             : num [1:79830] 1 3 1 1 1 1 2 1 1 1 ...

b) con `inspectdf::inspect_types()`

inspectdf::inspect_types(df)   #- muestra de que tipo son las variables
#> # A tibble: 4 × 4
#>   type             cnt  pcnt col_name    
#>   <chr>          <int> <dbl> <named list>
#> 1 factor             6 40    <chr [6]>   
#> 2 numeric            5 33.3  <chr [5]>   
#> 3 character          3 20    <chr [3]>   
#> 4 POSIXct POSIXt     1  6.67 <chr [1]>

3.3 Valores únicos de las variables

Conviene saber qué valores (únicos) tienen las distintas variables de una tabla de datos; por ejemplo, es obvio, pero hay que saber cual es el periodo muestral de nuestros datos y el hecho de que nuestros datos los trabajemos habitualmente en formato largo, hace que saber esto no sea tan-tan sencillo

pjpv2020.01::pjp_f_unique_values() 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
fecha_parto	24	2016-01-01 - 2016-02-01 - 2016-03-01 - 2016-04-
prov_inscrip.n	3	Alicante/Alacant - Castellón/Castelló - Valenci
tamanyo_muni_inscri.f	6	[ 0 - 10.000] - [10.001 - 20.000] - [20.001 - 5
parto_cesarea.f	2	Sin cesárea - Con cesárea
parto_normal.f	2	Parto normal - Parto con complicaciones
parto_nn_semanas	28	NA - 20 - 21 - 22 - 23 - 24 - 25 - 26 - 27 - 28
sexo_bebe.f	2	bebita - bebito
peso_bebe	1974	NA - 390 - 400 - 430 - 450 - 500 - 510 - 520 -
pais_nacionalidad_madre.n	64	NA - Alemania - Angola - Argelia - Argentina -
edad_madre.1	1657	13.5824777549624 - 13.8288843258042 - 13.919233
estudios_madre.ff	4	NA - Primarios - Medios - Universidad
pais_nacionalidad_padre.n	64	NA - Alemania - Angola - Argelia - Argentina -
edad_padre.1	2034	NA - 14.8281998631075 - 15.0006844626968 - 15.2
estudios_padre.ff	4	NA - Primarios - Medios - Universidad
nn_hijos_tot	12	1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9 - 10 - 11 - 1

pjpv2020.01::pjp_f_unique_values(). Parecida a la función anterior, pero devuelve los valores únicos de cada variable en una columna: a veces esto es mejor para poder visualizarlos

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

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

names(df)
#>  [1] "fecha_parto"               "prov_inscrip.n"           
#>  [3] "tamanyo_muni_inscri.f"     "parto_cesarea.f"          
#>  [5] "parto_normal.f"            "parto_nn_semanas"         
#>  [7] "sexo_bebe.f"               "peso_bebe"                
#>  [9] "pais_nacionalidad_madre.n" "edad_madre.1"             
#> [11] "estudios_madre.ff"         "pais_nacionalidad_padre.n"
#> [13] "edad_padre.1"              "estudios_padre.ff"        
#> [15] "nn_hijos_tot"
zz <- df %>% distinct(prov_inscrip.n, parto_cesarea.f)
zz %>% gt::gt()

prov_inscrip.n	parto_cesarea.f
Valencia/València	Sin cesárea
Castellón/Castelló	Sin cesárea
Valencia/València	Con cesárea
Castellón/Castelló	Con cesárea
Alicante/Alacant	Con cesárea
Alicante/Alacant	Sin cesárea

3.4 Estadísticos básicos de las variables

a) con `pjpv2020.01::pjp_f_estadisticos_basicos()`

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)   #- estadísticos básicos del df
zz <- pjpv2020.01::pjp_f_unique_values(df)          #- valores únicos de cada variable de df

b) con `summarytools::dfSummary()`

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::skim()`

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

skimr::skim(df)

d) con `gtsummary::tbl_summary()`

gtsummary::tbl_summary(df)

Además esta función permite especificar variables de agrupación y realizar algunos contrastes. Para ver todas las posibilidades que ofrece, la vignette está aquí

df %>% select(peso_bebe, parto_nn_semanas, estudios_madre.f, parto_cesarea.f) %>% 
  gtsummary::tbl_summary(by = "parto_cesarea.f") %>% 
  gtsummary::add_p()
names(df)

3.5 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 …

Unas slides de @allison_horst que ayudan a entender los distintos tipos de NA’s. Además aquí está el video y un tutorial interactivo aquí

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

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 79.830 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)

dlookr::plot_na_intersect() proporciona un gráfico muy útil. Nos ofrece las variables que tienen NA’s y la co-ocurrencia de NA’s entre ellas.

dlookr::plot_na_intersect(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(estudios_madre.f)

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

¿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 15 variables, nos quedaríamos con 62.348 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 peso_bebe y/o en la variable parto_nn_semanas, lo haríamos así:

zz <- df %>% tidyr::drop_na(c(peso_bebe, parto_nn_semanas))   #- quito filas con NA en peso_bebe 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

Si por el contrario decides que lo adecuado es imputar valores a los NA’s, puedes utilizar el paquete mice, amelia II o dlookr::imputate_na(). Si decides tratar de predecir los valores de los NA’s, aquí tienes un post. Otra opción es usar

3.6 Variables numéricas

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

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

df_numeric <- df %>% select_if(is.numeric)      #- antigua sintaxis
df_numeric <- df %>% select(where(is.numeric))  #- nueva API

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: 79830  
#> 
#>                     edad_madre.1   edad_padre.1   nn_hijos_tot   parto_nn_semanas   peso_bebe
#> ----------------- -------------- -------------- -------------- ------------------ -----------
#>              Mean          33.02          36.09           1.72              38.85     3202.81
#>           Std.Dev           5.58           6.12           0.87               1.98      550.87
#>               Min          13.58          14.83           1.00              20.00      390.00
#>                Q1          29.67          32.33           1.00              38.00     2910.00
#>            Median          33.50          36.08           2.00              39.00     3240.00
#>                Q3          36.91          39.75           2.00              40.00     3550.00
#>               Max          53.92          73.49          12.00              46.00     6335.00
#>               MAD           5.31           5.56           1.48               1.48      474.43
#>               IQR           7.25           7.42           1.00               2.00      640.00
#>                CV           0.17           0.17           0.51               0.05        0.17
#>          Skewness          -0.40           0.19           1.80              -2.43       -0.74
#>       SE.Skewness           0.01           0.01           0.01               0.01        0.01
#>          Kurtosis           0.10           1.02           6.54              11.04        2.12
#>           N.Valid       79830.00       77204.00       79830.00           72995.00    76814.00
#>         Pct.Valid         100.00          96.71         100.00              91.44       96.22

Hay muchas otras alternativas, por ejemplo: dlookr::describe(), dlookr::univar_numeric(), dlookr::diagnose_numeric(), SmartEDA::ExpNumStat(), `summarytools::

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)
#>                  parto_nn_semanas peso_bebe edad_madre.1 edad_padre.1
#> parto_nn_semanas             1.00      0.62        -0.06        -0.03
#> peso_bebe                    0.62      1.00        -0.01         0.01
#> edad_madre.1                -0.06     -0.01         1.00         0.63
#> edad_padre.1                -0.03      0.01         0.63         1.00
#> nn_hijos_tot                -0.10      0.01         0.21         0.22
#>                  nn_hijos_tot
#> parto_nn_semanas        -0.10
#> peso_bebe                0.01
#> edad_madre.1             0.21
#> edad_padre.1             0.22
#> nn_hijos_tot             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	parto_nn_semanas	peso_bebe	edad_madre.1	edad_padre.1	nn_hijos_tot
parto_nn_semanas	NA	0.62	-0.06	-0.03	-0.10
peso_bebe	0.62	NA	-0.01	0.01	0.01
edad_madre.1	-0.06	-0.01	NA	0.63	0.21
edad_padre.1	-0.03	0.01	0.63	NA	0.22
nn_hijos_tot	-0.10	0.01	0.21	0.22	NA

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

df %>% inspectdf::inspect_cor()
#> # A tibble: 10 × 7
#>    col_1        col_2               corr   p_value    lower    upper pcnt_nna
#>    <chr>        <chr>              <dbl>     <dbl>    <dbl>    <dbl>    <dbl>
#>  1 edad_padre.1 edad_madre.1      0.635  0          0.631    0.639       96.7
#>  2 peso_bebe    parto_nn_semanas  0.617  0          0.612    0.622       89.5
#>  3 nn_hijos_tot edad_padre.1      0.220  0          0.214    0.227       96.7
#>  4 nn_hijos_tot edad_madre.1      0.206  0          0.200    0.213      100  
#>  5 nn_hijos_tot parto_nn_semanas -0.0994 7.71e-159 -0.107   -0.0922      91.4
#>  6 edad_madre.1 parto_nn_semanas -0.0558 2.49e- 51 -0.0630  -0.0486      91.4
#>  7 edad_padre.1 parto_nn_semanas -0.0317 3.55e- 17 -0.0390  -0.0243      88.6
#>  8 edad_madre.1 peso_bebe        -0.0140 1.06e-  4 -0.0211  -0.00691     96.2
#>  9 edad_padre.1 peso_bebe         0.0135 2.36e-  4  0.00630  0.0207      93.2
#> 10 nn_hijos_tot peso_bebe         0.0126 5.03e-  4  0.00548  0.0196      96.2

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")
#> # Correlation Matrix (pearson-method)
#> 
#> Parameter1       |   Parameter2 |     r |         95% CI |      t |    df |         p
#> -------------------------------------------------------------------------------------
#> parto_nn_semanas |    peso_bebe |  0.62 | [ 0.61,  0.62] | 209.62 | 71464 | < .001***
#> parto_nn_semanas | edad_madre.1 | -0.06 | [-0.06, -0.05] | -15.10 | 72993 | < .001***
#> parto_nn_semanas | edad_padre.1 | -0.03 | [-0.04, -0.02] |  -8.43 | 70752 | < .001***
#> parto_nn_semanas | nn_hijos_tot | -0.10 | [-0.11, -0.09] | -26.99 | 72993 | < .001***
#> peso_bebe        | edad_madre.1 | -0.01 | [-0.02, -0.01] |  -3.88 | 76812 | < .001***
#> peso_bebe        | edad_padre.1 |  0.01 | [ 0.01,  0.02] |   3.68 | 74416 | < .001***
#> peso_bebe        | nn_hijos_tot |  0.01 | [ 0.01,  0.02] |   3.48 | 76812 | < .001***
#> edad_madre.1     | edad_padre.1 |  0.63 | [ 0.63,  0.64] | 228.38 | 77202 | < .001***
#> edad_madre.1     | nn_hijos_tot |  0.21 | [ 0.20,  0.21] |  59.56 | 79828 | < .001***
#> edad_padre.1     | nn_hijos_tot |  0.22 | [ 0.21,  0.23] |  62.80 | 77202 | < .001***
#> 
#> p-value adjustment method: Holm (1979)
#> Observations: 70754-79830

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 estudios_madre.ff:

DataExplorer::plot_boxplot(df, by = "estudios_madre.ff")

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

my_vv <- names(df)[3]
my_vv <- "parto_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(edad_madre.1, estudios_madre.ff, target = parto_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(sexo_bebe.f, edad_madre.1, peso_bebe,  estudios_madre.ff, parto_cesarea.f) %>% explore::explore_all(target = parto_cesarea.f)

Veamos la relación entre “estudios_madre.ff” y CESAREA.f:

df %>% janitor::tabyl(estudios_madre.ff, parto_cesarea.f) %>% janitor::adorn_percentages()
#>  estudios_madre.ff Sin cesárea Con cesárea
#>          Primarios   0.7322450   0.2677550
#>             Medios   0.7162494   0.2837506
#>        Universidad   0.6801212   0.3198788
#>               <NA>   0.7010455   0.2989545

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: peso_bebe

my_vv <- "peso_bebe"
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.factor) %>% names()      #- old fashion
df %>% select(where(is.factor)) %>% names()

En realidad para seleccionar a las variables no-numéricas habría que hacer:

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)

DataExplorer::plot_bar(df, by = "sexo_bebe.f")

O con el paquete SmartEDA

SmartEDA::ExpCatViz(df, Page = c(3,3))

3.8 Outliers

Un outlier es una observación que difieres, que está alejada, del resto de observaciones. ¿Qué hacer con los outliers? pues como todo depende … Aquí y aquí tienes dos buenos posts sobre el tema.

Una forma rápida de identificar los outliers es usar performance::check_outliers() o dlookr::diagnose_outlier(). Para imputar un valor a los outliers se puede utilizar dlookr::imputate_outlier()

3.9 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(parto_cesarea.f)
df %>% explore::explore(edad_madre.1, estudios_madre.ff, target = parto_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$sexo_bebe.f, style = "rmarkdown")

Frequencies

df$sexo_bebe.f

Type: Factor

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
bebita	38582	48.33	48.33	48.33	48.33
bebito	41248	51.67	100.00	51.67	100.00
<NA>	0			0.00	100.00
Total	79830	100.00	100.00	100.00	100.00

tabulación cruzada entre dos variables categóricas:

summarytools::ctable(df$sexo_bebe.f, df$parto_cesarea.f)

Cross-Tabulation, Row Proportions
sexo_bebe.f * parto_cesarea.f
Data Frame: df

	parto_cesarea.f	Sin cesárea	Con cesárea	Total
sexo_bebe.f
bebita		27599 (71.5%)	10983 (28.5%)	38582 (100.0%)
bebito		28380 (68.8%)	12868 (31.2%)	41248 (100.0%)
Total		55979 (70.1%)	23851 (29.9%)	79830 (100.0%)

Incluso se pueden hacer test chi-cuadrado

summarytools::ctable(df$estudios_madre.ff, df$parto_normal.f, chisq = TRUE)

Cross-Tabulation, Row Proportions
estudios_madre.ff * parto_normal.f
Data Frame: df

	parto_normal.f	Parto normal	Parto con complicaciones	Total
estudios_madre.ff
Primarios		6167 (86.0%)	1000 (14.0%)	7167 (100.0%)
Medios		28058 (84.9%)	4971 (15.1%)	33029 (100.0%)
Universidad		29656 (87.2%)	4335 (12.8%)	33991 (100.0%)
		4830 (85.6%)	813 (14.4%)	5643 (100.0%)
Total		68711 (86.1%)	11119 (13.9%)	79830 (100.0%)

Chi.squared	df	p.value
73.9143	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(parto_cesarea.f) %>% gt::gt()

parto_cesarea.f	n	percent
Sin cesárea	55979	0.7012276
Con cesárea	23851	0.2987724

df %>% janitor::tabyl(parto_cesarea.f, estudios_madre.ff) %>% gt::gt()

parto_cesarea.f	Primarios	Medios	Universidad	NA_
Sin cesárea	5248	23657	23118	3956
Con cesárea	1919	9372	10873	1687

df %>% janitor::tabyl(parto_cesarea.f, estudios_madre.ff, sexo_bebe.f)  #- xq no podemos usar gt::gt()
#> $bebita
#>  parto_cesarea.f Primarios Medios Universidad  NA_
#>      Sin cesárea      2639  11639       11448 1873
#>      Con cesárea       884   4286        5077  736
#> 
#> $bebito
#>  parto_cesarea.f Primarios Medios Universidad  NA_
#>      Sin cesárea      2609  12018       11670 2083
#>      Con cesárea      1035   5086        5796  951

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$parto_cesarea.f, df$estudios_madre.ff, useNA = "always") 
#>              
#>               Primarios Medios Universidad  <NA>
#>   Sin cesárea      5248  23657       23118  3956
#>   Con cesárea      1919   9372       10873  1687
#>   <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$parto_cesarea.f, df$estudios_madre.ff) 
my_tabla %>% as.data.frame() %>% 
            pivot_wider(names_from = 2, values_from = 3) %>% 
            gt::gt()

Var1	Primarios	Medios	Universidad
Sin cesárea	5248	23657	23118
Con cesárea	1919	9372	10873

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$peso_bebe, mu = 3250)
#> 
#>  One Sample t-test
#> 
#> data:  df$peso_bebe
#> t = -23.744, df = 76813, p-value < 0.00000000000000022
#> alternative hypothesis: true mean is not equal to 3250
#> 95 percent confidence interval:
#>  3198.911 3206.702
#> sample estimates:
#> mean of x 
#>  3202.807
t.test(df$peso_bebe ~ df$parto_cesarea.f)
#> 
#>  Welch Two Sample t-test
#> 
#> data:  df$peso_bebe by df$parto_cesarea.f
#> t = 22.581, df = 34951, p-value < 0.00000000000000022
#> alternative hypothesis: true difference in means between group Sin cesárea and group Con cesárea is not equal to 0
#> 95 percent confidence interval:
#>   99.23847 118.10410
#> sample estimates:
#> mean in group Sin cesárea mean in group Con cesárea 
#>                  3235.470                  3126.799

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(estudios_madre.ff) %>% 
  summarise(t_test = list(t.test(peso_bebe))) %>% 
  mutate(tidied = map(t_test, tidy)) %>% 
  tidyr::unnest(tidied) %>% 
  ggplot(aes(estimate, estudios_madre.ff)) + 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$parto_cesarea.f, df$sexo_bebe.f)
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  df$parto_cesarea.f and df$sexo_bebe.f
#> X-squared = 70.787, df = 1, p-value < 0.00000000000000022
chisq.test(df$parto_cesarea.f, df$estudios_madre.ff) 
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  df$parto_cesarea.f and df$estudios_madre.ff
#> X-squared = 140.77, 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(peso_bebe, parto_nn_semanas, sexo_bebe.f, parto_cesarea.f, edad_madre.1, edad_padre.1, estudios_madre.ff, estudios_padre.ff)
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 (peso_bebe) en función de todas las demás variables en df_m.

mod_1 <- lm(peso_bebe ~ . , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ ., data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2684.2  -272.3   -13.0   260.1  4039.2 
#> 
#> Coefficients:
#>                                Estimate Std. Error t value             Pr(>|t|)
#> (Intercept)                  -3690.5358    37.5808 -98.203 < 0.0000000000000002
#> parto_nn_semanas               174.4397     0.8868 196.710 < 0.0000000000000002
#> sexo_bebe.fbebito              135.8581     3.3925  40.046 < 0.0000000000000002
#> parto_cesarea.fCon cesárea       8.3805     3.7848   2.214             0.026816
#> edad_madre.1                    -0.3357     0.4456  -0.753             0.451262
#> edad_padre.1                     2.7483     0.3799   7.234    0.000000000000475
#> estudios_madre.ffMedios        -32.8690     6.7554  -4.866    0.000001143931279
#> estudios_madre.ffUniversidad   -23.6951     7.1633  -3.308             0.000941
#> estudios_padre.ffMedios        -24.1311     6.2499  -3.861             0.000113
#> estudios_padre.ffUniversidad   -14.3470     6.9183  -2.074             0.038104
#>                                 
#> (Intercept)                  ***
#> parto_nn_semanas             ***
#> sexo_bebe.fbebito            ***
#> parto_cesarea.fCon cesárea   *  
#> edad_madre.1                    
#> edad_padre.1                 ***
#> estudios_madre.ffMedios      ***
#> estudios_madre.ffUniversidad ***
#> estudios_padre.ffMedios      ***
#> estudios_padre.ffUniversidad *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 424.1 on 62688 degrees of freedom
#> Multiple R-squared:  0.3956, Adjusted R-squared:  0.3955 
#> F-statistic:  4559 on 9 and 62688 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$sexo_bebe.f)
#> [1] "bebita" "bebito"
levels(df$estudios_madre.ff)
#> [1] "Primarios"   "Medios"      "Universidad"

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

zz <- forcats::fct_relevel(df$sexo_bebe.f, "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(peso_bebe ~ parto_nn_semanas, data = df_m)
mod_1 <- lm(peso_bebe ~ parto_nn_semanas + sexo_bebe.f + parto_cesarea.f + edad_madre.1 , data = df_m)
mod_1 <- lm(log(peso_bebe) ~ log(parto_nn_semanas), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = log(peso_bebe) ~ log(parto_nn_semanas), data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.67055 -0.08264  0.00416  0.09018  1.63927 
#> 
#> Coefficients:
#>                       Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)           -1.37399    0.03976  -34.55 <0.0000000000000002 ***
#> log(parto_nn_semanas)  2.57713    0.01087  237.14 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1465 on 62696 degrees of freedom
#> Multiple R-squared:  0.4728, Adjusted R-squared:  0.4728 
#> F-statistic: 5.624e+04 on 1 and 62696 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(peso_bebe ~ parto_nn_semanas + parto_nn_semanas:edad_madre.1, data = df_m)
mod_1 <- lm(peso_bebe ~ parto_nn_semanas + I(parto_nn_semanas*edad_madre.1), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ parto_nn_semanas + I(parto_nn_semanas * 
#>     edad_madre.1), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2663.1  -277.3   -14.2   263.2  4067.8 
#> 
#> Coefficients:
#>                                        Estimate   Std. Error  t value
#> (Intercept)                        -3525.030516    34.468503 -102.268
#> parto_nn_semanas                     171.637088     0.912705  188.053
#> I(parto_nn_semanas * edad_madre.1)     0.046324     0.008384    5.526
#>                                                Pr(>|t|)    
#> (Intercept)                        < 0.0000000000000002 ***
#> parto_nn_semanas                   < 0.0000000000000002 ***
#> I(parto_nn_semanas * edad_madre.1)          0.000000033 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 429.9 on 62695 degrees of freedom
#> Multiple R-squared:  0.379,  Adjusted R-squared:  0.3789 
#> F-statistic: 1.913e+04 on 2 and 62695 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(peso_bebe ~ parto_nn_semanas + edad_madre.1 + parto_nn_semanas:edad_madre.1, data = df_m)
mod_1 <- lm(peso_bebe ~ parto_nn_semanas + edad_madre.1 + I(parto_nn_semanas*edad_madre.1), data = df_m)
mod_1 <- lm(peso_bebe ~ parto_nn_semanas*edad_madre.1, data = df_m)
mod_1 <- lm(peso_bebe ~ parto_nn_semanas*estudios_madre.ff, data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ parto_nn_semanas * estudios_madre.ff, 
#>     data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2662.5  -278.4   -14.4   263.6  4077.9 
#> 
#> Coefficients:
#>                                                 Estimate Std. Error t value
#> (Intercept)                                   -3454.8834   114.2206 -30.247
#> parto_nn_semanas                                172.0883     2.9322  58.689
#> estudios_madre.ffMedios                         -74.5179   124.9941  -0.596
#> estudios_madre.ffUniversidad                    -52.2686   125.2082  -0.417
#> parto_nn_semanas:estudios_madre.ffMedios          0.9322     3.2094   0.290
#> parto_nn_semanas:estudios_madre.ffUniversidad     0.7607     3.2149   0.237
#>                                                          Pr(>|t|)    
#> (Intercept)                                   <0.0000000000000002 ***
#> parto_nn_semanas                              <0.0000000000000002 ***
#> estudios_madre.ffMedios                                     0.551    
#> estudios_madre.ffUniversidad                                0.676    
#> parto_nn_semanas:estudios_madre.ffMedios                    0.771    
#> parto_nn_semanas:estudios_madre.ffUniversidad               0.813    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 429.9 on 62692 degrees of freedom
#> Multiple R-squared:  0.3791, Adjusted R-squared:  0.379 
#> F-statistic:  7655 on 5 and 62692 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(peso_bebe ~ parto_nn_semanas + I(parto_nn_semanas*parto_nn_semanas), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ parto_nn_semanas + I(parto_nn_semanas * 
#>     parto_nn_semanas), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2634.3  -276.6   -14.3   262.8  4317.7 
#> 
#> Coefficients:
#>                                          Estimate Std. Error t value
#> (Intercept)                            -5768.0193   228.5009  -25.24
#> parto_nn_semanas                         297.0664    12.4745   23.81
#> I(parto_nn_semanas * parto_nn_semanas)    -1.6984     0.1702   -9.98
#>                                                   Pr(>|t|)    
#> (Intercept)                            <0.0000000000000002 ***
#> parto_nn_semanas                       <0.0000000000000002 ***
#> I(parto_nn_semanas * parto_nn_semanas) <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 429.7 on 62695 degrees of freedom
#> Multiple R-squared:  0.3796, Adjusted R-squared:  0.3796 
#> F-statistic: 1.918e+04 on 2 and 62695 DF,  p-value: < 0.00000000000000022

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

mod_1 <- lm(peso_bebe ~ poly(parto_nn_semanas, degree = 3), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ poly(parto_nn_semanas, degree = 3), 
#>     data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2663.2  -273.2   -12.8   256.8  3608.3 
#> 
#> Coefficients:
#>                                       Estimate Std. Error t value
#> (Intercept)                           3205.032      1.703 1881.65
#> poly(parto_nn_semanas, degree = 3)1  84052.196    426.500  197.07
#> poly(parto_nn_semanas, degree = 3)2  -4287.939    426.500  -10.05
#> poly(parto_nn_semanas, degree = 3)3 -13032.791    426.500  -30.56
#>                                                Pr(>|t|)    
#> (Intercept)                         <0.0000000000000002 ***
#> poly(parto_nn_semanas, degree = 3)1 <0.0000000000000002 ***
#> poly(parto_nn_semanas, degree = 3)2 <0.0000000000000002 ***
#> poly(parto_nn_semanas, degree = 3)3 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 426.5 on 62694 degrees of freedom
#> Multiple R-squared:  0.3888, Adjusted R-squared:  0.3887 
#> F-statistic: 1.329e+04 on 3 and 62694 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(peso_bebe ~ parto_nn_semanas + edad_madre.1 + sexo_bebe.f , 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(edad_madre.1, sexo_bebe.f, parto_nn_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(peso_bebe ~ parto_nn_semanas, data = df_m)
mod_1_ee <- estimatr::lm_robust(peso_bebe ~ parto_nn_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	conf.low	conf.high
(Intercept)	3205.03	1.7	1881.65	3201.69	3208.37
poly(parto_nn_semanas, degree = 3)1	84052.20	426.5	197.07	83216.25	84888.14
poly(parto_nn_semanas, degree = 3)2	-4287.94	426.5	-10.05	-5123.88	-3452.00
poly(parto_nn_semanas, degree = 3)3	-13032.79	426.5	-30.56	-13868.73	-12196.85

mod_1 %>% broom::glance() %>% select(adj.r.squared, p.value)
#> # A tibble: 1 × 2
#>   adj.r.squared p.value
#>           <dbl>   <dbl>
#> 1         0.389       0
broom::glance(mod_1)
#> # A tibble: 1 × 12
#>   r.squared adj.r.sq…¹ sigma stati…² p.value    df  logLik    AIC    BIC devia…³
#>       <dbl>      <dbl> <dbl>   <dbl>   <dbl> <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
#> 1     0.389      0.389  427.  13291.       0     3 -4.69e5 9.37e5 9.37e5 1.14e10
#> # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
#> #   variable names ¹adj.r.squared, ²statistic, ³deviance

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)
#> # A tibble: 4 × 5
#>   term                                estimate std.error statistic   p.value
#>   <chr>                                  <dbl>     <dbl>     <dbl>     <dbl>
#> 1 (Intercept)                            3205.      1.70    1882.  0        
#> 2 poly(parto_nn_semanas, degree = 3)1   84052.    427.       197.  0        
#> 3 poly(parto_nn_semanas, degree = 3)2   -4288.    427.       -10.1 9.22e- 24
#> 4 poly(parto_nn_semanas, degree = 3)3  -13033.    427.       -30.6 1.42e-203

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

ggplot(data = df_m, mapping = aes(x = edad_madre.1, y = peso_bebe,  color = estudios_madre.ff)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = edad_madre.1, y = peso_bebe,  color = sexo_bebe.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = edad_madre.1, y = peso_bebe,  color = sexo_bebe.f)) +
      geom_point(alpha = 0.1) +  geom_smooth()

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

mod_1 <- lm(peso_bebe ~ edad_madre.1 + sexo_bebe.f , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = peso_bebe ~ edad_madre.1 + sexo_bebe.f, data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2785.27  -279.03    37.17   340.48  2642.20 
#> 
#> Coefficients:
#>                    Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)       3219.9816    14.0263 229.567 < 0.0000000000000002 ***
#> edad_madre.1        -2.2483     0.4102  -5.481         0.0000000423 ***
#> sexo_bebe.fbebito  116.0561     4.3340  26.778 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 542.3 on 62695 degrees of freedom
#> Multiple R-squared:  0.01178,    Adjusted R-squared:  0.01175 
#> F-statistic: 373.7 on 2 and 62695 DF,  p-value: < 0.00000000000000022

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

td_mod_1 %>% ggplot(mapping = aes(x = edad_madre.1, y = peso_bebe, color = sexo_bebe.f)) +
                geom_point(alpha = 0.1) +
                geom_line(aes(y = .fitted, group = sexo_bebe.f))

td_mod_1 %>% ggplot(mapping = aes(x = edad_madre.1, y = .fitted,  color = sexo_bebe.f)) +
                geom_line(aes(group = sexo_bebe.f))

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

library(broom)
df_m %>% group_by(sexo_bebe.f) %>% do(tidy(lm(peso_bebe ~ parto_nn_semanas, .)))
#> # A tibble: 4 × 6
#> # Groups:   sexo_bebe.f [2]
#>   sexo_bebe.f term             estimate std.error statistic p.value
#>   <fct>       <chr>               <dbl>     <dbl>     <dbl>   <dbl>
#> 1 bebita      (Intercept)        -3385.     48.8      -69.4       0
#> 2 bebita      parto_nn_semanas     168.      1.25     134.        0
#> 3 bebito      (Intercept)        -3696.     47.4      -78.0       0
#> 4 bebito      parto_nn_semanas     179.      1.22     147.        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 = parto_nn_semanas , y = .fitted, color = sexo_bebe.f)) +
   geom_point(mapping = aes(y = peso_bebe), alpha = 0.1) +
   geom_line()

o este:

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = parto_nn_semanas , y = .fitted, color = sexo_bebe.f)) +
   geom_point(mapping = aes(y = peso_bebe), alpha = 0.1) +
   geom_line() +
  facet_wrap(vars(estudios_madre.ff))

Comparación de modelos

mod_1 <- lm(peso_bebe ~ parto_nn_semanas + edad_madre.1,  data = df_m)
mod_2 <- lm(peso_bebe ~ parto_nn_semanas + edad_madre.1 + I(parto_nn_semanas*edad_madre.1), data = df_m)

anova(mod_1, mod_2)
#> Analysis of Variance Table
#> 
#> Model 1: peso_bebe ~ parto_nn_semanas + edad_madre.1
#> Model 2: peso_bebe ~ parto_nn_semanas + edad_madre.1 + I(parto_nn_semanas * 
#>     edad_madre.1)
#>   Res.Df         RSS Df Sum of Sq      F         Pr(>F)    
#> 1  62695 11587395452                                       
#> 2  62694 11581005094  1   6390357 34.594 0.000000004081 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(mod_1, mod_2)
#>       df      AIC
#> mod_1  4 938282.1
#> mod_2  5 938249.5

lmtest::lrtest(mod_1, mod_2)    
#> Likelihood ratio test
#> 
#> Model 1: peso_bebe ~ parto_nn_semanas + edad_madre.1
#> Model 2: peso_bebe ~ parto_nn_semanas + edad_madre.1 + I(parto_nn_semanas * 
#>     edad_madre.1)
#>   #Df  LogLik Df  Chisq     Pr(>Chisq)    
#> 1   4 -469137                             
#> 2   5 -469120  1 34.587 0.000000004076 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(!!!) 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(peso_bebe ~ parto_nn_semanas + edad_madre.1,  data = df_m),
                mod_2 <- lm(peso_bebe ~ parto_nn_semanas + edad_madre.1 + I(parto_nn_semanas*edad_madre.1), 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.3792746	0.3792449	429.7938	12769.06	0	3	-469119.7	938249.5	938294.7	11581005094	62694	62698
1	0.3789321	0.3789123	429.9089	19126.05	0	2	-469137.0	938282.1	938318.2	11587395452	62695	62698

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
(Intercept)	-3581.53	36.78	-97.37
parto_nn_semanas	173.17	0.89	195.46
edad_madre.1	1.70	0.33	5.22

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(parto_cesarea.f ~ peso_bebe + parto_nn_semanas + edad_madre.1 + estudios_madre.ff, family = binomial(link = "logit"), data = df_m)

summary(mod_logit)
#> 
#> Call:
#> glm(formula = parto_cesarea.f ~ peso_bebe + parto_nn_semanas + 
#>     edad_madre.1 + estudios_madre.ff, family = binomial(link = "logit"), 
#>     data = df_m)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.3656  -0.8545  -0.7289   1.3027   2.2622  
#> 
#> Coefficients:
#>                                 Estimate  Std. Error z value
#> (Intercept)                   4.48309769  0.20719072  21.638
#> peso_bebe                     0.00007551  0.00002072   3.644
#> parto_nn_semanas             -0.19393520  0.00592605 -32.726
#> edad_madre.1                  0.05964899  0.00184411  32.346
#> estudios_madre.ffMedios      -0.05914193  0.03568428  -1.657
#> estudios_madre.ffUniversidad -0.06553869  0.03595874  -1.823
#>                                          Pr(>|z|)    
#> (Intercept)                  < 0.0000000000000002 ***
#> peso_bebe                                0.000269 ***
#> parto_nn_semanas             < 0.0000000000000002 ***
#> edad_madre.1                 < 0.0000000000000002 ***
#> estudios_madre.ffMedios                  0.097445 .  
#> estudios_madre.ffUniversidad             0.068363 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 76719  on 62697  degrees of freedom
#> Residual deviance: 73750  on 62692  degrees of freedom
#> AIC: 73762
#> 
#> 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
#>                  edad_madre.1  0.0119 0.0004  33.2356 0.0000  0.0112  0.0127
#>       estudios_madre.ffMedios -0.0120 0.0073  -1.6452 0.0999 -0.0262  0.0023
#>  estudios_madre.ffUniversidad -0.0132 0.0073  -1.8080 0.0706 -0.0276  0.0011
#>              parto_nn_semanas -0.0388 0.0012 -33.7301 0.0000 -0.0411 -0.0366
#>                     peso_bebe  0.0000 0.0000   3.6445 0.0003  0.0000  0.0000

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

mod_logit %>% margins::margins(at = list(parto_nn_semanas = c(25, 35), estudios_madre.ff = c("Primarios", "Medios", "Universidad")),  variables = "edad_madre.1" )

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

margins::cplot(mod_logit, x = "estudios_madre.ff", dx = "edad_madre.1", what = "effect", drop = TRUE)
margins::cplot(mod_logit, x = "estudios_madre.ff", dx = "edad_madre.1")

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(peso_bebe ~ parto_nn_semanas + sexo_bebe.f + edad_madre.1 + edad_padre.1, data = df)
m2 <- lm(peso_bebe ~ parto_nn_semanas + sexo_bebe.f + edad_madre.1 + edad_padre.1. + estudios_madre.ff + estudios_padre.ff, data = df)
sjPlot::tab_model(m1)
sjPlot::plot_model(m1, sort.est = TRUE)
sjPlot::tab_model(m1, m2)

Con `stargazer`

mod_1 <- lm(peso_bebe ~ parto_nn_semanas + sexo_bebe.f , data = df_m)
stargazer::stargazer(mod_1, type = "html")


	Dependent variable:

	peso_bebe

parto_nn_semanas	173.917^***
	(0.874)

sexo_bebe.fbebito	135.968^***
	(3.395)

Constant	-3,623.934^***
	(34.094)


Observations	62,698
R²	0.394
Adjusted R²	0.394
Residual Std. Error	424.605 (df = 62695)
F Statistic	20,394.880^*** (df = 2; 62695)

Note:	p<0.1; p<0.05; p<0.01

Con `modelsummary`

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

mys_modelitos <- list()
mys_modelitos[["peso_bebe:  OLS 1"]]    <-  lm(peso_bebe    ~ parto_nn_semanas + sexo_bebe.f,            df_m)
mys_modelitos[["peso_bebe:  OLS 2"]]    <-  lm(peso_bebe    ~ parto_nn_semanas + sexo_bebe.f + edad_madre.1 , df_m)
mys_modelitos[["CESAREA: Logit 1"]]  <- glm(parto_cesarea.f ~ parto_nn_semanas + sexo_bebe.f , data = df_m, family = binomial(link = "logit"))

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

Resultados de estimación
	peso_bebe: OLS 1	peso_bebe: OLS 2	CESAREA: Logit 1
(Intercept)	−3623.934	−3693.200	6.311
	(34.094)	(36.427)	(0.179)
parto_nn_semanas	173.917	174.210	−0.186
	(0.874)	(0.875)	(0.005)
sexo_bebe.fbebito	135.968	136.014	0.121
	(3.395)	(3.394)	(0.018)
edad_madre.1		1.735
		(0.322)
Num.Obs.	62698	62698	62698
R2	0.394	0.394
R2 Adj.	0.394	0.394
AIC	936725.4	936698.3	74897.6
BIC	936761.6	936743.6	74924.8
Log.Lik.	−468358.705	−468344.169	−37445.810
F	20394.883	13612.370	844.944
RMSE	424.59	424.50	0.45

Con `reports`

informes con el paquete reports

library(report) #- devtools::install_github("neuropsychology/report")
my_model <- lm(peso_bebe ~ parto_nn_semanas + sexo_bebe.f, 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 scikit-learn para selección de modelos.
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

Seguro que se me olvidan muchos; además el ecosistema R está en constante evolución, así que seguro que surgen mejoras↩︎
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↩︎
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↩︎

EDA y modelos (datos de bebes)

Pedro J. Pérez (pedro.j.perez@uv.es). Universitat de València Web del curso: https://perezp44.github.io/intro-ds-22-23-web/

Septiembre de 2020 (actualizado el 07 septiembre 2022)