EDA y modelos (datos de bebes)
Pedro J. Pérez (pedro.j.perez@uv.es). Universitat de València
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_allya 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))- Tabla con
DT
- Tabla con
DT::datatable(dicc_show_1, filter = 'top', extensions = "Scroller",
options = list(autoWidth = TRUE,deferRender = TRUE,
scroller = TRUE, scrollY = 450 ))- Tabla con
reactable
- Tabla con
# 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, estos1:
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ónjanitor::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 dfb) 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 informec) 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 filasPero 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-RSi 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íasSi 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 APIa) 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.22Hay 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-79830d) 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.2989545e) 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)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.
summarytoolsis 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 9514.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
gtprimero las hemos de convertir a data.frame y casi seguro que habrá que hacer uso depivot_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(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 poco2) 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.00000000000000022Generalmente 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 optionEspecificació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.00000000000000022Si 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.00000000000000022Si 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.00000000000000022Recuerda 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.00000000000000022Tambié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.00000000000000022Resultados 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 corcheteAfortunadamente 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)
Esto es un título de primer nivel
1. Esto es un título de segundo nivel
2. Formateo de texto
Esto es texo normal
Esto se mostrará en negrita
Esto se mostrará en cursiva
Esto se mostrará tachado
Esto se mostrará como una cita
3. Listas
Listas con números:
- Primero
- Segundo
- Tercero
Listas con bullets:
- Primero
- segundo
- Tercero
Mas listas:
- Primero
- 1.1
- 1.2
- Segundo
- a
- b
- Tercero
Sub-indices y ..
Un superíndice: Y2
Un subíndice: Y2
Imágenes
Para mostrar una imágen basta con: 
También podríamos haberlo hecho así (da un poco más flexibilidad):
Formulas
Las formulas hay que escribirlas en Latex. Yo no sé mucho Latex, pero hay editores que facilitan mucho el trabajo: preguntad al profe!!
Una fórmula dentro del texto se escribiría así: \(E = mc^{2}\) ; es decir, enmarcándola por $
Si quieres la formula ella sola hay que poner dos $$:
\[E = mc^{2}\]
Hyperlinks
Esto es un hyperlink: http://www.wikipedia.es
Pero es mejor ponerlo así: enlace a la Wikipedia
Notas al pie de página
Para poner notas al pie has de poner 1
Tablas
Habrá un tutorial especifico para tablas, pero lo más básico y sencillo es con la función kable():
knitr::kable(summary(iris))| Sepal.Length | Sepal.Width | Petal.Length | Petal.Width | Species | |
|---|---|---|---|---|---|
| Min. :4.300 | Min. :2.000 | Min. :1.000 | Min. :0.100 | setosa :50 | |
| 1st Qu.:5.100 | 1st Qu.:2.800 | 1st Qu.:1.600 | 1st Qu.:0.300 | versicolor:50 | |
| Median :5.800 | Median :3.000 | Median :4.350 | Median :1.300 | virginica :50 | |
| Mean :5.843 | Mean :3.057 | Mean :3.758 | Mean :1.199 | NA | |
| 3rd Qu.:6.400 | 3rd Qu.:3.300 | 3rd Qu.:5.100 | 3rd Qu.:1.800 | NA | |
| Max. :7.900 | Max. :4.400 | Max. :6.900 | Max. :2.500 | NA |
5. Aún más cosas
Hay más cosas, como referencias bibliográficas, pero ya será para el segundo curso de R.
A pesar de que Markdown permite formatear el texto, en ciertos sentidos es limitado; pero si quieres aún más flexibilidad en el output, tendrás que aprender html y mejor html+CSS.
Algunos ejemplos:
- Si quieres introducir un párrafo en otro color tendrás que hacerlo en html
Esto se mostrará en ROJO!!
- si quieres centrar un párrafo.
- Una pagina web
- Un vídeo
esto es una nota al pie.↩