Regresión Multinomial
Taller de R - Bloque 05/Random Forest: Variables Numericas y Categoricas
Espacio Colectivo Mx
2024
Introducción
Con la intención de generar un análisis reproducible, fácil y entendible, se elabora el siguiente código para estudiar los modelos multinomiales. Se recupera el código de las sesiones de Julia Silge.
Sin antes, se ha de mencionar que, al ser un análisis sencillo no se profundiza teóricamente el diseño de muestra y en los supuestos de los modelos, por lo que es tarea de la persona lectora de este código buscar la información complementaría.
Se agradece cualquier tipo de comentario y se advierte que el código presente se mantiene en actualización.
Putos clave para la discusión de un modelo
- Objetivo del Modelo
- Estructura del Modelo
- Supuesto del Modelo
- Ajuste del Modelo
- Selección del Modelo
Librerías de trabajo
A continuación, se enlistas las librerías a utilizar, mismas que serán mencionadas a lo largo del código para poder identificar el uso que tienen.
El siguiente código instala las paqueterías necesarias en caso de no contar previamente con ellas y posteriormente las carga.
if(!require('pacman')) install.packages('pacman')
pacman::p_load(tidyverse,#manipulación de datos
janitor,#limpieza y tablas de frecuencia
tidymodels,#modelos
themis,
vip)#Notación cientifica
options(scipen = 999)Descargamos datos
Los datos provienen de las bases de tidytuesday, para este caso la data corresponde a información sobre volvanes en el mundo.
#Cargamos datos
volcano_raw <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-05-12/volcano.csv")Objetivos del Modelo
A partir de la variable tipo de volcanes, se construye un modelo multiclase con random forest classifier para predecir el tipo de volcano basada en ciertas caracteristicas.
Estructura de Datos
#Primeros cinco datos
head(volcano_raw, 3)## # A tibble: 3 × 26
## volcano_number volcano_name primary_volcano_type last_eruption_year country
## <dbl> <chr> <chr> <chr> <chr>
## 1 283001 Abu Shield(s) -6850 Japan
## 2 355096 Acamarachi Stratovolcano Unknown Chile
## 3 342080 Acatenango Stratovolcano(es) 1972 Guatemala
## # ℹ 21 more variables: region <chr>, subregion <chr>, latitude <dbl>,
## # longitude <dbl>, elevation <dbl>, tectonic_settings <chr>,
## # evidence_category <chr>, major_rock_1 <chr>, major_rock_2 <chr>,
## # major_rock_3 <chr>, major_rock_4 <chr>, major_rock_5 <chr>,
## # minor_rock_1 <chr>, minor_rock_2 <chr>, minor_rock_3 <chr>,
## # minor_rock_4 <chr>, minor_rock_5 <chr>, population_within_5_km <dbl>,
## # population_within_10_km <dbl>, population_within_30_km <dbl>, …#Estructura de datos
glimpse(volcano_raw)## Rows: 958
## Columns: 26
## $ volcano_number <dbl> 283001, 355096, 342080, 213004, 321040, 28317…
## $ volcano_name <chr> "Abu", "Acamarachi", "Acatenango", "Acigol-Ne…
## $ primary_volcano_type <chr> "Shield(s)", "Stratovolcano", "Stratovolcano(…
## $ last_eruption_year <chr> "-6850", "Unknown", "1972", "-2080", "950", "…
## $ country <chr> "Japan", "Chile", "Guatemala", "Turkey", "Uni…
## $ region <chr> "Japan, Taiwan, Marianas", "South America", "…
## $ subregion <chr> "Honshu", "Northern Chile, Bolivia and Argent…
## $ latitude <dbl> 34.500, -23.292, 14.501, 38.537, 46.206, 37.6…
## $ longitude <dbl> 131.600, -67.618, -90.876, 34.621, -121.490, …
## $ elevation <dbl> 641, 6023, 3976, 1683, 3742, 1728, 1733, 1250…
## $ tectonic_settings <chr> "Subduction zone / Continental crust (>25 km)…
## $ evidence_category <chr> "Eruption Dated", "Evidence Credible", "Erupt…
## $ major_rock_1 <chr> "Andesite / Basaltic Andesite", "Dacite", "An…
## $ major_rock_2 <chr> "Basalt / Picro-Basalt", "Andesite / Basaltic…
## $ major_rock_3 <chr> "Dacite", " ", " ", "Basalt / Picro-Basalt", …
## $ major_rock_4 <chr> " ", " ", " ", "Andesite / Basaltic Andesite"…
## $ major_rock_5 <chr> " ", " ", " ", " ", " ", " ", " ", " ", " ", …
## $ minor_rock_1 <chr> " ", " ", "Basalt / Picro-Basalt", " ", "Daci…
## $ minor_rock_2 <chr> " ", " ", " ", " ", " ", "Basalt / Picro-Basa…
## $ minor_rock_3 <chr> " ", " ", " ", " ", " ", " ", " ", "Andesite …
## $ minor_rock_4 <chr> " ", " ", " ", " ", " ", " ", " ", " ", " ", …
## $ minor_rock_5 <chr> " ", " ", " ", " ", " ", " ", " ", " ", " ", …
## $ population_within_5_km <dbl> 3597, 0, 4329, 127863, 0, 428, 101, 51, 0, 98…
## $ population_within_10_km <dbl> 9594, 7, 60730, 127863, 70, 3936, 485, 6042, …
## $ population_within_30_km <dbl> 117805, 294, 1042836, 218469, 4019, 717078, 1…
## $ population_within_100_km <dbl> 4071152, 9092, 7634778, 2253483, 393303, 5024…#Variable dependiente
volcano_raw %>%
count(primary_volcano_type,#datos a contar
sort = TRUE)#ordenar datos## # A tibble: 26 × 2
## primary_volcano_type n
## <chr> <int>
## 1 Stratovolcano 353
## 2 Stratovolcano(es) 107
## 3 Shield 85
## 4 Volcanic field 71
## 5 Pyroclastic cone(s) 70
## 6 Caldera 65
## 7 Complex 46
## 8 Shield(s) 33
## 9 Submarine 27
## 10 Lava dome(s) 26
## # ℹ 16 more rows#NA's
volcano_raw %>%
summarise(across(everything(), ~ sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "variable", values_to = "na_count")## # A tibble: 26 × 2
## variable na_count
## <chr> <int>
## 1 volcano_number 0
## 2 volcano_name 0
## 3 primary_volcano_type 0
## 4 last_eruption_year 0
## 5 country 0
## 6 region 0
## 7 subregion 0
## 8 latitude 0
## 9 longitude 0
## 10 elevation 0
## # ℹ 16 more rows- Los valores faltantes no estan registrados como NA’s
- Se cuenta con variables numericas y de tipo caracter
- Existen diversas clases para la variable dependiente b
dt.volcano <- volcano_raw %>%
#Creamos nuevo conjunto de datos con el conjunto volcano_raw
#Reclasificamos nuestra variable de interes
transmute(volcanotype = case_when(str_detect(primary_volcano_type,"Stratovolcano") ~ "Stratovolcano",
str_detect(primary_volcano_type, "Shield") ~ "Shield",
TRUE ~ "Other"),
volcano_number, latitude, longitude, elevation, tectonic_settings, major_rock_1) %>%
#convertimos a factor
mutate_if(is.character,factor)La falta de datos para dividir la muestra en datos de entrenamiento entorpece el modelo, por lo que realizamos un remuestreo con boostraps, con la cual entrenaremos los datos y posteriormente constrastaremos. El paquete boot está asociado al libro Bootstrap Methods and Their Applications (Davison and Hinkley (Bootstrap Methods and Their Application, 1997)), entre las ventajas de los procedimientos bootstrap:
- No necesitamos la condición de normalidad en la que se basan las técnicas de inferencia paramétricas.
- En la mayoría de casos, la distribución bootstrap de un estadístico tiene aproximadamente la misma forma y dispersión que la distribución del estadístico que se obtiene al seleccionar muchas muestras de la población, pero su centro es el valor que toma el estadístico en la muestra original, en lugar del verdadero valor del parámetro que intentamos estimar.
- La técnica bootstrap nos permite calcular el error típico de cualquier estadístico sin necesidad de recurrir a la teoría estadística.
#library(tidymodels)
volcano_boot <- bootstraps(dt.volcano)
#Aceso a los datos
#La columna splits tiene información de las muestras seleccionadas
as.tibble(volcano_boot$splits[[1]])## # A tibble: 958 × 7
## volcanotype volcano_number latitude longitude elevation tectonic_settings
## <fct> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 Stratovolcano 263100 -7.32 108. 2665 Subduction zone / …
## 2 Other 241010 -35.3 174. 388 Subduction zone / …
## 3 Other 371050 64.1 -21.4 803 Rift zone / Oceani…
## 4 Stratovolcano 263040 -6.75 107. 1699 Subduction zone / …
## 5 Other 390050 -65.0 -60.0 368 Intraplate / Conti…
## 6 Other 290010 43.8 146. 535 Subduction zone / …
## 7 Stratovolcano 263270 -7.81 112. 2563 Subduction zone / …
## 8 Stratovolcano 372050 63.9 -19.1 1280 Rift zone / Oceani…
## 9 Stratovolcano 221110 13.1 40.9 1250 Rift zone / Interm…
## 10 Stratovolcano 300270 56.7 161. 3283 Subduction zone / …
## # ℹ 948 more rows
## # ℹ 1 more variable: major_rock_1 <fct>Balance de datos para las clases
#Frecuencias
tabyl(dt.volcano$volcanotype)## dt.volcano$volcanotype n percent
## Other 379 0.3956159
## Shield 118 0.1231733
## Stratovolcano 461 0.4812109Hay un problema de datos balanceadtos.
#library(themis)
volcano.recipe <- recipe(volcanotype ~ ., data = dt.volcano) %>%
#Mantenemos volcano_number aunque no se usara
update_role(volcano_number, new_role = "Id") %>%
#agrupamos y generamos nuevas categorias con las de menos representación
step_other(tectonic_settings) %>%
step_other(major_rock_1) %>%
#Variables dummy
step_dummy(tectonic_settings, major_rock_1) %>%
#Eliminamos variables sin varianza
step_zv(all_predictors()) %>%
#Normalizamos variables
step_normalize(all_predictors()) %>%
#Balancea datos
step_smote(volcanotype)
#Aplicar prepraciones
volcano.prep <- prep(volcano.recipe)
#Resultados del cocinamiento
juice(volcano.prep) %>%
count(volcanotype)## # A tibble: 3 × 2
## volcanotype n
## <fct> <int>
## 1 Other 461
## 2 Shield 461
## 3 Stratovolcano 461rf_spec <- rand_forest(trees = 1000) %>%
#Queremos clasificación
set_mode("classification") %>%
#Multinomial
set_engine("ranger")volcano.wf <- workflow() %>%
add_recipe(volcano.recipe) %>%
add_model(rf_spec)
volcano.wf## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 6 Recipe Steps
##
## • step_other()
## • step_other()
## • step_dummy()
## • step_zv()
## • step_normalize()
## • step_smote()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
##
## Main Arguments:
## trees = 1000
##
## Computational engine: rangervolcano.results <- fit_resamples(
volcano.wf,
resamples = volcano_boot,
control = control_resamples(save_pred = TRUE,
verbose = TRUE)
)volcano.results %>%
collect_metrics()## # A tibble: 3 × 6
## .metric .estimator mean n std_err .config
## <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 accuracy multiclass 0.645 25 0.00479 Preprocessor1_Model1
## 2 brier_class multiclass 0.242 25 0.00156 Preprocessor1_Model1
## 3 roc_auc hand_till 0.788 25 0.00417 Preprocessor1_Model1volcano.results %>%
#Valores predichos
collect_predictions() %>%
#Matriz de confuciión
conf_mat(volcanotype, .pred_class)## Truth
## Prediction Other Shield Stratovolcano
## Other 1944 307 834
## Shield 341 570 266
## Stratovolcano 1221 171 3189#library(vip)
#Impacto en el modelo por variable
rf_spec %>%
set_engine("ranger", importance = "permutation") %>%
fit(
volcanotype ~ .,
data = juice(volcano.prep) %>%
janitor::clean_names() %>%
dplyr::select(-volcano_number)
) %>%
vip(geom = "point")
- Las variables que mayormente influyen en el modelo son la longitud, latitud y longitud
volcano.results %>%
collect_predictions() %>%
mutate(correcto = volcanotype == .pred_class) %>%
count(correcto) %>%
adorn_totals() %>%
adorn_percentages("col")## correcto n
## FALSE 0.3550831
## TRUE 0.6449169
## Total 1.0000000volcano.results %>%
collect_predictions() %>%
mutate(correcto = volcanotype == .pred_class) %>%
left_join(dt.volcano %>%
mutate(.row = row_number()))## # A tibble: 8,843 × 15
## .pred_class .pred_Other .pred_Shield .pred_Stratovolcano id .row
## <fct> <dbl> <dbl> <dbl> <chr> <int>
## 1 Other 0.487 0.229 0.285 Bootstrap01 1
## 2 Stratovolcano 0.191 0.128 0.681 Bootstrap01 6
## 3 Stratovolcano 0.456 0.0417 0.502 Bootstrap01 13
## 4 Other 0.651 0.0125 0.336 Bootstrap01 14
## 5 Other 0.600 0.144 0.256 Bootstrap01 15
## 6 Stratovolcano 0.208 0.255 0.537 Bootstrap01 17
## 7 Stratovolcano 0.429 0.0703 0.501 Bootstrap01 20
## 8 Other 0.398 0.239 0.363 Bootstrap01 26
## 9 Shield 0.346 0.500 0.154 Bootstrap01 29
## 10 Shield 0.0796 0.806 0.115 Bootstrap01 31
## # ℹ 8,833 more rows
## # ℹ 9 more variables: volcanotype <fct>, .config <chr>, correcto <lgl>,
## # volcano_number <dbl>, latitude <dbl>, longitude <dbl>, elevation <dbl>,
## # tectonic_settings <fct>, major_rock_1 <fct>