# Instalar si es necesario
paquetes <- c("tidyverse", "tidymodels", "ranger", "vip", "pdp", "xgboost")

# require(): intenta cargar el paquete y devuelve TRUE/FALSE
# Si no está instalado (FALSE), lo instala y carga
for (pkg in paquetes) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}

# Cargar tidymodels (carga varios paquetes a la vez:
# rsample, parsnip, recipes, workflows, tune, yardstick, etc.)
library(tidymodels)
library(ranger)      # Motor rápido para Random Forest
library(vip)         # Variable Importance Plots
library(pdp)         # Partial Dependence Plots

# Configurar semilla para reproducibilidad
set.seed(2026)

# Cargar el dataset (generado con datos/crear-datos.R)
datos <- read_csv("datos/latinobarometro_sim.csv", show_col_types = FALSE)

# Convertir variables categóricas a factores
# El primer nivel de voto es la clase positiva para yardstick
datos <- datos |>
  mutate(
    pais = factor(pais),
    zona = factor(zona),
    genero = factor(genero),
    uso_internet = factor(uso_internet, levels = c("nunca", "semanal", "diario")),
    voto = factor(voto, levels = c("si", "no"))
  )

# Estructura y distribución del voto
glimpse(datos)

datos |> count(voto) |> mutate(proporcion = n / sum(n))

# initial_split(): divide los datos aleatoriamente en train y test
# prop = 0.75: 75% para entrenamiento, 25% para test
# strata = voto: estratificar para mantener la proporción de clases
division <- initial_split(datos, prop = 0.75, strata = voto)

datos_train <- training(division)
datos_test <- testing(division)

# Verificar proporciones
cat("Proporción en train:\n")

prop.table(table(datos_train$voto))

cat("\nProporción en test:\n")

prop.table(table(datos_test$voto))

# recipe(): define el preprocesamiento como una "receta de cocina"
# voto ~ ...: la fórmula indica variable objetivo ~ predictores
receta <- recipe(voto ~ edad + educacion_anios + ingreso_hogar + zona +
                 genero + confianza_gobierno + confianza_justicia +
                 satisfaccion_democracia + percepcion_economia +
                 uso_internet + interes_politica,
                 data = datos_train) |>
  # step_dummy(): convierte categóricas a variables indicadoras (0/1)
  step_dummy(all_nominal_predictors()) |>
  # step_normalize(): centra (media = 0) y escala (sd = 1) las numéricas
  step_normalize(all_numeric_predictors()) |>
  # step_zv(): elimina variables con varianza cero (constantes)
  step_zv(all_predictors())

# prep(): estima los parámetros de la receta (ej: medias para normalizar)
# juice(): aplica la receta y devuelve los datos transformados
receta |> prep() |> juice() |> glimpse()

# logistic_reg(): modelo de regresión logística
# set_engine("glm"): usar el motor glm de R base
# set_mode("classification"): tarea de clasificación (no regresión)
modelo_logit <- logistic_reg() |>
  set_engine("glm") |>
  set_mode("classification")

# workflow(): combina preprocesamiento (receta) + modelo en un solo objeto
wf_logit <- workflow() |>
  add_recipe(receta) |>    # agregar la receta de preprocesamiento
  add_model(modelo_logit)  # agregar el modelo

# fit(): ajustar el workflow completo a los datos de entrenamiento
ajuste_logit <- fit(wf_logit, data = datos_train)

# predict(): genera predicciones de clase ("si"/"no")
# type = "prob": genera probabilidades para cada clase (.pred_no, .pred_si)
# bind_cols(): une las columnas de predicciones con los datos reales
pred_logit <- predict(ajuste_logit, datos_test) |>
  bind_cols(predict(ajuste_logit, datos_test, type = "prob")) |>
  bind_cols(datos_test |> select(voto))

# metrics(): calcula múltiples métricas de evaluación
# truth: variable real, estimate: predicción de clase, .pred_si: probabilidades
metricas_logit <- pred_logit |>
  metrics(truth = voto, estimate = .pred_class, .pred_si)

metricas_logit

# conf_mat(): construye la tabla de predicho vs. real
# autoplot(type = "heatmap"): visualización como mapa de calor
conf_mat(pred_logit, truth = voto, estimate = .pred_class) |>
  autoplot(type = "heatmap") +
  scale_fill_gradient(low = "white", high = "#2d4563") +
  labs(title = "Matriz de confusión - Regresión logística")

# rand_forest(): modelo de Random Forest
# tune(): marcador especial que indica "optimizar este valor automáticamente"
modelo_rf_tune <- rand_forest(
  mtry = tune(),       # Número de variables a considerar en cada split
  trees = 500,         # Número de árboles (fijo, no se optimiza)
  min_n = tune()       # Mínimo de observaciones en nodo terminal
) |>
  # importance = "impurity": calcular importancia con reducción de impureza (Gini)
  set_engine("ranger", importance = "impurity") |>
  set_mode("classification")

# workflow(): combina receta + modelo
wf_rf_tune <- workflow() |>
  add_recipe(receta) |>
  add_model(modelo_rf_tune)

# extract_parameter_set_dials(): muestra los hiperparámetros marcados con tune()
modelo_rf_tune |> extract_parameter_set_dials()

# grid_regular(): crea una grilla con valores uniformemente espaciados
# Cada parámetro tiene un rango definido; levels indica cuántos valores probar
grilla_rf <- grid_regular(
  mtry(range = c(2, 8)),    # De 2 a 8 variables por split
  min_n(range = c(5, 30)),  # De 5 a 30 obs mínimas por nodo
  levels = c(4, 4)          # 4 valores de cada uno = 4 x 4 = 16 combinaciones
)

# Ver la grilla
grilla_rf

# Alternativa: grilla aleatoria (más eficiente para muchos hiperparámetros)
# grilla_rf <- grid_random(
#   mtry(range = c(2, 8)),
#   min_n(range = c(5, 30)),
#   size = 20  # 20 combinaciones aleatorias
# )

# vfold_cv(): divide los datos en v grupos (folds) para validación cruzada
# strata = voto: estratificar para que cada fold mantenga la proporción de clases
folds <- vfold_cv(datos_train, v = 5, strata = voto)

# Ver los folds
folds

# Cada fold tiene ~300 obs en training y ~75 en assessment

# tune_grid(): evalúa cada combinación de hiperparámetros con CV
# Para cada combinación de la grilla, entrena el modelo en cada fold
# y calcula las métricas especificadas
resultados_tune <- tune_grid(
  wf_rf_tune,               # el workflow con tune() pendientes
  resamples = folds,         # los folds de validación cruzada
  grid = grilla_rf,          # la grilla de combinaciones a probar
  # metric_set(): define qué métricas calcular en cada evaluación
  metrics = metric_set(accuracy, roc_auc, f_meas),
  control = control_grid(verbose = FALSE)  # no imprimir progreso
)

# collect_metrics(): extrae los resultados promediados de todos los folds
resultados_tune |>
  collect_metrics() |>
  filter(.metric == "roc_auc") |>
  arrange(desc(mean))

# autoplot(): método genérico que sabe graficar objetos de tidymodels
autoplot(resultados_tune) +
  theme_minimal() +
  labs(title = "Resultados del tuning de Random Forest")

# select_best(): elige la combinación con el mejor valor de la métrica
mejor_rf <- select_best(resultados_tune, metric = "roc_auc")
mejor_rf

# finalize_workflow(): reemplaza los tune() por los valores óptimos
# Luego ajustamos con todos los datos de entrenamiento
wf_rf_final <- finalize_workflow(wf_rf_tune, mejor_rf)
ajuste_rf <- fit(wf_rf_final, data = datos_train)

# Predicciones y métricas en test
pred_rf <- predict(ajuste_rf, datos_test) |>
  bind_cols(predict(ajuste_rf, datos_test, type = "prob")) |>
  bind_cols(datos_test |> select(voto))

metricas_rf <- pred_rf |>
  metrics(truth = voto, estimate = .pred_class, .pred_si)

cat("Regresión logística:\n"); print(metricas_logit)

cat("\nRandom Forest:\n");     print(metricas_rf)

# roc_curve(): calcula sensibilidad y especificidad para cada umbral
# truth: la variable real, .pred_si: probabilidades predichas
roc_logit <- pred_logit |>
  roc_curve(truth = voto, .pred_si) |>
  mutate(modelo = "Regresión logística")

roc_rf <- pred_rf |>
  roc_curve(truth = voto, .pred_si) |>
  mutate(modelo = "Random Forest")

# Combinar y graficar
bind_rows(roc_logit, roc_rf) |>
  ggplot(aes(x = 1 - specificity, y = sensitivity, color = modelo)) +
  geom_path(linewidth = 1.2) +
  geom_abline(linetype = "dashed", color = "gray50") +
  coord_equal() +
  labs(title = "Comparación de curvas ROC",
       x = "1 - Especificidad (Tasa de falsos positivos)",
       y = "Sensibilidad (Tasa de verdaderos positivos)",
       color = "Modelo") +
  theme_minimal()

# extract_fit_parsnip(): extrae el modelo ajustado del workflow
# (devuelve un objeto parsnip, no el objeto nativo del motor)
modelo_extraido <- extract_fit_parsnip(ajuste_rf)

# vip(): gráfico de importancia de variables
# num_features: cuántas variables mostrar (las más importantes)
vip(modelo_extraido, num_features = 15) +
  labs(title = "Importancia de variables (Gini)",
       subtitle = "Random Forest para predicción de voto") +
  theme_minimal()

# extract_fit_engine(): extrae el objeto nativo del motor (ranger)
modelo_ranger <- extract_fit_engine(ajuste_rf)

# bake(): aplica la receta ya estimada a los datos de entrenamiento
datos_prep <- bake(prep(receta), new_data = datos_train)

# partial() calcula el efecto marginal sobre P(voto = "si")
# which.class = 1: primera clase del factor, es decir, "si"
pdp_edad <- partial(modelo_ranger, pred.var = "edad",
                    train = datos_prep, prob = TRUE, which.class = 1)

pdp_interes <- partial(modelo_ranger, pred.var = "interes_politica",
                       train = datos_prep, prob = TRUE, which.class = 1)

# Graficar lado a lado con patchwork (o cowplot)
library(patchwork)
p1 <- autoplot(pdp_edad) +
  labs(title = "PDP: Edad", x = "Edad (normalizada)",
       y = "P(voto = sí)") + theme_minimal()

p2 <- autoplot(pdp_interes) +
  labs(title = "PDP: Interés en política", x = "Interés (normalizado)",
       y = "P(voto = sí)") + theme_minimal()

p1 + p2

# Probar thresholds de 0.3 a 0.7 y guardar los resultados
thresholds <- seq(0.3, 0.7, by = 0.05)
resultados <- data.frame(threshold = numeric(), f1 = numeric())

for (t in thresholds) {
  pred_nuevo <- pred_logit |>
    mutate(.pred_class_nuevo = factor(
      ifelse(.pred_si > t, "si", "no"),
      levels = c("si", "no")))

  f1 <- f_meas(pred_nuevo, truth = voto,
               estimate = .pred_class_nuevo)
  resultados <- rbind(resultados,
                      data.frame(threshold = t, f1 = f1$.estimate))
}

# Ver todos los resultados, ordenados de mejor a peor
resultados |> arrange(desc(f1))

# Visualizar
ggplot(resultados, aes(x = threshold, y = f1)) +
  geom_line(linewidth = 1.2, color = "#2d4563") +
  geom_point(size = 3, color = "#2d4563") +
  labs(title = "F1-score según el threshold de clasificación",
       x = "Threshold", y = "F1-score") +
  theme_minimal()

# Filtrar datos de Uruguay
datos_uruguay <- datos |> filter(pais == "Uruguay")
cat("Observaciones en Uruguay:", nrow(datos_uruguay), "\n")

# Con pocos datos (~28 obs), no dividimos en train/test
# Entrenamos con todo el subconjunto solo para comparar VIP
receta_uy <- recipe(voto ~ edad + educacion_anios + ingreso_hogar + zona +
                    genero + confianza_gobierno + confianza_justicia +
                    satisfaccion_democracia + percepcion_economia +
                    uso_internet + interes_politica,
                    data = datos_uruguay) |>
  step_dummy(all_nominal_predictors()) |>
  step_normalize(all_numeric_predictors()) |>
  step_zv(all_predictors())

modelo_rf_uy <- rand_forest(trees = 500, mtry = 4, min_n = 3) |>
  set_engine("ranger", importance = "impurity") |>
  set_mode("classification")

ajuste_rf_uy <- workflow() |>
  add_recipe(receta_uy) |>
  add_model(modelo_rf_uy) |>
  fit(data = datos_uruguay)

# Comparar importancia de variables
vip(extract_fit_parsnip(ajuste_rf_uy), num_features = 10) +
  labs(title = "Importancia de variables - solo Uruguay",
       subtitle = paste0("n = ", nrow(datos_uruguay), " observaciones")) +
  theme_minimal()

# Modelo de árbol con costo de complejidad a ajustar
modelo_arbol <- decision_tree(
  cost_complexity = tune(),
  tree_depth = 10,
  min_n = 10
) |>
  set_engine("rpart") |>
  set_mode("classification")

# Workflow con la misma receta
wf_arbol <- workflow() |>
  add_recipe(receta) |>
  add_model(modelo_arbol)

# Grilla de búsqueda para cost_complexity
# range en escala log10: 10^-4 = 0.0001 a 10^-1 = 0.1
grilla_arbol <- grid_regular(
  cost_complexity(range = c(-4, -1)),
  levels = 10
)

# Tuning con validación cruzada (reutilizamos folds)
resultados_arbol <- tune_grid(
  wf_arbol,
  resamples = folds,
  grid = grilla_arbol,
  metrics = metric_set(roc_auc),
  control = control_grid(verbose = FALSE)
)

# Mejor modelo
mejor_arbol <- select_best(resultados_arbol, metric = "roc_auc")
cat("Mejor cost_complexity:", mejor_arbol$cost_complexity, "\n")

# Ajustar y evaluar
wf_arbol_final <- finalize_workflow(wf_arbol, mejor_arbol)
ajuste_arbol <- fit(wf_arbol_final, data = datos_train)

pred_arbol <- predict(ajuste_arbol, datos_test) |>
  bind_cols(predict(ajuste_arbol, datos_test, type = "prob")) |>
  bind_cols(datos_test |> select(voto))

# Comparar métricas
cat("Árbol de decisión:\n")

print(pred_arbol |> metrics(truth = voto, estimate = .pred_class, .pred_si))

cat("\nRandom Forest:\n")

print(metricas_rf)

# Visualizar el tuning
autoplot(resultados_arbol) +
  theme_minimal() +
  labs(title = "Tuning del árbol de decisión")

# boost_tree(): Gradient Boosting (árboles secuenciales)
# tree_depth y learn_rate a optimizar; min_n fijo
modelo_xgb <- boost_tree(
  trees = 500,
  tree_depth = tune(),
  learn_rate = tune(),
  min_n = 10
) |>
  set_engine("xgboost") |>
  set_mode("classification")

wf_xgb <- workflow() |>
  add_recipe(receta) |>
  add_model(modelo_xgb)

# learn_rate en escala log10: 10^-3 a 10^-1
grilla_xgb <- grid_regular(
  tree_depth(range = c(3, 8)),
  learn_rate(range = c(-3, -1)),
  levels = c(3, 3)
)

resultados_xgb <- tune_grid(
  wf_xgb, resamples = folds, grid = grilla_xgb,
  metrics = metric_set(roc_auc),
  control = control_grid(verbose = FALSE)
)

mejor_xgb <- select_best(resultados_xgb, metric = "roc_auc")
wf_xgb_final <- finalize_workflow(wf_xgb, mejor_xgb)
ajuste_xgb <- fit(wf_xgb_final, data = datos_train)

# Predicciones y comparación final
pred_xgb <- predict(ajuste_xgb, datos_test) |>
  bind_cols(predict(ajuste_xgb, datos_test, type = "prob")) |>
  bind_cols(datos_test |> select(voto))

bind_rows(
  pred_logit |> metrics(truth = voto, estimate = .pred_class, .pred_si) |>
    mutate(modelo = "Regresión logística"),
  pred_rf |> metrics(truth = voto, estimate = .pred_class, .pred_si) |>
    mutate(modelo = "Random Forest"),
  pred_xgb |> metrics(truth = voto, estimate = .pred_class, .pred_si) |>
    mutate(modelo = "XGBoost")
) |>
  select(modelo, .metric, .estimate) |>
  pivot_wider(names_from = .metric, values_from = .estimate) |>
  arrange(desc(roc_auc))