Generando tablas de modelos y contrastes

tablas estadísticas

Post en proceso de elaboración

Autor

Pedro J. Pérez

Fecha de Publicación

5 de febrero de 2023

Fecha de modificación

25 de junio de 2023

1 Intro

De momento, hasta que tenga tiempo, simplemente voy a dejar 2 links a materiales sobre el tema que hice para un curso de introducción a la Ciencia de datos con R:

Tablas para modelos
Un poco de modelos. Una sección, dentro de un tutorial sobre EDA (análisis exploratorio de datos) sobre modelos estadísticos con R

Además, por si interesa el tema, voy a dejar aquí 3 scripts que he utilizado en clase para practicar sobre modelos y técnicas estadísticas:

Script (tablas y tablas para modelos)

#- script para seguir el tutorial nº 8: Tablas con .qmd y .Rmd
#- en el tutorial hay más contenido
library(tidyverse)
options(scipen = 999) #- para quitar la notación científica


#- cargamos datos para generar algunos resultados para mostrarlos en tablas
my_url <- "https://raw.githubusercontent.com/perezp44/iris_data/master/data/PIAAC_data_small.csv"
df_original <- read_csv(my_url)

df <- df_original %>% select(Country, Gender, Education, Wage_month, Wage_hour, Numeracy_score)

#- remotes::install_github("perezp44/pjpv2020.01")
df_aa <- pjpv2020.01::pjp_f_estadisticos_basicos(df) #- estadísticos básicos del df
df_bb <- pjpv2020.01::pjp_f_unique_values(df)        #- valores únicos de cada variable de df




#- tt_6 ------------------------------------------------------------------------
#- Calculamos la media, el mínimo, el máximo y la desviación típica de Wage_month
df_tt_6 <- df %>% 
    group_by(Country) %>% 
    summarise(W_medio  = mean(Wage_month, na.rm = TRUE) ,
              W_minimo = min(Wage_month, na.rm = TRUE)  ,
              W_maximo = max(Wage_month, na.rm = TRUE)  ,
              W_sd = sd(Wage_month, na.rm = TRUE) ) %>% 
              ungroup()


knitr::kable(df_tt_6)      #- funcionará ok en ficheros .qmd y .Rmd
knitr::kable(df_tt_6, format = "pipe") #- solo para verlo ahora


#- tt_9 ------------------------------------------------------------------------
#- ¿Cuanto más cobran los hombres (en %)?
df_tt_9 <- df %>% 
  group_by(Country, Gender) %>% 
  summarise(W_mes_medio = mean(Wage_month, na.rm = TRUE)) %>% 
  ungroup() %>% 
  pivot_wider(names_from = Gender, values_from = W_mes_medio) %>% 
  mutate(dif_W = Male-Female, dif_percent_W = dif_W/Female)

knitr::kable(df_tt_9, format = "pipe")


#- tt_10 ------------------------------------------------------------------------
#- Numeracy Score por país y nivel de estudios. La tabla nos va a salir alargada
df_tt_10 <- df %>% 
  group_by(Country, Education) %>% 
  summarise(Numeracy_media = mean(Numeracy_score, na.rm = TRUE)) %>% 
  ungroup() 

knitr::kable(df_tt_10, format = "pipe")


#- tt_11 -----------------------------------------------------------------------
#- Hagamos la anterior tabla más ancha (Una columna para cada país)
df_tt_11 <- df_tt_10 %>% 
            pivot_wider(names_from = Education, 
                        values_from = Numeracy_media)

knitr::kable(df_tt_11, format = "pipe")


#- quiero ordenar la tabla de menor a mayor nivel educativo
df <- df %>% 
  mutate(Education.f = forcats::as_factor(Education), .after = Education)

levels(df$Education.f)
df %>% count(Education)


#- renombrando los levels de los factores
df <- df %>% 
  mutate(Education.f = forcats::fct_recode(Education.f,
                    "Primaria"         = "Primary",
                    "Secundaria"       = "Secondary", 
                    "Secundaria_post"  = "Upper_second",
                    "Terciaria"        = "Tertiary" )) 
levels(df$Education.f)
df %>% count(Education.f)


#- reordenando los levels de un factor
df <- df %>% 
  mutate(Education.f = forcats::fct_relevel(Education.f, 
                                  "Primaria", "Secundaria", "Secundaria_post")) 
levels(df$Education.f)
df %>% count(Education.f)


#- tt_12 -----------------------------------------------------------------------
df_tt_12 <- df %>% 
  group_by(Country, Education.f) %>% 
  summarise(Numeracy_media = mean(Numeracy_score, na.rm = TRUE)) %>% 
  ungroup() %>%  
  pivot_wider(names_from = Education.f, 
              values_from = Numeracy_media) 

knitr::kable(df_tt_12, format = "pipe")

#- kable() más cosas -----------------------------------------------------------
#- kable() admite dar formato a algunos elementos de la tabla
knitr::kable(df_tt_12, format = "pipe",
             align = "c", 
             caption = "Numeracy Score por país",
             digits = 2, 
             format.args = list(decimal.mark = ",", big.mark = "."))


#- kableExtra pkg --------------------------------------------------------------
knitr::kable(df_tt_12) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))


#- scroll box
df_tt_12 %>%
  knitr::kable() %>%
  kableExtra::kable_styling(font_size = 11) %>%
  kableExtra::scroll_box(width = "50%", height = "60%")


#- colorines, girar ...
knitr::kable(df_tt_12) %>%
       kableExtra::kable_styling(full_width = F) %>%
       kableExtra::column_spec(1, bold = T, border_right = T) %>%
       kableExtra::column_spec(3, width = "20em", background = "yellow") %>% 
       kableExtra::row_spec(3:4, bold = T, color = "white", background = "#D7261E") %>% 
       kableExtra::row_spec(0, angle = 10)

no_quitar <- c("df", "df_tt_6")
rm(list=ls()[! ls() %in% no_quitar])

#- pkg gt: un ejemplo ----------------------------------------------------------
library(gt) #-remotes::install_github("rstudio/gt")

gt::gt(df_tt_6)


#- la tabla, al igual q un ggplot, se guarda en un objeto de clase "list"
gt_tbl <- df_tt_6 %>% gt()
gt_tbl

gt_tbl <- gt_tbl %>% tab_header(title = md("**Genero y nivel educativo**"),
                      subtitle = md("Porcentaje de *H y M* en cada nivel educativo"))
gt_tbl


#- pkg DT: un ejemplo ----------------------------------------------------------
DT::datatable(iris)
DT::datatable(df_tt_6)



DT::datatable(iris, filter = 'top', 
              options = list(pageLength = 5, autoWidth = TRUE ))

#- pkg reactable: un ejemplo ----------------------------------------------------------
reactable::reactable(iris)


#- pkg rpivotTable: un ejemplo ----------------------------------------------------------

library(rpivotTable) #- remotes::install_github(c("ramnathv/htmlwidgets", "smartinsightsfromdata/rpivotTable"))
rpivotTable(df, rows = "Gender", cols = c("Country"), width = "100%", height = "400px")


#- tablas para modelos ---------------------------------------------------------
urla <- "https://raw.githubusercontent.com/perezp44/iris_data/master/data/PIAAC_data_small.csv"
df <- read_csv(urla)


#- estimamos un modelo lineal
my_model <- lm(Wage_hour ~ Numeracy_score + Gender , data = df)
my_model

summary(my_model)


#- pkg stargazer ---------------------------------------------------------------
stargazer::stargazer(my_model, type = "html")
stargazer::stargazer(my_model, type = "text")


#- creo una variable dicotómica (para estimar un Logit)
df <- df %>% 
  mutate(Numeracy_score_b = 
                      ifelse(Numeracy_score > mean(Numeracy_score, na.rm = TRUE), 1, 0)) 

#- estimamos varios modelos y los almacenamos en una lista
my_models <- list()
my_models[['W:  OLS 1']]   <- lm( Wage_hour         ~ Numeracy_score + Gender , df)
my_models[['Nu: OLS 2']]   <- lm( Numeracy_score    ~ Education + Gender , df)
my_models[['Nu: Logit 1']] <- glm( Numeracy_score_b ~ Education + Gender , df, family = binomial())


stargazer::stargazer(my_models, type = "html", title = "Results", align = TRUE)
stargazer::stargazer(my_models, type = "text", title = "Results", align = TRUE)


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

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



#- más pkgs para tablas descriptivas -------------------------------------------

#- pkg janitor
df %>% janitor::tabyl(Education, Country, Gender)

Script (EDA)

#- ejemplo para ilustrar algunas ideas y pkgs para EDA
#- más cosas en el tutorial
library(tidyverse)

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

df_dicc <- pjpv.curso.R.2022::pjp_dicc(df)
df_uniques <- pjpv.curso.R.2022::pjp_valores_unicos(df)

#- INFORMES DESCRIPTIVOS -------------------------------------------------------
#- Hay pkgs q hacen informes descriptivos completos

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


#- ej: un nuevo paquete con una tabla chula: https://github.com/agstn/dataxray
#devtools::install_github("glin/reactable")
#devtools::install_github("agstn/dataxray")
library(dataxray)

df %>% make_xray() %>% view_xray()



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

#- co-ocurrencias de NA's
naniar::gg_miss_upset(df)


#- ¿QUÉ HACEMOS con los NA's? --------------------------------------------------
#- quitamos las filas q tengan NA en cualquiera de las variables/columnas
zz1 <- df %>% tidyr::drop_na()

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


#- SELECCIONAR las vv. NUMÉRICAS -----------------------------------------------
df_numeric <- df %>% select_if(is.numeric)      #- antigua sintaxis
df_numeric <- df %>% select(where(is.numeric))  #- nueva API

summarytools::descr(df)

DataExplorer::plot_histogram(df, ncol = 2)
DataExplorer::plot_density(df, ncol = 2)

corrr::correlate(df_numeric)

#- tabla con matriz de correlaciones, con gt::gt()
df_numeric %>% corrr::correlate() %>%
  gt::gt() %>%
  gt::fmt_number(decimals = 1, sep_mark = ".",  dec_mark = ",")

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

#- https://r-coder.com/correlation-plot-r/
#- https://albert-rapp.de/posts/ggplot2-tips/13_alternative_corrplots/13_alternative_corrplots.html


#- Boxplots ------------------
#DataExplorer::plot_boxplot(df, by = "estudios_madre.ff")
DataExplorer::plot_boxplot(df, by = "parto_cesarea.f")

#- boxplots (3 v.)
df %>% explore::explore(edad_madre.1, estudios_madre.ff, target = parto_cesarea.f)

#- explore_all()
df %>% select(sexo_bebe.f, edad_madre.1, peso_bebe,  estudios_madre.ff, parto_cesarea.f) %>%
  explore::explore_all(target = parto_cesarea.f)

#- janitor:: tabulaciones cruzadas
df %>% janitor::tabyl(estudios_madre.ff, parto_cesarea.f) %>%
       janitor::adorn_percentages() %>%
       gt::gt()


#- VARIABLES CATEGÓRICAS -------------------------------------------------------
#- no pasa nada si no lo entendéis del todo (siempre se puede hacer "a mano")
df %>% select(where(is.factor)) %>% names()     #- se nos escapan las q son character.

vv_numericas <- df %>% select(where(is.numeric)) %>% names()
vv_numericas

zz <- df %>% select(!vv_numericas)         #- deprecated but funciona
names(zz)
zz <- df %>% select(!all_of(vv_numericas)) #- forma correcta


#- porcentajes
inspectdf::inspect_cat(zz) %>% inspectdf::show_plot(high_cardinality = 1)

#- graficos de barras
DataExplorer::plot_bar(zz)
DataExplorer::plot_bar(zz, by = "sexo_bebe.f")
DataExplorer::plot_bar(zz, by = "parto_cesarea.f")
DataExplorer::plot_bar(zz, by = "estudios_madre.ff")



#- TABLAS ----------------------------------------------------------------------
summarytools::freq(df$sexo_bebe.f, style = "rmarkdown")

#- cross-tabulations
summarytools::ctable(df$sexo_bebe.f, df$parto_cesarea.f)
summarytools::ctable(df$sexo_bebe.f, df$parto_cesarea.f, chisq = TRUE)


#- JANITOR
df %>% janitor::tabyl(parto_cesarea.f) %>% gt::gt()

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

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

#- ademas, como janitor almacena las tablas en df, podemos usar gt()
df %>% janitor::tabyl(parto_cesarea.f, estudios_madre.ff) %>% gt::gt()


#- CONTRASTES ------------------------------------------------------------------
t.test(df$peso_bebe, mu = 3250)
t.test(df$peso_bebe ~ df$parto_cesarea.f)


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

#- frente a todas las v. del df
mod_1 <- lm(peso_bebe ~ . , data = df_m)
summary(mod_1)


mod_1 <- lm(peso_bebe ~ parto_nn_semanas + edad_madre.1 + sexo_bebe.f , data = df_m)

summary(mod_1)                   #- tabla resumen

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

sjPlot::tab_model(m1)

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

sjPlot::tab_model(m1, m2, show.ci   = FALSE)


#- mas tablas de modelos
stargazer::stargazer(m1, type = "text")  #- latex. html

stargazer::stargazer(m1, m2,  type = "text")  #- latex. html


library(apaTables) #- install.packages("apaTables",dep = TRUE)
apaTables::apa.reg.table(m1)

apaTables::
#- en el tutorial de tablas vimos más posibilidades

Script (easystats package)

#- easystats: (https://easystats.github.io/easystats/)
#- An R Framework for Easy Statistical Modeling, Visualization, and Reporting
#- install.packages("easystats")

#- report: https://easystats.github.io/report/ ----------------
library(tidyverse)
library(easystats)

report(iris)

iris %>% select(-starts_with("Sepal")) %>% group_by(Species) %>%
  report() %>% summary()

report(t.test(mtcars$mpg ~ mtcars$am))

cor.test(iris$Sepal.Length, iris$Sepal.Width) %>%
  report() %>%  as.data.frame()

aov(Sepal.Length ~ Species, data = iris) %>% report()


model <- lm(Sepal.Length ~ Species, data = iris)
summary(model)
report(model)

model <- glm(vs ~ mpg * drat, data = mtcars, family = "binomial")
summary(model)
report(model)


library(lme4)
model <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
summary(model)
report(model)


#- parameters: https://easystats.github.io/parameters/

model <- lm(Sepal.Width ~ Petal.Length * Species + Petal.Width, data = iris)
model_parameters(model)
model_parameters(model, standardize = "refit")


#- see: https://easystats.github.io/see/

results <- summary(correlation(iris))
plot(results, show_data = "points")


model <- lm(wt ~ am * cyl, data = mtcars)
summary(model)
plot(parameters(model))



model <- lm(mpg ~ factor(cyl) * wt, data = mtcars)
summary(model)
results <- fortify(model)

# step-3
ggplot(results) +
  geom_point(aes(x = wt, y = mpg, color = `factor(cyl)`)) +
  geom_line(aes(x = wt, y = .fitted, color = `factor(cyl)`))

Extensión: Algunas referencias para cuando retome el tema

Esto puede ayudar a escribir ecuaciones/modelos estimados
Quarto con gt de A. Rapp
Making APA Tables with the gt Package
Tablas de MRLM con pkg finalfit
paquete rempsyc: Convenience functions for psychology. Hace tablas APA
Tablas side by side
Copy Table in Excel and Paste as a Markdown Table

Reutilizar

https://creativecommons.org/licenses/by/4.0/deed.es