# ============================================================================
# PRONÓSTICO SÍSMICO A 30 DÍAS: CHILE Y CALIFORNIA
# ============================================================================
# Flujo general:
# 1. Configuración y paquetes
# 2. Funciones auxiliares
# 3. Descarga y control de calidad
# 4. Homogeneización y categorización
# 5. Preparación de catálogos y ventanas de 30 días
# 6. Diagnóstico descriptivo de las ventanas
# 7. Modelo base histórico y regresión logística
# 8. Diagnósticos de las regresiones logísticas
# 9. Modelos de conteo Poisson/binomial negativa
# 10. Consolidación, gráficos y visualización de resultados
# ============================================================================
# ============================================================================
# 1. CONFIGURACIÓN Y PAQUETES
# ============================================================================
paquetes_necesarios <- c(
"tidyverse",
"lubridate",
"pROC",
"MASS",
"AER",
"car"
)
paquetes_faltantes <- paquetes_necesarios[
!vapply(paquetes_necesarios, requireNamespace, logical(1), quietly = TRUE)
]
if (length(paquetes_faltantes) > 0) {
install.packages(paquetes_faltantes)
}
suppressPackageStartupMessages({
library(tidyverse)
library(lubridate)
library(pROC)
})
options(dplyr.summarise.inform = FALSE)
# Parámetros generales --------------------------------------------------------
FECHA_INICIO <- as.POSIXct("2000-01-01 00:00:00", tz = "UTC")
FECHA_FIN <- as.POSIXct("2026-01-01 00:00:00", tz = "UTC")
PROPORCION_ENTRENAMIENTO <- 0.80
EPSILON <- 1e-15
DESCARGAR_DESDE_USGS <- TRUE
GUARDAR_BASES_CRUDAS <- TRUE
GENERAR_GRAFICOS <- TRUE
ABRIR_VISTAS <- interactive()
ARCHIVO_CHILE <- "chile_limpia.csv"
ARCHIVO_CALIFORNIA <- "california_limpia.csv"
# ============================================================================
# 2. FUNCIONES AUXILIARES GENERALES
# ============================================================================
imprimir_seccion <- function(titulo) {
cat(
"\n",
paste(rep("=", 76), collapse = ""), "\n",
titulo, "\n",
paste(rep("=", 76), collapse = ""), "\n",
sep = ""
)
}
division_segura <- function(numerador, denominador) {
ifelse(denominador > 0, numerador / denominador, NA_real_)
}
calcular_auc_seguro <- function(observado, probabilidad) {
observado <- as.integer(observado)
probabilidad <- as.numeric(probabilidad)
validos <- is.finite(observado) & is.finite(probabilidad)
observado <- observado[validos]
probabilidad <- probabilidad[validos]
if (dplyr::n_distinct(observado) < 2) {
return(NA_real_)
}
if (dplyr::n_distinct(probabilidad) < 2) {
return(0.5)
}
curva <- pROC::roc(
response = observado,
predictor = probabilidad,
levels = c(0, 1),
direction = "<",
quiet = TRUE
)
as.numeric(pROC::auc(curva))
}
calcular_youden <- function(observado, probabilidad) {
observado <- as.integer(observado)
probabilidad <- as.numeric(probabilidad)
validos <- is.finite(observado) & is.finite(probabilidad)
observado <- observado[validos]
probabilidad <- probabilidad[validos]
auc <- calcular_auc_seguro(observado, probabilidad)
# No existe un punto de corte identificable cuando falta una clase
# o todas las probabilidades son iguales.
if (
dplyr::n_distinct(observado) < 2 ||
dplyr::n_distinct(probabilidad) < 2
) {
return(list(
disponible = FALSE,
roc = NULL,
corte = NA_real_,
sensibilidad = NA_real_,
especificidad = NA_real_,
indice_youden = NA_real_,
auc = auc
))
}
curva <- pROC::roc(
response = observado,
predictor = probabilidad,
levels = c(0, 1),
direction = "<",
quiet = TRUE
)
punto <- as.data.frame(
pROC::coords(
curva,
x = "best",
best.method = "youden",
ret = c("threshold", "sensitivity", "specificity"),
transpose = FALSE
)
)
corte <- as.numeric(punto$threshold[1])
sensibilidad <- as.numeric(punto$sensitivity[1])
especificidad <- as.numeric(punto$specificity[1])
disponible <- all(is.finite(c(corte, sensibilidad, especificidad)))
list(
disponible = disponible,
roc = curva,
corte = if (disponible) corte else NA_real_,
sensibilidad = if (disponible) sensibilidad else NA_real_,
especificidad = if (disponible) especificidad else NA_real_,
indice_youden = if (disponible) sensibilidad + especificidad - 1 else NA_real_,
auc = as.numeric(pROC::auc(curva))
)
}
graficar_roc_youden <- function(resultado_youden, titulo) {
if (!isTRUE(resultado_youden$disponible) || is.null(resultado_youden$roc)) {
message("No se grafica la curva ROC de '", titulo,
"' porque no existe un punto de corte de Youden identificable.")
return(invisible(NULL))
}
curva <- resultado_youden$roc
datos_roc <- data.frame(
uno_menos_especificidad = 1 - curva$specificities,
sensibilidad = curva$sensitivities
) %>%
arrange(uno_menos_especificidad, sensibilidad)
x_youden <- 1 - resultado_youden$especificidad
y_youden <- resultado_youden$sensibilidad
plot(
x = datos_roc$uno_menos_especificidad,
y = datos_roc$sensibilidad,
type = "l",
lwd = 2,
xlim = c(0, 1),
ylim = c(0, 1),
xaxs = "i",
yaxs = "i",
asp = 1,
main = titulo,
xlab = "1 - Especificidad",
ylab = "Sensibilidad"
)
abline(a = 0, b = 1, lty = 2)
segments(
x0 = x_youden, y0 = 0,
x1 = x_youden, y1 = y_youden,
lty = 3
)
segments(
x0 = 0, y0 = y_youden,
x1 = x_youden, y1 = y_youden,
lty = 3
)
points(
x = x_youden,
y = y_youden,
pch = 21,
bg = "white",
lwd = 2,
cex = 1.6
)
text(
x = x_youden,
y = y_youden,
labels = "Punto óptimo",
pos = 4,
offset = 0.7,
cex = 0.8
)
legend(
"bottomright",
legend = c(
paste0("AUC = ", round(resultado_youden$auc, 3)),
paste0("Corte = ", round(resultado_youden$corte, 4)),
paste0("Sensibilidad = ", round(resultado_youden$sensibilidad, 3)),
paste0("Especificidad = ", round(resultado_youden$especificidad, 3)),
paste0("Youden J = ", round(resultado_youden$indice_youden, 3))
),
bty = "n",
cex = 0.8
)
invisible(NULL)
}
calcular_metricas_probabilisticas <- function(
predicciones,
nombre_analisis,
nombre_metodo
) {
auc <- calcular_auc_seguro(
predicciones$observado,
predicciones$probabilidad
)
predicciones %>%
mutate(
probabilidad_segura = pmin(
pmax(probabilidad, EPSILON),
1 - EPSILON
)
) %>%
summarise(
analisis = nombre_analisis,
metodo = nombre_metodo,
observaciones = n(),
ventanas_positivas = sum(observado == 1, na.rm = TRUE),
prevalencia = mean(observado, na.rm = TRUE),
probabilidad_promedio = mean(probabilidad, na.rm = TRUE),
brier = mean((probabilidad - observado)^2, na.rm = TRUE),
log_loss = -mean(
observado * log(probabilidad_segura) +
(1 - observado) * log(1 - probabilidad_segura),
na.rm = TRUE
),
auc = auc
)
}
calcular_matriz_confusion <- function(
predicciones,
nombre_analisis,
nombre_metodo,
punto_corte = NA_real_,
indice_youden = NA_real_
) {
clasificacion_disponible <- any(!is.na(predicciones$predicho))
if (!clasificacion_disponible) {
return(tibble(
analisis = nombre_analisis,
metodo = nombre_metodo,
clasificacion_disponible = FALSE,
punto_corte_youden = punto_corte,
indice_youden = indice_youden,
verdaderos_positivos = NA_integer_,
falsos_negativos = NA_integer_,
verdaderos_negativos = NA_integer_,
falsos_positivos = NA_integer_,
sensibilidad = NA_real_,
especificidad = NA_real_,
precision = NA_real_,
valor_predictivo_negativo = NA_real_,
tasa_falsas_alarmas = NA_real_,
exactitud = NA_real_,
sensibilidad_porcentaje = NA_real_,
especificidad_porcentaje = NA_real_,
precision_porcentaje = NA_real_,
valor_predictivo_negativo_porcentaje = NA_real_,
tasa_falsas_alarmas_porcentaje = NA_real_,
exactitud_porcentaje = NA_real_
))
}
predicciones %>%
summarise(
verdaderos_positivos = sum(predicho == 1 & observado == 1, na.rm = TRUE),
falsos_negativos = sum(predicho == 0 & observado == 1, na.rm = TRUE),
verdaderos_negativos = sum(predicho == 0 & observado == 0, na.rm = TRUE),
falsos_positivos = sum(predicho == 1 & observado == 0, na.rm = TRUE)
) %>%
mutate(
analisis = nombre_analisis,
metodo = nombre_metodo,
clasificacion_disponible = TRUE,
punto_corte_youden = punto_corte,
indice_youden = indice_youden,
sensibilidad = division_segura(
verdaderos_positivos,
verdaderos_positivos + falsos_negativos
),
especificidad = division_segura(
verdaderos_negativos,
verdaderos_negativos + falsos_positivos
),
precision = division_segura(
verdaderos_positivos,
verdaderos_positivos + falsos_positivos
),
valor_predictivo_negativo = division_segura(
verdaderos_negativos,
verdaderos_negativos + falsos_negativos
),
tasa_falsas_alarmas = division_segura(
falsos_positivos,
falsos_positivos + verdaderos_negativos
),
exactitud = division_segura(
verdaderos_positivos + verdaderos_negativos,
verdaderos_positivos + falsos_negativos +
verdaderos_negativos + falsos_positivos
),
sensibilidad_porcentaje = round(100 * sensibilidad, 2),
especificidad_porcentaje = round(100 * especificidad, 2),
precision_porcentaje = round(100 * precision, 2),
valor_predictivo_negativo_porcentaje = round(
100 * valor_predictivo_negativo,
2
),
tasa_falsas_alarmas_porcentaje = round(
100 * tasa_falsas_alarmas,
2
),
exactitud_porcentaje = round(100 * exactitud, 2),
.before = 1
)
}
# ============================================================================
# 3. DESCARGA Y CONTROL DE CALIDAD
# ============================================================================
# Descarga anual para evitar sobrepasar el límite de resultados de la API USGS.
descargar_usgs_anual <- function(
minlat,
maxlat,
minlon,
maxlon,
anio_inicio,
anio_fin,
minmag
) {
base_url <- "https://earthquake.usgs.gov/fdsnws/event/1/query"
lista_datos <- list()
for (anio in anio_inicio:anio_fin) {
starttime <- paste0(anio, "-01-01")
endtime <- paste0(anio + 1, "-01-01")
url <- paste0(
base_url,
"?format=csv",
"&starttime=", starttime,
"&endtime=", endtime,
"&minlatitude=", minlat,
"&maxlatitude=", maxlat,
"&minlongitude=", minlon,
"&maxlongitude=", maxlon,
"&minmagnitude=", minmag,
"&eventtype=earthquake",
"&orderby=time-asc"
)
cat("Descargando año:", anio, "\n")
datos_anio <- tryCatch(
readr::read_csv(url, show_col_types = FALSE),
error = function(e) {
message("Error en año ", anio, ": ", e$message)
NULL
}
)
if (!is.null(datos_anio) && nrow(datos_anio) > 0) {
lista_datos[[as.character(anio)]] <- datos_anio
}
}
dplyr::bind_rows(lista_datos)
}
if (DESCARGAR_DESDE_USGS) {
base_california <- descargar_usgs_anual(
minlat = 32,
maxlat = 42,
minlon = -125,
maxlon = -114,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 3
)
base_chile <- descargar_usgs_anual(
minlat = -46,
maxlat = -17,
minlon = -76,
maxlon = -66,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 3
)
if (GUARDAR_BASES_CRUDAS) {
readr::write_csv(base_chile, ARCHIVO_CHILE)
readr::write_csv(base_california, ARCHIVO_CALIFORNIA)
}
} else {
base_chile <- readr::read_csv(ARCHIVO_CHILE, show_col_types = FALSE)
base_california <- readr::read_csv(ARCHIVO_CALIFORNIA, show_col_types = FALSE)
}
# Alias compatibles con el código original.
Base_CHILE <- base_chile
Base_CA <- base_california
chile <- base_chile
california <- base_california
variables_clave <- c(
"id", "time", "latitude", "longitude", "place",
"depth", "mag", "magType", "type"
)
control_calidad <- tibble(
region = c("Chile", "California"),
observaciones = c(nrow(chile), nrow(california)),
filas_con_faltantes = c(
sum(!complete.cases(chile[, variables_clave])),
sum(!complete.cases(california[, variables_clave]))
),
duplicados_exactos = c(
sum(duplicated(chile)),
sum(duplicated(california))
)
)
imprimir_seccion("CONTROL DE CALIDAD DE LOS CATÁLOGOS")
print(control_calidad, n = Inf)
resumen_tipos_magnitud <- bind_rows(
chile %>%
count(magType, name = "n") %>%
mutate(region = "Chile"),
california %>%
count(magType, name = "n") %>%
mutate(region = "California")
) %>%
group_by(region) %>%
mutate(porcentaje = round(100 * n / sum(n), 2)) %>%
ungroup() %>%
select(region, magType, n, porcentaje)
imprimir_seccion("TIPOS DE MAGNITUD REGISTRADOS")
print(resumen_tipos_magnitud, n = Inf)
# ============================================================================
# 4. HOMOGENEIZACIÓN Y CATEGORIZACIÓN
# ============================================================================
# Chile: conversión a magnitud momento homogeneizada.
chile <- chile %>%
mutate(
magType = tolower(trimws(magType)),
Mw_hom = case_when(
magType %in% c("mw", "mwb", "mwc", "mwr", "mww") ~ mag,
magType == "ml" & depth <= 50 ~ 0.80 * mag + 1.15,
magType == "ml" & depth > 50 ~ 0.94 * mag + 0.30,
magType == "ms" ~ 0.74 * mag + 1.60,
magType == "mb" ~ 1.04 * mag - 0.02,
magType %in% c("m", "md", "mc") ~ mag,
TRUE ~ NA_real_
)
)
# California: se mantiene la magnitud preferida informada por USGS.
clasificar_magnitud <- function(magnitud) {
case_when(
magnitud < 4 ~ "Menor",
magnitud < 5 ~ "Ligera",
magnitud < 6 ~ "Moderada",
magnitud < 7 ~ "Fuerte",
magnitud < 8 ~ "Mayor",
TRUE ~ "Gran Terremoto"
)
}
clasificar_profundidad <- function(profundidad) {
case_when(
profundidad <= 20 ~ "Superficial Superior",
profundidad <= 70 ~ "Superficial Inferior",
profundidad <= 300 ~ "Intermedio",
TRUE ~ "Profundo"
)
}
chile <- chile %>%
mutate(
mag_class = clasificar_magnitud(Mw_hom),
depth_class = clasificar_profundidad(depth)
)
california <- california %>%
mutate(
mag_class = clasificar_magnitud(mag),
depth_class = clasificar_profundidad(depth)
)
# ============================================================================
# 5. CATÁLOGOS, ZONAS Y VENTANAS DE 30 DÍAS
# ============================================================================
asignar_zona_chile <- function(latitude) {
case_when(
latitude >= -25 ~ "1. Norte",
latitude >= -35 ~ "2. Norte chico",
latitude >= -45 ~ "3. Centro-sur",
latitude >= -52 ~ "4. Patagonia norte",
TRUE ~ "5. Austral"
)
}
asignar_zona_california <- function(latitude) {
case_when(
latitude < 35 ~ "1. Sur",
latitude < 38 ~ "2. Centro",
TRUE ~ "3. Norte"
)
}
preparar_catalogo <- function(
datos,
region_nombre,
variable_magnitud,
funcion_zona,
fecha_inicio = FECHA_INICIO,
fecha_fin = FECHA_FIN
) {
datos %>%
mutate(
fecha = lubridate::ymd_hms(time, tz = "UTC", quiet = TRUE),
magnitud = .data[[variable_magnitud]],
region = region_nombre,
zona = as.character(funcion_zona(latitude))
) %>%
filter(
!is.na(fecha),
!is.na(magnitud),
!is.na(zona),
fecha >= fecha_inicio,
fecha < fecha_fin,
tolower(trimws(type)) == "earthquake"
) %>%
arrange(zona, fecha)
}
chile_modelo <- preparar_catalogo(
datos = chile,
region_nombre = "Chile",
variable_magnitud = "Mw_hom",
funcion_zona = asignar_zona_chile
)
california_modelo <- preparar_catalogo(
datos = california,
region_nombre = "California",
variable_magnitud = "mag",
funcion_zona = asignar_zona_california
)
resumen_eventos_zona <- bind_rows(
chile_modelo %>% count(region, zona, name = "eventos"),
california_modelo %>% count(region, zona, name = "eventos")
)
imprimir_seccion("EVENTOS UTILIZABLES POR ZONA")
print(resumen_eventos_zona, n = Inf)
crear_base_30dias_zonas <- function(
datos,
umbral_predictor,
umbral_objetivo,
fecha_inicio = FECHA_INICIO,
fecha_fin = FECHA_FIN
) {
inicios <- seq(
from = fecha_inicio,
to = fecha_fin - lubridate::days(30),
by = "30 days"
)
bloques <- tibble(
bloque = seq_along(inicios),
inicio = inicios,
fin = inicios + lubridate::days(30)
)
fecha_final_util <- max(bloques$fin)
zonas <- sort(unique(datos$zona))
bloques_zonas <- tidyr::crossing(
bloques,
zona = zonas
)
eventos_bloque <- datos %>%
filter(
fecha >= fecha_inicio,
fecha < fecha_final_util
) %>%
mutate(
bloque = findInterval(
as.numeric(fecha),
as.numeric(inicios)
)
) %>%
filter(
bloque >= 1,
bloque <= nrow(bloques)
) %>%
group_by(zona, bloque) %>%
summarise(
x_actual = sum(magnitud >= umbral_predictor, na.rm = TRUE),
y_conteo_actual = sum(magnitud > umbral_objetivo, na.rm = TRUE),
.groups = "drop"
)
bloques_zonas %>%
left_join(eventos_bloque, by = c("zona", "bloque")) %>%
mutate(
x_actual = replace_na(x_actual, 0L),
y_conteo_actual = replace_na(y_conteo_actual, 0L)
) %>%
arrange(zona, bloque) %>%
group_by(zona) %>%
mutate(
y_conteo_futuro = lead(y_conteo_actual),
y_binario_futuro = as.integer(y_conteo_futuro > 0),
x_lag1 = lag(x_actual),
y_binario_actual = as.integer(y_conteo_actual > 0),
y_lag1 = lag(y_binario_actual),
log_x_actual = log1p(x_actual),
log_x_lag1 = log1p(x_lag1)
) %>%
filter(!is.na(y_conteo_futuro)) %>%
ungroup()
}
chile_3_5_zonas <- crear_base_30dias_zonas(
datos = chile_modelo,
umbral_predictor = 3,
umbral_objetivo = 5
)
chile_4_6_zonas <- crear_base_30dias_zonas(
datos = chile_modelo,
umbral_predictor = 4,
umbral_objetivo = 6
)
california_3_5_zonas <- crear_base_30dias_zonas(
datos = california_modelo,
umbral_predictor = 3,
umbral_objetivo = 5
)
california_4_6_zonas <- crear_base_30dias_zonas(
datos = california_modelo,
umbral_predictor = 4,
umbral_objetivo = 6
)
analisis_configuracion <- list(
chile_35 = list(
base = chile_3_5_zonas,
nombre = "Chile: M >= 3 hacia M > 5",
region = "Chile",
predictor = 3,
objetivo = 5
),
chile_46 = list(
base = chile_4_6_zonas,
nombre = "Chile: M >= 4 hacia M > 6",
region = "Chile",
predictor = 4,
objetivo = 6
),
california_35 = list(
base = california_3_5_zonas,
nombre = "California: M >= 3 hacia M > 5",
region = "California",
predictor = 3,
objetivo = 5
),
california_46 = list(
base = california_4_6_zonas,
nombre = "California: M >= 4 hacia M > 6",
region = "California",
predictor = 4,
objetivo = 6
)
)
# ============================================================================
# 6. DIAGNÓSTICO DESCRIPTIVO DE LAS VENTANAS
# ============================================================================
ljung_seguro <- function(x, lag_maximo = 12) {
x <- x[is.finite(x)]
if (length(x) < 10 || length(unique(x)) < 2) {
return(NA_real_)
}
lag_prueba <- min(
lag_maximo,
floor(length(x) / 5),
length(x) - 1
)
if (lag_prueba < 1) {
return(NA_real_)
}
Box.test(
x,
lag = lag_prueba,
type = "Ljung-Box"
)$p.value
}
diagnosticar_por_zona <- function(base, nombre) {
base %>%
group_by(zona) %>%
summarise(
analisis = nombre,
bloques = n(),
ventanas_positivas = sum(y_binario_futuro == 1, na.rm = TRUE),
ventanas_negativas = sum(y_binario_futuro == 0, na.rm = TRUE),
prevalencia = mean(y_binario_futuro, na.rm = TRUE),
total_eventos_futuros = sum(y_conteo_futuro, na.rm = TRUE),
media_conteo = mean(y_conteo_futuro, na.rm = TRUE),
varianza_conteo = var(y_conteo_futuro, na.rm = TRUE),
indice_dispersion = if_else(
media_conteo > 0,
varianza_conteo / media_conteo,
NA_real_
),
p_ljung_x = ljung_seguro(x_actual),
p_ljung_y = ljung_seguro(y_conteo_actual),
.groups = "drop"
) %>%
select(analisis, zona, everything())
}
diagnosticos_zonas <- purrr::map_dfr(
analisis_configuracion,
~ diagnosticar_por_zona(.x$base, .x$nombre)
)
imprimir_seccion("DIAGNÓSTICO DESCRIPTIVO POR ZONA")
print(diagnosticos_zonas, n = Inf, width = Inf)
graficar_acf_zonas <- function(base, variable, titulo) {
zonas <- sort(unique(base$zona))
numero_columnas <- 2
numero_filas <- ceiling(length(zonas) / numero_columnas)
configuracion_anterior <- par(no.readonly = TRUE)
on.exit(par(configuracion_anterior))
par(
mfrow = c(numero_filas, numero_columnas),
mar = c(4, 4, 4, 1)
)
for (zona_actual in zonas) {
valores <- base %>%
filter(zona == zona_actual) %>%
arrange(bloque) %>%
pull(all_of(variable))
valores <- valores[is.finite(valores)]
if (length(valores) < 2 || length(unique(valores)) < 2) {
plot.new()
title(main = paste(titulo, zona_actual, sep = "\n"))
text(0.5, 0.5, "Serie constante o no evaluable")
} else {
acf(
valores,
main = paste(titulo, zona_actual, sep = "\n")
)
}
}
invisible(NULL)
}
# ============================================================================
# 7. MODELO BASE HISTÓRICO Y REGRESIÓN LOGÍSTICA
# ============================================================================
preparar_datos_binarios <- function(base) {
base %>%
mutate(
zona = factor(zona),
log_x_actual = log1p(x_actual),
log_x_lag1 = log1p(x_lag1),
y_binario_actual = as.integer(y_conteo_actual > 0)
) %>%
tidyr::drop_na(
bloque,
inicio,
fin,
zona,
log_x_actual,
log_x_lag1,
y_binario_actual,
y_binario_futuro,
y_conteo_futuro
) %>%
arrange(bloque, zona)
}
dividir_temporalmente <- function(
datos,
nombre_analisis,
proporcion_entrenamiento = PROPORCION_ENTRENAMIENTO
) {
bloques <- sort(unique(datos$bloque))
cantidad_train <- floor(proporcion_entrenamiento * length(bloques))
if (cantidad_train < 1 || cantidad_train >= length(bloques)) {
stop("No fue posible realizar la división temporal 80/20 en: ", nombre_analisis)
}
bloques_train <- bloques[seq_len(cantidad_train)]
bloques_test <- bloques[(cantidad_train + 1):length(bloques)]
train <- datos %>%
filter(bloque %in% bloques_train)
test <- datos %>%
filter(bloque %in% bloques_test) %>%
mutate(zona = factor(zona, levels = levels(train$zona)))
if (dplyr::n_distinct(train$y_binario_futuro) < 2) {
stop("El entrenamiento no contiene ambas clases en: ", nombre_analisis)
}
if (dplyr::n_distinct(test$y_binario_futuro) < 2) {
warning(
"El testeo no contiene ambas clases en: ", nombre_analisis,
". Algunas métricas no podrán calcularse."
)
}
resumen <- tibble(
analisis = nombre_analisis,
bloques_entrenamiento = length(bloques_train),
bloques_testeo = length(bloques_test),
observaciones_entrenamiento = nrow(train),
observaciones_testeo = nrow(test),
primer_bloque_entrenamiento = min(bloques_train),
ultimo_bloque_entrenamiento = max(bloques_train),
primer_bloque_testeo = min(bloques_test),
ultimo_bloque_testeo = max(bloques_test)
)
list(
train = train,
test = test,
bloques_train = bloques_train,
bloques_test = bloques_test,
resumen = resumen
)
}
ajustar_modelo_base <- function(train, test, nombre_analisis) {
probabilidad_global <- (
sum(train$y_binario_futuro == 1, na.rm = TRUE) + 1
) / (
nrow(train) + 2
)
modelo_zona <- train %>%
group_by(zona) %>%
summarise(
bloques = n(),
ventanas_positivas = sum(y_binario_futuro == 1, na.rm = TRUE),
ventanas_negativas = sum(y_binario_futuro == 0, na.rm = TRUE),
probabilidad_historica = (ventanas_positivas + 1) / (bloques + 2),
.groups = "drop"
)
predicciones <- test %>%
left_join(
modelo_zona %>%
select(zona, probabilidad = probabilidad_historica),
by = "zona"
) %>%
mutate(
probabilidad = coalesce(probabilidad, probabilidad_global),
inicio_predictor = inicio,
fin_predictor = fin,
inicio_objetivo = fin,
fin_objetivo = fin + lubridate::days(30),
observado = y_binario_futuro,
conteo_observado = y_conteo_futuro,
analisis = nombre_analisis,
metodo = "Modelo base histórico"
) %>%
select(
bloque,
inicio_predictor,
fin_predictor,
inicio_objetivo,
fin_objetivo,
zona,
observado,
conteo_observado,
probabilidad,
analisis,
metodo
)
metricas <- calcular_metricas_probabilisticas(
predicciones,
nombre_analisis,
"Modelo base histórico"
)
list(
probabilidad_global = probabilidad_global,
modelo_zona = modelo_zona,
predicciones = predicciones,
metricas = metricas
)
}
ajustar_modelo_logistico <- function(train, test, nombre_analisis) {
formula_completa <- y_binario_futuro ~
log_x_actual +
log_x_lag1 +
y_binario_actual +
zona
modelo_completo <- glm(
formula_completa,
data = train,
family = binomial(link = "logit")
)
modelo_nulo <- glm(
y_binario_futuro ~ 1,
data = train,
family = binomial(link = "logit")
)
modelo_step <- step(
modelo_completo,
scope = list(
lower = formula(modelo_nulo),
upper = formula(modelo_completo)
),
direction = "both",
trace = 0
)
if (!isTRUE(modelo_completo$converged)) {
warning("El modelo logístico completo no convergió en: ", nombre_analisis)
}
if (!isTRUE(modelo_step$converged)) {
warning("El modelo logístico stepwise no convergió en: ", nombre_analisis)
}
probabilidad_train <- as.numeric(
predict(modelo_step, newdata = train, type = "response")
)
youden <- calcular_youden(
train$y_binario_futuro,
probabilidad_train
)
probabilidad_test <- as.numeric(
predict(modelo_step, newdata = test, type = "response")
)
if (isTRUE(youden$disponible)) {
predicho <- as.integer(probabilidad_test >= youden$corte)
} else {
predicho <- rep(NA_integer_, length(probabilidad_test))
}
predicciones <- test %>%
mutate(
probabilidad = probabilidad_test,
punto_corte = youden$corte,
predicho = predicho,
observado = y_binario_futuro,
conteo_observado = y_conteo_futuro,
inicio_predictor = inicio,
fin_predictor = fin,
inicio_objetivo = fin,
fin_objetivo = fin + lubridate::days(30),
resultado = case_when(
is.na(predicho) ~ NA_character_,
predicho == 1 & observado == 1 ~ "Verdadero positivo",
predicho == 1 & observado == 0 ~ "Falso positivo",
predicho == 0 & observado == 0 ~ "Verdadero negativo",
predicho == 0 & observado == 1 ~ "Falso negativo",
TRUE ~ NA_character_
),
analisis = nombre_analisis,
metodo = "Regresión logística stepwise"
) %>%
select(
bloque,
inicio_predictor,
fin_predictor,
inicio_objetivo,
fin_objetivo,
zona,
observado,
conteo_observado,
probabilidad,
punto_corte,
predicho,
resultado,
analisis,
metodo
)
metricas <- calcular_metricas_probabilisticas(
predicciones,
nombre_analisis,
"Regresión logística stepwise"
)
matriz_confusion <- calcular_matriz_confusion(
predicciones,
nombre_analisis,
"Regresión logística stepwise",
punto_corte = youden$corte,
indice_youden = youden$indice_youden
)
odds_ratios <- tibble(
variable = names(coef(modelo_step)),
coeficiente = as.numeric(coef(modelo_step)),
odds_ratio = exp(as.numeric(coef(modelo_step)))
)
comparacion_nulo_step <- anova(
modelo_nulo,
modelo_step,
test = "Chisq"
)
list(
modelo_completo = modelo_completo,
modelo_nulo = modelo_nulo,
modelo_step = modelo_step,
odds_ratios = odds_ratios,
comparacion_nulo_step = comparacion_nulo_step,
probabilidad_train = probabilidad_train,
youden = youden,
predicciones = predicciones,
metricas = metricas,
matriz_confusion = matriz_confusion
)
}
ejecutar_analisis_binario <- function(base, nombre_analisis) {
imprimir_seccion(paste("ANÁLISIS BINARIO:", nombre_analisis))
datos <- preparar_datos_binarios(base)
division <- dividir_temporalmente(datos, nombre_analisis)
print(division$resumen, n = Inf)
base_historico <- ajustar_modelo_base(
division$train,
division$test,
nombre_analisis
)
logistico <- ajustar_modelo_logistico(
division$train,
division$test,
nombre_analisis
)
cat("\nProbabilidad histórica global:",
round(base_historico$probabilidad_global, 4), "\n")
print(base_historico$modelo_zona, n = Inf)
cat("\nModelo logístico completo:\n")
print(summary(logistico$modelo_completo))
cat("\nModelo logístico stepwise:\n")
print(summary(logistico$modelo_step))
cat("\nOdds ratios:\n")
print(logistico$odds_ratios, n = Inf)
cat("\nComparación del modelo nulo y stepwise:\n")
print(logistico$comparacion_nulo_step)
if (isTRUE(logistico$youden$disponible)) {
cat(
"\nPunto de corte de Youden:", round(logistico$youden$corte, 4),
"\nSensibilidad de entrenamiento:",
round(logistico$youden$sensibilidad, 4),
"\nEspecificidad de entrenamiento:",
round(logistico$youden$especificidad, 4),
"\nÍndice de Youden:",
round(logistico$youden$indice_youden, 4),
"\nAUC de entrenamiento:",
round(logistico$youden$auc, 4),
"\n"
)
} else {
cat(
"\nEl modelo entrega probabilidades constantes o no dispone de ambas clases;",
" no se calcula un punto de corte de Youden ni métricas de clasificación.\n"
)
}
list(
datos = datos,
train = division$train,
test = division$test,
bloques_train = division$bloques_train,
bloques_test = division$bloques_test,
resumen_division = division$resumen,
base = base_historico,
logistico = logistico
)
}
resultados_binarios <- purrr::map(
analisis_configuracion,
~ ejecutar_analisis_binario(.x$base, .x$nombre)
)
# Alias compatibles con el código original -----------------------------------
resultado_binario_chile_35 <- resultados_binarios$chile_35
resultado_binario_chile_46 <- resultados_binarios$chile_46
resultado_binario_california_35 <- resultados_binarios$california_35
resultado_binario_california_46 <- resultados_binarios$california_46
# Chile 3 -> 5
datos_chile_35 <- resultado_binario_chile_35$datos
train_chile_35 <- resultado_binario_chile_35$train
test_chile_35 <- resultado_binario_chile_35$test
p_global_chile_35 <- resultado_binario_chile_35$base$probabilidad_global
modelo_base_chile_35 <- resultado_binario_chile_35$base$modelo_zona
pred_base_chile_35 <- resultado_binario_chile_35$base$predicciones
metricas_base_chile_35 <- resultado_binario_chile_35$base$metricas
modelo_completo_chile_35 <- resultado_binario_chile_35$logistico$modelo_completo
modelo_nulo_chile_35 <- resultado_binario_chile_35$logistico$modelo_nulo
modelo_step_chile_35 <- resultado_binario_chile_35$logistico$modelo_step
odds_ratios_chile_35 <- resultado_binario_chile_35$logistico$odds_ratios
roc_train_chile_35 <- resultado_binario_chile_35$logistico$youden$roc
corte_youden_chile_35 <- resultado_binario_chile_35$logistico$youden$corte
pred_logistico_chile_35 <- resultado_binario_chile_35$logistico$predicciones
metricas_logistico_chile_35 <- resultado_binario_chile_35$logistico$metricas
matriz_chile_35 <- resultado_binario_chile_35$logistico$matriz_confusion
# Chile 4 -> 6
datos_chile_46 <- resultado_binario_chile_46$datos
train_chile_46 <- resultado_binario_chile_46$train
test_chile_46 <- resultado_binario_chile_46$test
p_global_chile_46 <- resultado_binario_chile_46$base$probabilidad_global
modelo_base_chile_46 <- resultado_binario_chile_46$base$modelo_zona
pred_base_chile_46 <- resultado_binario_chile_46$base$predicciones
metricas_base_chile_46 <- resultado_binario_chile_46$base$metricas
modelo_completo_chile_46 <- resultado_binario_chile_46$logistico$modelo_completo
modelo_nulo_chile_46 <- resultado_binario_chile_46$logistico$modelo_nulo
modelo_step_chile_46 <- resultado_binario_chile_46$logistico$modelo_step
odds_ratios_chile_46 <- resultado_binario_chile_46$logistico$odds_ratios
roc_train_chile_46 <- resultado_binario_chile_46$logistico$youden$roc
corte_youden_chile_46 <- resultado_binario_chile_46$logistico$youden$corte
pred_logistico_chile_46 <- resultado_binario_chile_46$logistico$predicciones
metricas_logistico_chile_46 <- resultado_binario_chile_46$logistico$metricas
matriz_chile_46 <- resultado_binario_chile_46$logistico$matriz_confusion
# California 3 -> 5
datos_california_35 <- resultado_binario_california_35$datos
train_california_35 <- resultado_binario_california_35$train
test_california_35 <- resultado_binario_california_35$test
p_global_california_35 <- resultado_binario_california_35$base$probabilidad_global
modelo_base_california_35 <- resultado_binario_california_35$base$modelo_zona
pred_base_california_35 <- resultado_binario_california_35$base$predicciones
metricas_base_california_35 <- resultado_binario_california_35$base$metricas
modelo_completo_california_35 <- resultado_binario_california_35$logistico$modelo_completo
modelo_nulo_california_35 <- resultado_binario_california_35$logistico$modelo_nulo
modelo_step_california_35 <- resultado_binario_california_35$logistico$modelo_step
odds_ratios_california_35 <- resultado_binario_california_35$logistico$odds_ratios
roc_train_california_35 <- resultado_binario_california_35$logistico$youden$roc
corte_youden_california_35 <- resultado_binario_california_35$logistico$youden$corte
pred_logistico_california_35 <- resultado_binario_california_35$logistico$predicciones
metricas_logistico_california_35 <- resultado_binario_california_35$logistico$metricas
matriz_california_35 <- resultado_binario_california_35$logistico$matriz_confusion
# California 4 -> 6
datos_california_46 <- resultado_binario_california_46$datos
train_california_46 <- resultado_binario_california_46$train
test_california_46 <- resultado_binario_california_46$test
p_global_california_46 <- resultado_binario_california_46$base$probabilidad_global
modelo_base_california_46 <- resultado_binario_california_46$base$modelo_zona
pred_base_california_46 <- resultado_binario_california_46$base$predicciones
metricas_base_california_46 <- resultado_binario_california_46$base$metricas
modelo_completo_california_46 <- resultado_binario_california_46$logistico$modelo_completo
modelo_nulo_california_46 <- resultado_binario_california_46$logistico$modelo_nulo
modelo_step_california_46 <- resultado_binario_california_46$logistico$modelo_step
odds_ratios_california_46 <- resultado_binario_california_46$logistico$odds_ratios
roc_train_california_46 <- resultado_binario_california_46$logistico$youden$roc
corte_youden_california_46 <- resultado_binario_california_46$logistico$youden$corte
pred_logistico_california_46 <- resultado_binario_california_46$logistico$predicciones
metricas_logistico_california_46 <- resultado_binario_california_46$logistico$metricas
matriz_california_46 <- resultado_binario_california_46$logistico$matriz_confusion
# ============================================================================
# 8. DIAGNÓSTICOS DE LAS REGRESIONES LOGÍSTICAS STEPWISE
# ============================================================================
obtener_residuos_pearson <- function(modelo, datos, nombre_modelo) {
residuos <- residuals(modelo, type = "pearson")
if (length(residuos) != nrow(datos)) {
stop(
"La cantidad de residuos no coincide con las observaciones en: ",
nombre_modelo
)
}
datos %>%
mutate(residuo_pearson = as.numeric(residuos)) %>%
arrange(zona, bloque)
}
ljung_box_residuos_seguro <- function(residuos, lag_maximo = 30) {
residuos <- residuos[is.finite(residuos)]
n_residuos <- length(residuos)
if (n_residuos < 10 || length(unique(residuos)) < 2) {
return(tibble(
observaciones = n_residuos,
rezagos = NA_integer_,
estadistico = NA_real_,
grados_libertad = NA_real_,
p_valor = NA_real_,
conclusion = "No fue posible realizar la prueba"
))
}
rezagos <- min(
lag_maximo,
floor(n_residuos / 5),
n_residuos - 1
)
if (rezagos < 1) {
return(tibble(
observaciones = n_residuos,
rezagos = NA_integer_,
estadistico = NA_real_,
grados_libertad = NA_real_,
p_valor = NA_real_,
conclusion = "No fue posible realizar la prueba"
))
}
prueba <- Box.test(
residuos,
lag = rezagos,
type = "Ljung-Box"
)
p_valor <- as.numeric(prueba$p.value)
tibble(
observaciones = n_residuos,
rezagos = rezagos,
estadistico = as.numeric(prueba$statistic),
grados_libertad = as.numeric(prueba$parameter),
p_valor = p_valor,
conclusion = if_else(
p_valor < 0.05,
"Se detecta autocorrelación temporal significativa",
"No se detecta autocorrelación temporal significativa"
)
)
}
modelos_step <- purrr::imap(
resultados_binarios,
~ list(
nombre = analisis_configuracion[[.y]]$nombre,
modelo = .x$logistico$modelo_step,
datos = .x$train
)
)
tabla_ljung <- purrr::map_dfr(
modelos_step,
function(elemento) {
datos_residuos <- obtener_residuos_pearson(
elemento$modelo,
elemento$datos,
elemento$nombre
)
datos_residuos %>%
group_by(zona) %>%
group_modify(
~ ljung_box_residuos_seguro(.x$residuo_pearson, lag_maximo = 30)
) %>%
ungroup() %>%
mutate(analisis = elemento$nombre, .before = 1)
}
)
imprimir_seccion("LJUNG-BOX SOBRE RESIDUOS DE PEARSON")
print(tabla_ljung, n = Inf, width = Inf)
graficar_acf_residuos <- function(
modelo,
datos,
nombre_modelo,
lag_maximo = 30
) {
datos_residuos <- obtener_residuos_pearson(
modelo,
datos,
nombre_modelo
)
zonas <- sort(unique(datos_residuos$zona))
numero_columnas <- 2
numero_filas <- ceiling(length(zonas) / numero_columnas)
configuracion_anterior <- par(no.readonly = TRUE)
on.exit(par(configuracion_anterior))
par(
mfrow = c(numero_filas, numero_columnas),
mar = c(4, 4, 4, 1)
)
for (zona_actual in zonas) {
residuos_zona <- datos_residuos %>%
filter(zona == zona_actual) %>%
arrange(bloque) %>%
pull(residuo_pearson)
residuos_zona <- residuos_zona[is.finite(residuos_zona)]
if (
length(residuos_zona) < 2 ||
length(unique(residuos_zona)) < 2
) {
plot.new()
title(main = paste(nombre_modelo, zona_actual, sep = "\n"))
text(0.5, 0.5, "Serie de residuos no evaluable")
} else {
acf(
residuos_zona,
lag.max = min(lag_maximo, length(residuos_zona) - 1),
main = paste(nombre_modelo, zona_actual, sep = "\n"),
xlab = "Rezago",
ylab = "Autocorrelación"
)
}
}
invisible(NULL)
}
resultados_vif <- purrr::map(
modelos_step,
function(elemento) {
predictores <- attr(terms(elemento$modelo), "term.labels")
if (length(predictores) <= 1) {
return(list(
calculado = FALSE,
motivo = if (length(predictores) == 0) {
"Modelo compuesto solo por el intercepto"
} else {
paste("Modelo con un solo predictor:", predictores)
},
resultado = NULL
))
}
list(
calculado = TRUE,
motivo = NA_character_,
resultado = car::vif(elemento$modelo)
)
}
)
calcular_mcfadden <- function(modelo_step, nombre_modelo) {
modelo_nulo <- update(modelo_step, formula = . ~ 1)
loglik_step <- as.numeric(logLik(modelo_step))
loglik_nulo <- as.numeric(logLik(modelo_nulo))
tibble(
analisis = nombre_modelo,
observaciones = nobs(modelo_step),
logLik_nulo = loglik_nulo,
logLik_step = loglik_step,
devianza_nula = deviance(modelo_nulo),
devianza_step = deviance(modelo_step),
R2_McFadden = 1 - loglik_step / loglik_nulo
)
}
resultados_mcfadden <- purrr::map_dfr(
modelos_step,
~ calcular_mcfadden(.x$modelo, .x$nombre)
)
imprimir_seccion("PSEUDO-R2 DE MCFADDEN")
print(resultados_mcfadden, n = Inf, width = Inf)
imprimir_seccion("VIF / GVIF DE LOS MODELOS STEPWISE")
purrr::iwalk(
resultados_vif,
function(resultado, clave) {
cat("\n", modelos_step[[clave]]$nombre, "\n", sep = "")
if (resultado$calculado) {
print(resultado$resultado)
} else {
cat(resultado$motivo, "\n")
}
}
)
# ============================================================================
# 9. MODELOS DE CONTEO: POISSON O BINOMIAL NEGATIVA
# ============================================================================
ajustar_modelo_conteo <- function(train, test, nombre_analisis) {
train <- train %>%
mutate(
zona = factor(zona),
log_y_actual = log1p(y_conteo_actual)
)
test <- test %>%
mutate(
zona = factor(zona, levels = levels(train$zona)),
log_y_actual = log1p(y_conteo_actual)
)
formula_completa <- y_conteo_futuro ~
log_x_actual +
log_x_lag1 +
log_y_actual +
zona
modelo_poisson_completo <- glm(
formula_completa,
data = train,
family = poisson(link = "log")
)
dispersion_pearson <- sum(
residuals(modelo_poisson_completo, type = "pearson")^2
) / df.residual(modelo_poisson_completo)
prueba_dispersion <- tryCatch(
AER::dispersiontest(
modelo_poisson_completo,
alternative = "greater"
),
error = function(e) e
)
if (inherits(prueba_dispersion, "error")) {
p_dispersion <- NA_real_
estadistico_dispersion <- NA_real_
warning(
"No fue posible ejecutar dispersiontest en ", nombre_analisis,
": ", prueba_dispersion$message
)
} else {
p_dispersion <- as.numeric(prueba_dispersion$p.value)
estadistico_dispersion <- as.numeric(prueba_dispersion$statistic)
}
existe_sobredispersion <- is.finite(p_dispersion) && p_dispersion < 0.05
modelo_nb_completo <- NULL
if (existe_sobredispersion) {
modelo_nb_completo <- MASS::glm.nb(
formula_completa,
data = train
)
modelo_completo_elegido <- modelo_nb_completo
modelo_nulo_elegido <- MASS::glm.nb(
y_conteo_futuro ~ 1,
data = train
)
distribucion_elegida <- "Binomial negativa"
usar_binomial_negativa <- TRUE
} else {
modelo_completo_elegido <- modelo_poisson_completo
modelo_nulo_elegido <- glm(
y_conteo_futuro ~ 1,
data = train,
family = poisson(link = "log")
)
distribucion_elegida <- "Poisson"
usar_binomial_negativa <- FALSE
}
modelo_step <- step(
modelo_completo_elegido,
scope = list(
lower = formula(modelo_nulo_elegido),
upper = formula(modelo_completo_elegido)
),
direction = "both",
trace = 0
)
conteo_train <- as.numeric(
predict(modelo_step, newdata = train, type = "response")
)
if (usar_binomial_negativa) {
theta <- modelo_step$theta
probabilidad_train <- 1 - dnbinom(
0,
size = theta,
mu = conteo_train
)
} else {
theta <- NA_real_
probabilidad_train <- 1 - dpois(
0,
lambda = conteo_train
)
}
youden <- calcular_youden(
train$y_binario_futuro,
probabilidad_train
)
conteo_predicho <- as.numeric(
predict(modelo_step, newdata = test, type = "response")
)
if (usar_binomial_negativa) {
probabilidad_evento <- 1 - dnbinom(
0,
size = theta,
mu = conteo_predicho
)
} else {
probabilidad_evento <- 1 - dpois(
0,
lambda = conteo_predicho
)
}
probabilidad_evento <- pmin(
pmax(probabilidad_evento, EPSILON),
1 - EPSILON
)
if (isTRUE(youden$disponible)) {
predicho <- as.integer(probabilidad_evento >= youden$corte)
} else {
predicho <- rep(NA_integer_, length(probabilidad_evento))
}
nombre_metodo <- paste("Modelo de conteo", distribucion_elegida)
predicciones <- test %>%
transmute(
bloque,
inicio_predictor = inicio,
fin_predictor = fin,
inicio_objetivo = fin,
fin_objetivo = fin + lubridate::days(30),
zona,
conteo_observado = y_conteo_futuro,
conteo_predicho = conteo_predicho,
observado = y_binario_futuro,
ocurrencia_observada = y_binario_futuro,
probabilidad = probabilidad_evento,
punto_corte = youden$corte,
predicho = predicho,
resultado = case_when(
is.na(predicho) ~ NA_character_,
predicho == 1 & observado == 1 ~ "Verdadero positivo",
predicho == 1 & observado == 0 ~ "Falso positivo",
predicho == 0 & observado == 0 ~ "Verdadero negativo",
predicho == 0 & observado == 1 ~ "Falso negativo",
TRUE ~ NA_character_
),
analisis = nombre_analisis,
metodo = nombre_metodo
)
metricas_conteo <- predicciones %>%
summarise(
analisis = first(analisis),
metodo = first(metodo),
observaciones = n(),
media_observada = mean(conteo_observado, na.rm = TRUE),
media_predicha = mean(conteo_predicho, na.rm = TRUE),
sesgo = mean(conteo_predicho - conteo_observado, na.rm = TRUE),
mae = mean(abs(conteo_predicho - conteo_observado), na.rm = TRUE),
rmse = sqrt(
mean((conteo_predicho - conteo_observado)^2, na.rm = TRUE)
)
)
metricas_probabilidad <- calcular_metricas_probabilisticas(
predicciones,
nombre_analisis,
nombre_metodo
)
matriz_confusion <- calcular_matriz_confusion(
predicciones,
nombre_analisis,
nombre_metodo,
punto_corte = youden$corte,
indice_youden = youden$indice_youden
)
diagnostico <- tibble(
analisis = nombre_analisis,
dispersion_pearson = dispersion_pearson,
estadistico_dispersiontest = estadistico_dispersion,
p_dispersiontest = p_dispersion,
sobredispersion_significativa = existe_sobredispersion,
modelo_elegido = distribucion_elegida,
aic_modelo_completo = AIC(modelo_completo_elegido),
aic_modelo_step = AIC(modelo_step),
formula_completa = paste(
deparse(formula(modelo_completo_elegido)),
collapse = ""
),
formula_step = paste(
deparse(formula(modelo_step)),
collapse = ""
),
punto_corte_youden = youden$corte,
indice_youden = youden$indice_youden,
auc_entrenamiento = youden$auc,
theta = theta
)
imprimir_seccion(paste("MODELO DE CONTEO:", nombre_analisis))
cat("\nModelo Poisson completo:\n")
print(summary(modelo_poisson_completo))
cat("\nDispersión de Pearson:", round(dispersion_pearson, 4), "\n")
if (!inherits(prueba_dispersion, "error")) {
print(prueba_dispersion)
}
cat("\nDistribución elegida:", distribucion_elegida, "\n")
cat("\nModelo completo elegido:\n")
print(summary(modelo_completo_elegido))
cat("\nModelo stepwise elegido:\n")
print(summary(modelo_step))
cat("\nFórmula seleccionada:\n")
print(formula(modelo_step))
if (isTRUE(youden$disponible)) {
cat(
"\nCorte de Youden:", round(youden$corte, 4),
"\nSensibilidad de entrenamiento:", round(youden$sensibilidad, 4),
"\nEspecificidad de entrenamiento:", round(youden$especificidad, 4),
"\nÍndice de Youden:", round(youden$indice_youden, 4),
"\nAUC de entrenamiento:", round(youden$auc, 4),
"\n"
)
} else {
cat(
"\nNo se calcula un punto de corte de Youden porque las",
" probabilidades de entrenamiento son constantes o falta una clase.\n"
)
}
list(
modelo = modelo_step,
modelo_completo = modelo_completo_elegido,
modelo_poisson_completo = modelo_poisson_completo,
modelo_binomial_negativa_completo = modelo_nb_completo,
prueba_dispersion = prueba_dispersion,
distribucion = distribucion_elegida,
youden = youden,
roc_train = youden$roc,
corte_youden = youden$corte,
predicciones = predicciones,
matriz_confusion = matriz_confusion,
metricas_conteo = metricas_conteo,
metricas_probabilidad = metricas_probabilidad,
diagnostico = diagnostico
)
}
resultados_conteo <- purrr::imap(
resultados_binarios,
~ ajustar_modelo_conteo(
train = .x$train,
test = .x$test,
nombre_analisis = analisis_configuracion[[.y]]$nombre
)
)
# Alias compatibles con el código original -----------------------------------
resultado_conteo_chile_35 <- resultados_conteo$chile_35
resultado_conteo_chile_46 <- resultados_conteo$chile_46
resultado_conteo_california_35 <- resultados_conteo$california_35
resultado_conteo_california_46 <- resultados_conteo$california_46
pred_conteo_chile_35 <- resultado_conteo_chile_35$predicciones
pred_conteo_chile_46 <- resultado_conteo_chile_46$predicciones
pred_conteo_california_35 <- resultado_conteo_california_35$predicciones
pred_conteo_california_46 <- resultado_conteo_california_46$predicciones
# ============================================================================
# 10. CONSOLIDACIÓN Y PRESENTACIÓN DE RESULTADOS
# ============================================================================
metricas_base_logistico <- purrr::map_dfr(
resultados_binarios,
~ bind_rows(.x$base$metricas, .x$logistico$metricas)
)
metricas_conteo_todos <- purrr::map_dfr(
resultados_conteo,
"metricas_conteo"
)
metricas_prob_conteo_todos <- purrr::map_dfr(
resultados_conteo,
"metricas_probabilidad"
)
diagnostico_conteo_todos <- purrr::map_dfr(
resultados_conteo,
"diagnostico"
)
matrices_logisticas_todas <- purrr::map_dfr(
resultados_binarios,
~ .x$logistico$matriz_confusion
)
matrices_conteo_todas <- purrr::map_dfr(
resultados_conteo,
"matriz_confusion"
)
predicciones_logisticas_todas <- purrr::map_dfr(
resultados_binarios,
~ .x$logistico$predicciones
) %>%
arrange(analisis, bloque, zona)
predicciones_conteo_todas <- purrr::map_dfr(
resultados_conteo,
"predicciones"
) %>%
arrange(analisis, bloque, zona)
comparacion_final_todos <- bind_rows(
metricas_base_logistico,
metricas_prob_conteo_todos
) %>%
arrange(analisis, brier)
# Resultados por región -------------------------------------------------------
metricas_comparacion_chile <- metricas_base_logistico %>%
filter(str_starts(analisis, "Chile"))
metricas_comparacion_california <- metricas_base_logistico %>%
filter(str_starts(analisis, "California"))
metricas_prob_conteo_chile <- metricas_prob_conteo_todos %>%
filter(str_starts(analisis, "Chile"))
metricas_prob_conteo_california <- metricas_prob_conteo_todos %>%
filter(str_starts(analisis, "California"))
metricas_conteo_chile <- metricas_conteo_todos %>%
filter(str_starts(analisis, "Chile"))
metricas_conteo_california <- metricas_conteo_todos %>%
filter(str_starts(analisis, "California"))
diagnostico_conteo_chile <- diagnostico_conteo_todos %>%
filter(str_starts(analisis, "Chile"))
diagnostico_conteo_california <- diagnostico_conteo_todos %>%
filter(str_starts(analisis, "California"))
matriz_confusion_chile <- matrices_logisticas_todas %>%
filter(str_starts(analisis, "Chile"))
matriz_confusion_california <- matrices_logisticas_todas %>%
filter(str_starts(analisis, "California"))
matrices_conteo_chile <- matrices_conteo_todas %>%
filter(str_starts(analisis, "Chile"))
matrices_conteo_california <- matrices_conteo_todas %>%
filter(str_starts(analisis, "California"))
predicciones_logisticas_chile <- predicciones_logisticas_todas %>%
filter(str_starts(analisis, "Chile"))
predicciones_logisticas_california <- predicciones_logisticas_todas %>%
filter(str_starts(analisis, "California"))
comparacion_final_chile <- comparacion_final_todos %>%
filter(str_starts(analisis, "Chile"))
comparacion_final_california <- comparacion_final_todos %>%
filter(str_starts(analisis, "California"))
imprimir_seccion("COMPARACIÓN PROBABILÍSTICA FINAL")
print(comparacion_final_todos, n = Inf, width = Inf)
imprimir_seccion("MÉTRICAS DEL PRONÓSTICO DE CONTEOS")
print(metricas_conteo_todos, n = Inf, width = Inf)
imprimir_seccion("MATRICES DE CONFUSIÓN: REGRESIONES LOGÍSTICAS")
print(matrices_logisticas_todas, n = Inf, width = Inf)
imprimir_seccion("MATRICES DE CONFUSIÓN: MODELOS DE CONTEO")
print(matrices_conteo_todas, n = Inf, width = Inf)
imprimir_seccion("DIAGNÓSTICOS DE LOS MODELOS DE CONTEO")
print(diagnostico_conteo_todos, n = Inf, width = Inf)
# Matrices clásicas observado vs predicho ------------------------------------
imprimir_matriz_clasica <- function(predicciones, titulo) {
imprimir_seccion(titulo)
if (all(is.na(predicciones$predicho))) {
cat("No existe un punto de corte de Youden; no se genera la matriz.\n")
} else {
print(
table(
Observado = predicciones$observado,
Predicho = predicciones$predicho,
useNA = "ifany"
)
)
}
}
purrr::iwalk(
resultados_binarios,
~ imprimir_matriz_clasica(
.x$logistico$predicciones,
paste("REGRESIÓN LOGÍSTICA:", analisis_configuracion[[.y]]$nombre)
)
)
purrr::iwalk(
resultados_conteo,
~ imprimir_matriz_clasica(
.x$predicciones,
paste("MODELO DE CONTEO:", analisis_configuracion[[.y]]$nombre)
)
)
# ============================================================================
# 11. GRÁFICOS
# ============================================================================
if (GENERAR_GRAFICOS) {
# ACF de las series utilizadas en el diagnóstico descriptivo original.
graficar_acf_zonas(
chile_3_5_zonas,
variable = "x_actual",
titulo = paste0(
"Chile: actividad sísmica predictora\n",
"Conteo de sismos M >= 3 por bloque de 30 días"
)
)
graficar_acf_zonas(
chile_4_6_zonas,
variable = "y_conteo_actual",
titulo = paste0(
"Chile: eventos sísmicos objetivo\n",
"Conteo de sismos M > 6 por bloque de 30 días"
)
)
graficar_acf_zonas(
california_3_5_zonas,
variable = "x_actual",
titulo = paste0(
"California: actividad sísmica predictora\n",
"Conteo de sismos M >= 3 por bloque de 30 días"
)
)
graficar_acf_zonas(
california_4_6_zonas,
variable = "y_conteo_actual",
titulo = paste0(
"California: eventos sísmicos objetivo\n",
"Conteo de sismos M > 6 por bloque de 30 días"
)
)
# ACF de residuos de Pearson de los modelos logísticos stepwise.
purrr::walk(
modelos_step,
~ graficar_acf_residuos(
modelo = .x$modelo,
datos = .x$datos,
nombre_modelo = .x$nombre,
lag_maximo = 30
)
)
# Curvas ROC de los modelos logísticos.
purrr::iwalk(
resultados_binarios,
~ graficar_roc_youden(
.x$logistico$youden,
paste0(
"Curva ROC - Entrenamiento\nRegresión logística: ",
analisis_configuracion[[.y]]$nombre
)
)
)
# Curvas ROC de los modelos de conteo.
purrr::iwalk(
resultados_conteo,
~ graficar_roc_youden(
.x$youden,
paste0(
"Curva ROC - Entrenamiento\nModelo de conteo: ",
analisis_configuracion[[.y]]$nombre
)
)
)
}
# ============================================================================
# 12. VISTAS INTERACTIVAS EN RSTUDIO
# ============================================================================
if (ABRIR_VISTAS) {
View(comparacion_final_todos)
View(metricas_conteo_todos)
View(matrices_logisticas_todas)
View(matrices_conteo_todas)
View(diagnostico_conteo_todos)
View(resultados_mcfadden)
View(tabla_ljung)
View(predicciones_logisticas_todas)
View(predicciones_conteo_todas)
}
# ============================================================================
# FIN DEL SCRIPT
# ============================================================================