#- script para clase_13: algunas cosas sobre  EDA y modelos
#- más cosas en el tutorial
library(tidyverse)
options(scipen = 999) #- quitar notación científica

#- DATOS (de bebes) ------------------------------------------------------------
url_bebes <- "https://github.com/perezp44/archivos_download.2/raw/main/df_bebes_EDA.rds"
df_orig <- rio::import(url_bebes)

#- una primera mirada a los datos
df_aa <- pjpv.curso.R.2022::pjp_dicc(df_orig)
#- valores únicos de las variables
df_bb <- pjpv.curso.R.2022::pjp_valores_unicos(df_orig)
df_cc <- pjpv.curso.R.2022::pjp_unique_values(df_orig)


#- vamos a seleccionar una muestra reducida; p.ej 1000 bebes (para ir rápido)
set.seed(1234) #- for reproducibility, para q nos salgan a todos mismos rtdos
df <- df_orig %>% slice_sample(n = 1000)

#- no quitar -------------------------------------------------------------------
no_quitar <- c("df_orig", "df", "df_aa")
rm(list=ls()[! ls() %in% no_quitar])




#- MODELOS ---------------------------------------------------------------------
df_m <- df %>% select(peso_bebe, parto_nn_semanas, sexo_bebe.f, parto_cesarea.f, edad_madre.1, estudios_madre.ff)
#df_m <- df_m %>% drop_na()


#- estimamos varios modelos

m1 <- lm(peso_bebe ~ 
           parto_nn_semanas + sexo_bebe.f + edad_madre.1 , 
         data = df_m)
m2 <- lm(peso_bebe ~ 
           parto_nn_semanas + sexo_bebe.f + edad_madre.1 + edad_padre.1 + estudios_madre.ff , 
         data = df)

#- vamos a estimar un logit
df_m <- df_m %>% 
  mutate(peso_bebe.f = ifelse(peso_bebe > mean(peso_bebe, na.rm = TRUE), 1, 0 ))

m3 <- glm(peso_bebe.f ~ parto_nn_semanas + sexo_bebe.f , data = df_m, family = binomial(link = "logit"))


mys_modelitos <- list("OLS 1" = m1, "OLS 2" = m2, "Logit 1" = m3)



#- Tablas para modelos ---------------------------------------------------------

#- con pkg stargazer
stargazer::stargazer(m1)                 #- por defecto en Latex
stargazer::stargazer(m1, type = "html")
stargazer::stargazer(m1, type = "text")
stargazer::stargazer(mys_modelitos, type = "text")


#- con sjPlot package
sjPlot::tab_model(m1)
sjPlot::tab_model(mys_modelitos)


sjPlot::plot_model(m1, sort.est = TRUE)


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

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


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


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



#- cocoon pkg
#- https://jeffreyrstevens.github.io/cocoon/articles/cocoon.html
