Se presenta el código fuente en R utilizado para la descarga, procesamiento y análisis de los datos sísmicos correspondientes a la Asesoría III:
############### Script Asesoria III ################
# LIBRERIAS ----
library(readr)
library(dplyr)
library(lubridate)
library(pROC)
library(randomForest)
library(sandwich)
library(lmtest)
library(MASS)
library(AER)
# BASES ----
## OBTENCION DE DATOS ----
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(
read_csv(url, show_col_types = FALSE),
error = function(e) {
message("Error en año ", anio, ": ", e$message)
return(NULL)
}
)
if (!is.null(datos_anio) && nrow(datos_anio) > 0) {
lista_datos[[as.character(anio)]] <- datos_anio
}
}
bind_rows(lista_datos)
}
chileM3 <- descargar_usgs_anual(
minlat = -46,
maxlat = -17,
minlon = -76,
maxlon = -66,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 3
)
chileM4 <- descargar_usgs_anual(
minlat = -46,
maxlat = -17,
minlon = -76,
maxlon = -66,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 4
)
californiaM3 <- descargar_usgs_anual(
minlat = 32,
maxlat = 42,
minlon = -125,
maxlon = -114,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 3
)
californiaM4 <- descargar_usgs_anual(
minlat = 32,
maxlat = 42,
minlon = -125,
maxlon = -114,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 4
)
# Guardado de bases
write.csv(chileM3,"chileM3.csv", row.names = FALSE)
write.csv(chileM4,"chileM4.csv", row.names = FALSE)
write.csv(californiaM3,"californiaM3.csv", row.names = FALSE)
write.csv(californiaM4,"californiaM4.csv", row.names = FALSE)
# Lectura de bases
californiaM3 <- read_csv("californiaM3.csv")
californiaM4 <- read_csv("californiaM4.csv")
chileM3 <- read_csv("chileM3.csv")
chileM4 <- read_csv("chileM4.csv")
californiaM3 <- californiaM3 %>% rename(magnitud = mag)
californiaM4 <- californiaM4 %>% rename(magnitud = mag)
## HOMOGENIZACION CHILE ----
table(chileM3$magType)
table(chileM4$magType)
homg_chile <- function(base)
{
base = base %>%
mutate(
magnitud = 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","mb_lg") ~ mag,
TRUE ~ NA_real_
)
)
}
chileM3=homg_chile(chileM3)
chileM4=homg_chile(chileM4)
# Formato fechas
chileM3$time_fecha <- as.POSIXct(chileM3$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
chileM4$time_fecha <- as.POSIXct(chileM4$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
californiaM3$time_fecha <- as.POSIXct(californiaM3$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
californiaM4$time_fecha <- as.POSIXct(californiaM4$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
## CREACION VENTANAS ----
crear_ventanas <- function(base, mag_min, mag_objetivo, step_dias, ventana_dias = 30) {
base <- base %>% arrange(time_fecha)
fecha_min <- min(base$time_fecha)
fecha_max <- max(base$time_fecha) - days(ventana_dias * 2)
# Secuencia de inicios de ventana X
inicios_X <- seq(fecha_min, fecha_max, by = paste(step_dias, "days"))
resultados <- lapply(inicios_X, function(t_inicio) {
t_fin_X <- t_inicio + days(ventana_dias)
t_fin_Y <- t_fin_X + days(ventana_dias)
# Filtrar sismos por ventana
sismos_X <- base %>% filter(time_fecha >= t_inicio & time_fecha < t_fin_X)
sismos_Y <- base %>% filter(time_fecha >= t_fin_X & time_fecha < t_fin_Y)
# CONSTRUCCIÓN DE PREDICTORES (Xt-1)
# Conteo de sismos menores
conteo_X <- sum(sismos_X$magnitud >= mag_min & sismos_X$magnitud < mag_objetivo, na.rm = TRUE)
prof_media_X <- ifelse(nrow(sismos_X) > 0, mean(sismos_X$depth, na.rm = TRUE), NA)
mag_max_X <- ifelse(nrow(sismos_X) > 0, max(sismos_X$magnitud, na.rm = TRUE), NA)
# Y_{t-1}: Ocurrencia del evento objetivo en la ventana de observación
Y_t_menos_1 <- ifelse(any(sismos_X$magnitud >= mag_objetivo, na.rm = TRUE), 1, 0)
# CONSTRUCCIÓN DE RESPUESTA (Yt)
# Y_t_binario: Ocurrencia del evento objetivo (clasificación)
Y_t_binario <- ifelse(any(sismos_Y$magnitud >= mag_objetivo, na.rm = TRUE), 1, 0)
# Y_t_conteo: Número exacto de sismos mayores (para modelos de conteo)
Y_t_conteo <- sum(sismos_Y$magnitud >= mag_objetivo, na.rm = TRUE)
data.frame(
inicio_obs = t_inicio,
fin_obs = t_fin_X,
fin_pred = t_fin_Y,
X_conteo_menores = conteo_X,
X_prof_media = prof_media_X,
X_mag_max = mag_max_X,
Y_previo = Y_t_menos_1,
Y_objetivo = Y_t_binario,
Y_conteo_objetivo = Y_t_conteo
)
})
return(bind_rows(resultados))
}
## PRUEBA EMPÍRICA: VENTANAS TRASLAPADAS VS NO TRASLAPADAS ---
### Caso M3->M5 ----
chileM3_5_traslapadas <- crear_ventanas(base = chileM3, mag_min = 3.0, mag_objetivo = 5.0, step_dias = 1)
#N ventanas:
nrow(chileM3_5_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(chileM3_5_traslapadas$X_conteo_menores >= 1))
#% ventanas con Y_objetivo = 1
round(100 * mean(chileM3_5_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(chileM3_5_traslapadas$X_conteo_menores)
chileM3_5_no_traslapadas <- crear_ventanas(base = chileM3, mag_min = 3.0, mag_objetivo = 5.0, step_dias = 30)
#N ventanas:
nrow(chileM3_5_no_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(chileM3_5_no_traslapadas$X_conteo_menores >= 1), 2)
#% ventanas con Y_objetivo = 1
round(100 * mean(chileM3_5_no_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(chileM3_5_no_traslapadas$X_conteo_menores)
# Autocorrelación de la variable de ocurrencia del evento objetivo (Y_objetivo)
par(mfrow = c(1, 2))
acf_traslapada <- acf(chileM3_5_traslapadas$Y_objetivo, lag.max = 35, plot = T, main = "ACF: Ventanas Traslapadas\n(Avance Diario)")
print(round(acf_traslapada$acf[1:31], 3))
acf_no_traslapada <- acf(chileM3_5_no_traslapadas$Y_objetivo, lag.max = 10, plot = T, main = "ACF: Ventanas No Traslapadas\n(Bloques 30 días)")
print(round(acf_no_traslapada$acf[1:6], 3))
californiaM3_5_traslapadas <- crear_ventanas(base = californiaM3, mag_min = 3.0, mag_objetivo = 5.0, step_dias = 1)
#N ventanas:
nrow(californiaM3_5_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(californiaM3_5_traslapadas$X_conteo_menores >= 1))
#% ventanas con Y_objetivo = 1
round(100 * mean(californiaM3_5_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(californiaM3_5_traslapadas$X_conteo_menores)
californiaM3_5_no_traslapadas <- crear_ventanas(base = californiaM3, mag_min = 3.0, mag_objetivo = 5.0, step_dias = 30)
#N ventanas:
nrow(californiaM3_5_no_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(californiaM3_5_no_traslapadas$X_conteo_menores >= 1), 2)
#% ventanas con Y_objetivo = 1
round(100 * mean(californiaM3_5_no_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(californiaM3_5_no_traslapadas$X_conteo_menores)
### Caso M4->M6 ----
chileM4_6_traslapadas <- crear_ventanas(base = chileM4, mag_min = 4.0, mag_objetivo = 6.0, step_dias = 1)
#N ventanas:
nrow(chileM4_6_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(chileM4_6_traslapadas$X_conteo_menores >= 1))
#% ventanas con Y_objetivo = 1
round(100 * mean(chileM4_6_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(chileM4_6_traslapadas$X_conteo_menores)
chileM4_6_no_traslapadas <- crear_ventanas(base = chileM4, mag_min = 4.0, mag_objetivo = 6.0, step_dias = 30)
#N ventanas:
nrow(chileM4_6_no_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(chileM4_6_no_traslapadas$X_conteo_menores >= 1), 2)
#% ventanas con Y_objetivo = 1
round(100 * mean(chileM4_6_no_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(chileM4_6_no_traslapadas$X_conteo_menores)
californiaM4_6_traslapadas <- crear_ventanas(base = californiaM4, mag_min = 4.0, mag_objetivo = 6.0, step_dias = 1)
#N ventanas:
nrow(californiaM4_6_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(californiaM4_6_traslapadas$X_conteo_menores >= 1))
#% ventanas con Y_objetivo = 1
round(100 * mean(californiaM4_6_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(californiaM4_6_traslapadas$X_conteo_menores)
californiaM4_6_no_traslapadas <- crear_ventanas(base = californiaM4, mag_min = 4.0, mag_objetivo = 6.0, step_dias = 30)
#N ventanas:
nrow(californiaM4_6_no_traslapadas)
#% ventanas con al menos 1 evento M3-4.99
round(100 * mean(californiaM4_6_no_traslapadas$X_conteo_menores >= 1), 2)
#% ventanas con Y_objetivo = 1
round(100 * mean(californiaM4_6_no_traslapadas$Y_objetivo), 2)
#Resumen de X_conteo_menores (conteo, no solo ocurrencia)
summary(californiaM4_6_no_traslapadas$X_conteo_menores)
# MODELOS CASO CHILE M4->M6 ----
base <- chileM4_6_no_traslapadas
base <- base[!is.na(base$Y_previo), ]
#N total de ventanas utilizables
nrow(base)
#% de ventanas con Y_objetivo = 1
round(100 * mean(base$Y_objetivo), 2)
## PARTICION ENTRENAMIENTO/PRUEBA 80-20 ----
base <- base[order(base$inicio_obs), ]
n_total <- nrow(base)
n_train <- floor(0.8 * n_total)
train <- base[1:n_train, ]
test <- base[(n_train + 1):n_total, ]
cat("Entrenamiento:", nrow(train), "ventanas (", as.character(min(train$inicio_obs)),
"a", as.character(max(train$inicio_obs)), ") -- % positivos:",
round(100 * mean(train$Y_objetivo), 2), "\n")
cat("Prueba: ", nrow(test), "ventanas (", as.character(min(test$inicio_obs)),
"a", as.character(max(test$inicio_obs)), ") -- % positivos:",
round(100 * mean(test$Y_objetivo), 2), "\n")
# Referencia climatológica: tasa base histórica
prob_climatologica <- mean(train$Y_objetivo)
# Brier Score de la referencia climatológica sobre el conjunto de prueba
brier_climatologia <- mean((prob_climatologica - test$Y_objetivo)^2)
brier_climatologia
## METODO 1: Regresión logística simple ----
modelo_M1 <- glm(Y_objetivo ~ X_conteo_menores, data = train,
family = binomial(link = "logit"))
summary(modelo_M1)
pred_M1_test <- predict(modelo_M1, newdata = test, type = "response")
## METODO 2: Modelo autológistico ----
modelo_M2 <- glm(Y_objetivo ~ X_conteo_menores + Y_previo, data = train,
family = binomial(link = "logit"))
summary(modelo_M2)
# Test de razón de verosimilitudes
print(anova(modelo_M1, modelo_M2, test = "LRT"))
pred_M2_test <- predict(modelo_M2, newdata = test, type = "response")
## METODO 3: Random Forest de clasificación ----
train$Y_objetivo_factor <- factor(train$Y_objetivo, levels = c(0, 1), labels = c("No", "Si"))
set.seed(2026) # reproducibilidad del bootstrap interno de Random Forest
modelo_M3 <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max,
data = train,
ntree = 500,
importance = TRUE,
na.action = na.omit
)
modelo_M3
#Importancia de variables (Random Forest)
importance(modelo_M3)
pred_M3_test <- predict(modelo_M3, newdata = test, type = "prob")[, "Si"]
## Métricas de validación retrospectiva: AUC, Brier Score,
tabla_metricas <- data.frame(
metodo = c("Climatologia", "M1: Logistica simple", "M2: Autologistico", "M3: Random Forest"),
AUC = c(NA,
as.numeric(auc(roc(test$Y_objetivo, pred_M1_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M2_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M3_test, quiet = TRUE)))),
Brier_Score = c(brier_climatologia,
mean((pred_M1_test - test$Y_objetivo)^2),
mean((pred_M2_test - test$Y_objetivo)^2),
mean((pred_M3_test - test$Y_objetivo)^2))
)
tabla_metricas$Brier_Skill_Score <- 1 - tabla_metricas$Brier_Score / brier_climatologia
# Comparación final: AUC, Brier Score y Brier Skill Score (conjunto de prueba)
print(tabla_metricas, row.names = FALSE)
## Chequeo de sensibilidad: Random Forest solo con X_conteo_menores,
set.seed(2026)
modelo_M3_simple <- randomForest(
Y_objetivo_factor ~ X_conteo_menores,
data = train,
ntree = 500,
importance = TRUE,
na.action = na.omit
)
# Random Forest (solo X_conteo_menores)
modelo_M3_simple
pred_M3_simple_test <- predict(modelo_M3_simple, newdata = test, type = "prob")[, "Si"]
auc_M3_simple <- as.numeric(auc(roc(test$Y_objetivo, pred_M3_simple_test, quiet = TRUE)))
brier_M3_simple <- mean((pred_M3_simple_test - test$Y_objetivo)^2)
bss_M3_simple <- 1 - brier_M3_simple / brier_climatologia
cat("\nRandom Forest (solo X_conteo_menores) -- AUC:", round(auc_M3_simple, 4),
"| Brier:", round(brier_M3_simple, 4),
"| BSS:", round(bss_M3_simple, 4), "\n")
cat("\nComparar con Random Forest completo (4 covariables): AUC =",
round(as.numeric(auc(roc(test$Y_objetivo, pred_M3_test, quiet = TRUE))), 4), "\n")
# Cambio de ajustes Random Forest
n_positivos_train <- sum(train$Y_objetivo_factor == "Si")
n_negativos_train <- sum(train$Y_objetivo_factor == "No")
cat("\nDistribución de clases en entrenamiento -- No:", n_negativos_train,
"| Si:", n_positivos_train, "\n")
## Ajuste 1: corrección de desbalance vía muestreo estratificado
set.seed(2026)
modelo_M3_balanceado <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max,
data = train,
ntree = 500,
importance = TRUE,
na.action = na.omit,
strata = train$Y_objetivo_factor,
sampsize = c("No" = n_positivos_train, "Si" = n_positivos_train),
nodesize = 10
)
modelo_M3_balanceado
pred_M3_bal_test <- predict(modelo_M3_balanceado, newdata = test, type = "prob")[, "Si"]
auc_M3_bal <- as.numeric(auc(roc(test$Y_objetivo, pred_M3_bal_test, quiet = TRUE)))
brier_M3_bal <- mean((pred_M3_bal_test - test$Y_objetivo)^2)
bss_M3_bal <- 1 - brier_M3_bal / brier_climatologia
cat("\nRandom Forest balanceado -- AUC:", round(auc_M3_bal, 4),
"| Brier:", round(brier_M3_bal, 4),
"| BSS:", round(bss_M3_bal, 4), "\n")
## Comparación resumen: original vs. balanceado
cat("\n=== Comparación: Random Forest original vs. con ajustes ===\n")
tabla_rf_comparacion <- data.frame(
configuracion = c("Original (defaults)", "Balanceado + nodesize=10"),
AUC = c(as.numeric(auc(roc(test$Y_objetivo, pred_M3_test, quiet = TRUE))), auc_M3_bal),
Brier_Score = c(mean((pred_M3_test - test$Y_objetivo)^2), brier_M3_bal),
Brier_Skill_Score = c(1 - mean((pred_M3_test - test$Y_objetivo)^2) / brier_climatologia, bss_M3_bal)
)
print(tabla_rf_comparacion, row.names = FALSE)
## USO DE log_X_conteo EN MODELOS ----
train$log_X_conteo <- log(train$X_conteo_menores + 1) # +1 para evitar log(0)
test$log_X_conteo <- log(test$X_conteo_menores + 1)
modelo_M1_log <- glm(Y_objetivo ~ log_X_conteo, data = train,
family = binomial(link = "logit"))
summary(modelo_M1_log)
pred_M1_log_test <- predict(modelo_M1_log, newdata = test, type = "response")
auc_M1_log <- as.numeric(auc(roc(test$Y_objetivo, pred_M1_log_test, quiet = TRUE)))
brier_M1_log <- mean((pred_M1_log_test - test$Y_objetivo)^2)
bss_M1_log <- 1 - brier_M1_log / brier_climatologia
cat("\nMétodo 1 (log del conteo) -- AUC:", round(auc_M1_log, 4),
"| Brier:", round(brier_M1_log, 4),
"| BSS:", round(bss_M1_log, 4), "\n")
cat("\nComparar con Método 1 (conteo sin transformar): AUC =",
round(as.numeric(auc(roc(test$Y_objetivo, pred_M1_test, quiet = TRUE))), 4),
"| p-valor de X_conteo_menores =",
round(summary(modelo_M1)$coefficients["X_conteo_menores", "Pr(>|z|)"], 4), "\n")
cat("\nAIC -- conteo sin transformar:", round(AIC(modelo_M1), 2),
"| AIC -- log del conteo:", round(AIC(modelo_M1_log), 2), "\n")
# Modelo autologistico
modelo_M2_log <- glm(Y_objetivo ~ log_X_conteo + Y_previo, data = train,
family = binomial(link = "logit"))
print(summary(modelo_M2_log)$coefficients)
print(anova(modelo_M1_log, modelo_M2_log, test = "LRT"))
pred_M2_log_test <- predict(modelo_M2_log, newdata = test, type = "response")
auc_M2_log <- as.numeric(auc(roc(test$Y_objetivo, pred_M2_log_test, quiet = TRUE)))
brier_M2_log <- mean((pred_M2_log_test - test$Y_objetivo)^2)
bss_M2_log <- 1 - brier_M2_log / brier_climatologia
cat("\nMétodo 2 (log del conteo) -- AUC:", round(auc_M2_log, 4),
"| Brier:", round(brier_M2_log, 4),
"| BSS:", round(bss_M2_log, 4), "\n")
cat("\n=== TABLA FINAL: Chile M4.0 -> M6.0 (predictor: log del conteo) ===\n")
tabla_final_chile <- data.frame(
metodo = c("Climatologia",
"M1: Logistica (log conteo)",
"M2: Autologistico (log conteo)",
"M3: Random Forest (balanceado)"),
AUC = c(NA, auc_M1_log, auc_M2_log, auc_M3_bal),
Brier_Score = c(brier_climatologia, brier_M1_log, brier_M2_log, brier_M3_bal),
Brier_Skill_Score = c(0, bss_M1_log, bss_M2_log, bss_M3_bal)
)
print(tabla_final_chile, row.names = FALSE)
# MODELOS CASO CALIFORNIA M3->M5 ----
base <- californiaM3_5_no_traslapadas
base <- base[!is.na(base$Y_previo), ]
#N total de ventanas utilizables
nrow(base)
#% de ventanas con Y_objetivo = 1
round(100 * mean(base$Y_objetivo), 2)
base$log_X_conteo <- log(base$X_conteo_menores + 1)
## PARTICION ENTRENAMIENTO/PRUEBA 80-20 ----
base <- base[order(base$inicio_obs), ]
n_total <- nrow(base)
n_train <- floor(0.8 * n_total)
train <- base[1:n_train, ]
test <- base[(n_train + 1):n_total, ]
cat("Entrenamiento:", nrow(train), "ventanas -- positivos:", sum(train$Y_objetivo),
"(", round(100 * mean(train$Y_objetivo), 2), "%)\n")
cat("Prueba: ", nrow(test), "ventanas -- positivos:", sum(test$Y_objetivo),
"(", round(100 * mean(test$Y_objetivo), 2), "%)\n")
# Referencia climatológica: tasa base histórica
prob_climatologica <- mean(train$Y_objetivo)
# Brier Score de la referencia climatológica sobre el conjunto de prueba
brier_climatologia <- mean((prob_climatologica - test$Y_objetivo)^2)
brier_climatologia
## MÉTODO 1: Regresión logística simple ----
modelo_M1 <- glm(Y_objetivo ~ log_X_conteo, data = train, family = binomial(link = "logit"))
summary(modelo_M1)
pred_M1_test <- predict(modelo_M1, newdata = test, type = "response")
## MÉTODO 2: Modelo autológistico ----
modelo_M2 <- glm(Y_objetivo ~ log_X_conteo + Y_previo, data = train, family = binomial(link = "logit"))
summary(modelo_M2)
print(anova(modelo_M1, modelo_M2, test = "LRT"))
pred_M2_test <- predict(modelo_M2, newdata = test, type = "response")
## MÉTODO 3a: Random Forest ----
train$Y_objetivo_factor <- factor(train$Y_objetivo, levels = c(0, 1), labels = c("No", "Si"))
set.seed(2026)
modelo_M3_default <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max,
data = train, ntree = 500, importance = TRUE, na.action = na.omit
)
print(modelo_M3_default)
pred_M3_default_test <- predict(modelo_M3_default, newdata = test, type = "prob")[, "Si"]
## MÉTODO 3b: Random Forest balanceado + nodesize=10 ----
n_positivos_train <- sum(train$Y_objetivo_factor == "Si")
set.seed(2026)
modelo_M3_balanceado <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max,
data = train, ntree = 500, importance = TRUE, na.action = na.omit,
strata = train$Y_objetivo_factor,
sampsize = c("No" = n_positivos_train, "Si" = n_positivos_train),
nodesize = 10
)
print(modelo_M3_balanceado)
pred_M3_bal_test <- predict(modelo_M3_balanceado, newdata = test, type = "prob")[, "Si"]
## Diagnóstico de observaciones influyentes ----
cooks_M1 <- cooks.distance(modelo_M1)
umbral_cook <- 4 / nrow(train)
influyentes <- which(cooks_M1 > umbral_cook)
cat("\n=== Diagnóstico de influencia (Método distancia de Cook) ===\n")
cat("N° de observaciones influyentes:", length(influyentes), "de", nrow(train), "\n")
if (length(influyentes) > 0) {
cat("\nVentanas más influyentes:\n")
tabla_influyentes <- data.frame(
inicio_obs = train$inicio_obs[influyentes],
X_conteo_menores = train$X_conteo_menores[influyentes],
Y_objetivo = train$Y_objetivo[influyentes],
cooks_dist = round(cooks_M1[influyentes], 4)
)
print(tabla_influyentes[order(-tabla_influyentes$cooks_dist), ], row.names = FALSE)
}
## Métricas de validación retrospectiva: AUC, Brier Score, BSS ----
tabla_metricas <- data.frame(
metodo = c("Climatologia",
"M1: Logistica (log conteo)",
"M2: Autologistico (log conteo)",
"M3a: Random Forest (default)",
"M3b: Random Forest (balanceado)"),
AUC = c(NA,
as.numeric(auc(roc(test$Y_objetivo, pred_M1_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M2_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M3_default_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M3_bal_test, quiet = TRUE)))),
Brier_Score = c(brier_climatologia,
mean((pred_M1_test - test$Y_objetivo)^2),
mean((pred_M2_test - test$Y_objetivo)^2),
mean((pred_M3_default_test - test$Y_objetivo)^2),
mean((pred_M3_bal_test - test$Y_objetivo)^2))
)
tabla_metricas$Brier_Skill_Score <- 1 - tabla_metricas$Brier_Score / brier_climatologia
cat("\n=== TABLA FINAL: California M3.0 -> M5.0 ===\n")
print(tabla_metricas, row.names = FALSE)
## Sensibilidad: reajuste de Método 1 excluyendo las observaciones influyentes ----
summary(modelo_M1)
train_sin_influyentes <- train[-influyentes, ]
modelo_M1_sin_influyentes <- glm(Y_objetivo ~ log_X_conteo, data = train_sin_influyentes,
family = binomial(link = "logit"))
summary(modelo_M1_sin_influyentes)
pred_M1_sin_influyentes_test <- predict(modelo_M1_sin_influyentes, newdata = test, type = "response")
auc_sin_influyentes <- as.numeric(auc(roc(test$Y_objetivo, pred_M1_sin_influyentes_test, quiet = TRUE)))
brier_sin_influyentes <- mean((pred_M1_sin_influyentes_test - test$Y_objetivo)^2)
bss_sin_influyentes <- 1 - brier_sin_influyentes / brier_climatologia
cat("\nMétodo 1 sin influyentes -- AUC:", round(auc_sin_influyentes, 4),
"| Brier:", round(brier_sin_influyentes, 4),
"| BSS:", round(bss_sin_influyentes, 4), "\n")
cat("\nComparar con Método 1 original -- AUC:", round(as.numeric(auc(roc(test$Y_objetivo, pred_M1_test, quiet = TRUE))), 4),
"| p-valor original:", round(summary(modelo_M1)$coefficients["log_X_conteo", "Pr(>|z|)"], 4), "\n")
## Reajuste excluyendo SOLO la observación más extrema
indice_extremo <- which(train$X_conteo_menores == max(train$X_conteo_menores))
train_sin_extremo <- train[-indice_extremo, ]
modelo_M1_sin_extremo <- glm(Y_objetivo ~ log_X_conteo, data = train_sin_extremo,
family = binomial(link = "logit"))
cat("\n=== Método 1, excluyendo SOLO la ventana más extrema (n =", nrow(train_sin_extremo), ") ===\n")
print(summary(modelo_M1_sin_extremo)$coefficients)
# MODELOS CASO CALIFORNIA M4.0->M6.0 ----
base <- californiaM4_6_traslapadas
base <- base[!is.na(base$Y_previo), ]
base <- base[order(base$inicio_obs), ]
base$log_X_conteo <- log(base$X_conteo_menores + 1)
#N total de ventanas (traslapadas, paso diario)
nrow(base)
#% de ventanas con Y_objetivo = 1
round(100 * mean(base$Y_objetivo), 2)
sum(base$Y_objetivo)
cat("\nVentanas sin ningún evento del umbral predictor (NA en X_prof_media/X_mag_max):",
sum(is.na(base$X_prof_media)), "de", nrow(base),
"(", round(100 * mean(is.na(base$X_prof_media)), 2), "%)\n")
base$sin_eventos_menores <- as.integer(is.na(base$X_prof_media))
base$X_prof_media[is.na(base$X_prof_media)] <- 0
base$X_mag_max[is.na(base$X_mag_max)] <- 0
cat("N total de ventanas (traslapadas, paso diario):", nrow(base), "\n")
cat("% de ventanas con Y_objetivo = 1:", round(100 * mean(base$Y_objetivo), 2),
"(", sum(base$Y_objetivo), "positivas )\n")
## Partición por calendario con diferencia de 60 días ----
rango_fechas <- range(base$inicio_obs)
span_dias <- as.numeric(difftime(rango_fechas[2], rango_fechas[1], units = "days"))
fecha_corte <- rango_fechas[1] + days(round(0.8 * span_dias))
purga_dias <- 60
train <- base[base$fin_pred <= fecha_corte, ]
test <- base[base$inicio_obs >= (fecha_corte + days(purga_dias)), ]
cat("\nFecha de corte (80% del rango temporal):", as.character(as.Date(fecha_corte)), "\n")
cat("Zona de purga:", purga_dias, "días a cada lado del corte\n")
cat("Entrenamiento:", nrow(train), "ventanas -- positivos:", sum(train$Y_objetivo),
"(", round(100 * mean(train$Y_objetivo), 2), "%)\n")
cat("Prueba: ", nrow(test), "ventanas -- positivos:", sum(test$Y_objetivo),
"(", round(100 * mean(test$Y_objetivo), 2), "%)\n")
cat("Ventanas descartadas por purga:", nrow(base) - nrow(train) - nrow(test), "\n")
train$bloque_calendario <- floor(as.numeric(difftime(train$inicio_obs, min(train$inicio_obs), units = "days")) / 30)
# Referencia climatológica (solo con entrenamiento)
prob_climatologica <- mean(train$Y_objetivo)
# Brier Score de la referencia climatológica sobre el conjunto de prueba
brier_climatologia <- mean((prob_climatologica - test$Y_objetivo)^2)
brier_climatologia
## MÉTODO 1: Regresión logística simple, con SE robusto por clúster ----
modelo_M1 <- glm(Y_objetivo ~ log_X_conteo, data = train, family = binomial(link = "logit"))
vcov_M1 <- vcovCL(modelo_M1, cluster = ~ bloque_calendario)
coefs_M1_robusto <- coeftest(modelo_M1, vcov = vcov_M1)
summary(modelo_M1)
coefs_M1_robusto
pred_M1_test <- predict(modelo_M1, newdata = test, type = "response")
## MÉTODO 2: Modelo autológistico, con SE robusto y test de Wald robusto para Y_previo ----
modelo_M2 <- glm(Y_objetivo ~ log_X_conteo + Y_previo, data = train, family = binomial(link = "logit"))
vcov_M2 <- vcovCL(modelo_M2, cluster = ~ bloque_calendario)
coefs_M2_robusto <- coeftest(modelo_M2, vcov = vcov_M2)
coefs_M2_robusto
pred_M2_test <- predict(modelo_M2, newdata = test, type = "response")
## MÉTODO 3: Random Forest (default y balanceado) ----
train$Y_objetivo_factor <- factor(train$Y_objetivo, levels = c(0, 1), labels = c("No", "Si"))
set.seed(2026)
modelo_M3_default <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max + sin_eventos_menores,
data = train, ntree = 500, importance = TRUE, na.action = na.omit
)
modelo_M3_default
pred_M3_default_test <- predict(modelo_M3_default, newdata = test, type = "prob")[, "Si"]
n_positivos_train <- sum(train$Y_objetivo_factor == "Si")
set.seed(2026)
modelo_M3_balanceado <- randomForest(
Y_objetivo_factor ~ X_conteo_menores + Y_previo + X_prof_media + X_mag_max + sin_eventos_menores,
data = train, ntree = 500, importance = TRUE, na.action = na.omit,
strata = train$Y_objetivo_factor,
sampsize = c("No" = n_positivos_train, "Si" = n_positivos_train),
nodesize = 10
)
modelo_M3_balanceado
pred_M3_bal_test <- predict(modelo_M3_balanceado, newdata = test, type = "prob")[, "Si"]
## Métricas de validación retrospectiva (conjunto de prueba purgado) ----
tabla_metricas <- data.frame(
metodo = c("Climatologia",
"M1: Logistica (log conteo)",
"M2: Autologistico (log conteo)",
"M3a: Random Forest (default)",
"M3b: Random Forest (balanceado)"),
AUC = c(NA,
as.numeric(auc(roc(test$Y_objetivo, pred_M1_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M2_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M3_default_test, quiet = TRUE))),
as.numeric(auc(roc(test$Y_objetivo, pred_M3_bal_test, quiet = TRUE)))),
Brier_Score = c(brier_climatologia,
mean((pred_M1_test - test$Y_objetivo)^2),
mean((pred_M2_test - test$Y_objetivo)^2),
mean((pred_M3_default_test - test$Y_objetivo)^2),
mean((pred_M3_bal_test - test$Y_objetivo)^2))
)
tabla_metricas$Brier_Skill_Score <- 1 - tabla_metricas$Brier_Score / brier_climatologia
cat("\n=== TABLA FINAL: California M4.0 -> M6.0 (ventanas traslapadas + purga) ===\n")
print(tabla_metricas, row.names = FALSE)
# Diagnostico adicional
summary(pred_M3_default_test)
cat("N° de valores únicos:", length(unique(pred_M3_default_test)), "de", length(pred_M3_default_test), "\n")
cat("% de predicciones iguales a 0:", round(100 * mean(pred_M3_default_test == 0), 2), "\n")
print(table(round(pred_M3_default_test, 2)))
# MODELOS DE CONTEO ----
## CHILE M4->M6 ----
base_conteo <- chileM4_6_no_traslapadas[!is.na(chileM4_6_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
modelo_pois1 <- glm(Y_conteo_objetivo ~ log_X_conteo,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois1))
modelo_pois2 <- glm(Y_conteo_objetivo ~ X_conteo_menores,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois2))
# Prueba formal de sobredispersión de Cameron & Trivedi
prueba_disp <- dispersiontest(modelo_pois1, trafo = 1)
print(prueba_disp)
# Lógica de decisión algorítmica:
if(prueba_disp$p.value < 0.05) {
cat("\nExiste sobredispersión significativa. Se procede a ajustar Binomial Negativa.\n")
modelo_base <- glm.nb(Y_conteo_objetivo ~ log_X_conteo, data = base_conteo)
} else {
cat("\nNo se detecta sobredispersión. El modelo Poisson es adecuado.\n")
modelo_base <- modelo_pois1
}
# Ajustar modelo completo incluyendo profundidad y magnitud máxima
if(prueba_disp$p.value < 0.05) {
modelo_completo <- glm.nb(Y_conteo_objetivo ~ log_X_conteo + X_prof_media + X_mag_max,
data = base_conteo)
} else {
modelo_completo <- glm(Y_conteo_objetivo ~ log_X_conteo + X_prof_media + X_mag_max,
data = base_conteo, family = poisson(link = "log"))
}
summary(modelo_completo)
# Comparación mediante test de Razón de Verosimilitudes (LRT)
print(anova(modelo_base, modelo_completo, test = "LRT"))
cat("\nResumen del Modelo Final Seleccionado:\n")
if(anova(modelo_base, modelo_completo, test = "LRT")$`Pr(>Chi)`[2] < 0.05) {
print(summary(modelo_completo))
} else {
print(summary(modelo_base))
}
## CHILE M3->M5 ----
base_conteo <- chileM3_5_no_traslapadas[!is.na(chileM3_5_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
modelo_pois1 <- glm(Y_conteo_objetivo ~ log_X_conteo,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois1))
modelo_pois2 <- glm(Y_conteo_objetivo ~ X_conteo_menores,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois2))
# Prueba formal de sobredispersión de Cameron & Trivedi
prueba_disp <- dispersiontest(modelo_pois2, trafo = 1)
print(prueba_disp)
# Lógica de decisión algorítmica:
if(prueba_disp$p.value < 0.05) {
cat("\nExiste sobredispersión significativa. Se procede a ajustar Binomial Negativa.\n")
modelo_base <- glm.nb(Y_conteo_objetivo ~ X_conteo_menores, data = base_conteo)
} else {
cat("\nNo se detecta sobredispersión. El modelo Poisson es adecuado.\n")
modelo_base <- modelo_pois2
}
# Ajustar modelo completo incluyendo profundidad y magnitud máxima
if(prueba_disp$p.value < 0.05) {
modelo_completo <- glm.nb(Y_conteo_objetivo ~ X_conteo_menores + X_prof_media + X_mag_max,
data = base_conteo)
} else {
modelo_completo <- glm(Y_conteo_objetivo ~ X_conteo_menores + X_prof_media + X_mag_max,
data = base_conteo, family = poisson(link = "log"))
}
summary(modelo_completo)
modelo_completo2 <- glm(Y_conteo_objetivo ~ X_conteo_menores + X_mag_max,
data = base_conteo, family = poisson(link = "log"))
summary(modelo_completo2)
# Comparación mediante test de Razón de Verosimilitudes (LRT)
print(anova(modelo_base, modelo_completo2, test = "LRT"))
cat("\nResumen del Modelo Final Seleccionado:\n")
# Si el p-valor del LRT < 0.05, el modelo completo aporta información significativa
if(anova(modelo_base, modelo_completo2, test = "LRT")$`Pr(>Chi)`[2] < 0.05) {
print(summary(modelo_completo2))
} else {
print(summary(modelo_base))
}
## CALIFORNIA M4->M6 ----
base_conteo <- californiaM4_6_no_traslapadas[!is.na(californiaM4_6_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
modelo_pois1 <- glm(Y_conteo_objetivo ~ log_X_conteo,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois1))
modelo_pois2 <- glm(Y_conteo_objetivo ~ X_conteo_menores,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois2))
# Prueba formal de sobredispersión de Cameron & Trivedi
prueba_disp <- dispersiontest(modelo_pois1, trafo = 1)
print(prueba_disp)
# Lógica de decisión algorítmica:
if(prueba_disp$p.value < 0.05) {
cat("\nExiste sobredispersión significativa. Se procede a ajustar Binomial Negativa.\n")
modelo_base <- glm.nb(Y_conteo_objetivo ~ log_X_conteo, data = base_conteo)
} else {
cat("\nNo se detecta sobredispersión. El modelo Poisson es adecuado.\n")
modelo_base <- modelo_pois1
}
# Ajustar modelo completo incluyendo profundidad y magnitud máxima
if(prueba_disp$p.value < 0.05) {
modelo_completo <- glm.nb(Y_conteo_objetivo ~ log_X_conteo + X_prof_media + X_mag_max,
data = base_conteo)
} else {
modelo_completo <- glm(Y_conteo_objetivo ~ log_X_conteo + X_prof_media + X_mag_max,
data = base_conteo, family = poisson(link = "log"))
}
summary(modelo_completo)
# Comparación mediante test de Razón de Verosimilitudes (LRT)
print(anova(modelo_base, modelo_completo, test = "LRT"))
cat("\nResumen del Modelo Final Seleccionado:\n")
# Si el p-valor del LRT < 0.05, el modelo completo aporta información significativa
if(anova(modelo_base, modelo_completo, test = "LRT")$`Pr(>Chi)`[2] < 0.05) {
print(summary(modelo_completo))
} else {
print(summary(modelo_base))
}
## CALIFORNIA M3->M5 ----
base_conteo <- californiaM3_5_no_traslapadas[!is.na(californiaM3_5_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
modelo_pois1 <- glm(Y_conteo_objetivo ~ log_X_conteo,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois1))
modelo_pois2 <- glm(Y_conteo_objetivo ~ X_conteo_menores,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois2))
modelo_pois3 <- glm(Y_conteo_objetivo ~ X_mag_max,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois3))
modelo_pois4 <- glm(Y_conteo_objetivo ~ X_prof_media,
data = base_conteo,
family = poisson(link = "log"))
print(summary(modelo_pois4))
# Prueba formal de sobredispersión de Cameron & Trivedi
prueba_disp <- dispersiontest(modelo_pois4, trafo = 1)
print(prueba_disp)
# Lógica de decisión algorítmica:
if(prueba_disp$p.value < 0.05) {
cat("\nExiste sobredispersión significativa. Se procede a ajustar Binomial Negativa.\n")
modelo_base <- glm.nb(Y_conteo_objetivo ~ X_prof_media, data = base_conteo)
} else {
cat("\nNo se detecta sobredispersión. El modelo Poisson es adecuado.\n")
modelo_base <- modelo_pois4
}
summary(modelo_base)
# Ajustar modelo completo incluyendo profundidad y magnitud máxima
if(prueba_disp$p.value < 0.05) {
modelo_completo <- glm.nb(Y_conteo_objetivo ~ X_conteo_menores + X_prof_media + X_mag_max,
data = base_conteo)
} else {
modelo_completo <- glm(Y_conteo_objetivo ~ X_conteo_menores + X_prof_media + X_mag_max,
data = base_conteo, family = poisson(link = "log"))
}
summary(modelo_completo)
# Comparación mediante test de Razón de Verosimilitudes (LRT)
print(anova(modelo_base, modelo_completo, test = "LRT"))
cat("\nResumen del Modelo Final Seleccionado:\n")
# Si el p-valor del LRT < 0.05, el modelo completo aporta información significativa
if(anova(modelo_base, modelo_completo, test = "LRT")$`Pr(Chi)`[2] < 0.05) {
print(summary(modelo_completo))
} else {
print(summary(modelo_base))
}
summary(modelo_base)
# COMPLEMENTO A LOS MODELOS DE CONTEO ----
probabilidad_alerta <- function(modelo, newdata) {
lambda_hat <- predict(modelo, newdata = newdata, type = "response")
if (inherits(modelo, "negbin")) {
theta_hat <- modelo$theta
prob <- 1 - (theta_hat / (theta_hat + lambda_hat))^theta_hat
} else {
prob <- 1 - exp(-lambda_hat)
}
return(prob)
}
## Validación retrospectiva de cada modelo final, partición cronológica 80/20 ----
### --- Chile M4->M6 ---
base_conteo <- chileM4_6_no_traslapadas[!is.na(chileM4_6_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo <- base_conteo[order(base_conteo$inicio_obs), ]
n_train <- floor(0.8 * nrow(base_conteo))
train_c <- base_conteo[1:n_train, ]
test_c <- base_conteo[(n_train + 1):nrow(base_conteo), ]
modelo_final_chileM4M6 <- glm(Y_conteo_objetivo ~ log_X_conteo, data = train_c,
family = poisson(link = "log"))
prob_clima_c <- mean(train_c$Y_conteo_objetivo >= 1)
brier_clima_c <- mean((prob_clima_c - (test_c$Y_conteo_objetivo >= 1))^2)
prob_modelo_c <- probabilidad_alerta(modelo_final_chileM4M6, test_c)
brier_modelo_c <- mean((prob_modelo_c - (test_c$Y_conteo_objetivo >= 1))^2)
cat("\n=== Validación retrospectiva -- Chile M4->M6 (modelo de conteo) ===\n")
cat("Brier Score climatología:", round(brier_clima_c, 4), "\n")
cat("Brier Score modelo: ", round(brier_modelo_c, 4), "\n")
cat("Brier Skill Score: ", round(1 - brier_modelo_c / brier_clima_c, 4), "\n")
### --- Chile M3->M5 ---
base_conteo <- chileM3_5_no_traslapadas[!is.na(chileM3_5_no_traslapadas$Y_previo), ]
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
base_conteo <- base_conteo[order(base_conteo$inicio_obs), ]
n_train <- floor(0.8 * nrow(base_conteo))
train_c <- base_conteo[1:n_train, ]
test_c <- base_conteo[(n_train + 1):nrow(base_conteo), ]
modelo_final_chileM3M5 <- glm(Y_conteo_objetivo ~ X_conteo_menores + X_mag_max, data = train_c,
family = poisson(link = "log"))
prob_clima_c <- mean(train_c$Y_conteo_objetivo >= 1)
brier_clima_c <- mean((prob_clima_c - (test_c$Y_conteo_objetivo >= 1))^2)
prob_modelo_c <- probabilidad_alerta(modelo_final_chileM3M5, test_c)
brier_modelo_c <- mean((prob_modelo_c - (test_c$Y_conteo_objetivo >= 1))^2)
cat("\n=== Validación retrospectiva -- Chile M3->M5 (modelo de conteo) ===\n")
cat("Brier Score climatología:", round(brier_clima_c, 4), "\n")
cat("Brier Score modelo: ", round(brier_modelo_c, 4), "\n")
cat("Brier Skill Score: ", round(1 - brier_modelo_c / brier_clima_c, 4), "\n")
## --- California M4->M6 ---
base_conteo <- californiaM4_6_no_traslapadas[!is.na(californiaM4_6_no_traslapadas$Y_previo), ]
base_conteo$log_X_conteo <- log(base_conteo$X_conteo_menores + 1)
base_conteo <- base_conteo[order(base_conteo$inicio_obs), ]
n_train <- floor(0.8 * nrow(base_conteo))
train_c <- base_conteo[1:n_train, ]
test_c <- base_conteo[(n_train + 1):nrow(base_conteo), ]
modelo_final_califM4M6 <- glm(Y_conteo_objetivo ~ log_X_conteo, data = train_c,
family = poisson(link = "log"))
prob_clima_c <- mean(train_c$Y_conteo_objetivo >= 1)
brier_clima_c <- mean((prob_clima_c - (test_c$Y_conteo_objetivo >= 1))^2)
prob_modelo_c <- probabilidad_alerta(modelo_final_califM4M6, test_c)
brier_modelo_c <- mean((prob_modelo_c - (test_c$Y_conteo_objetivo >= 1))^2)
cat("\n=== Validación retrospectiva -- California M4->M6 (modelo de conteo) ===\n")
cat("Brier Score climatología:", round(brier_clima_c, 4), "\n")
cat("Brier Score modelo: ", round(brier_modelo_c, 4), "\n")
cat("Brier Skill Score: ", round(1 - brier_modelo_c / brier_clima_c, 4), "\n")
## --- California M3->M5 ---
base_conteo <- californiaM3_5_no_traslapadas[!is.na(californiaM3_5_no_traslapadas$Y_previo), ]
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo <- base_conteo[order(base_conteo$inicio_obs), ]
n_train <- floor(0.8 * nrow(base_conteo))
train_c <- base_conteo[1:n_train, ]
test_c <- base_conteo[(n_train + 1):nrow(base_conteo), ]
modelo_final_califM3M5 <- glm.nb(Y_conteo_objetivo ~ X_prof_media, data = train_c)
prob_clima_c <- mean(train_c$Y_conteo_objetivo >= 1)
brier_clima_c <- mean((prob_clima_c - (test_c$Y_conteo_objetivo >= 1))^2)
prob_modelo_c <- probabilidad_alerta(modelo_final_califM3M5, test_c)
brier_modelo_c <- mean((prob_modelo_c - (test_c$Y_conteo_objetivo >= 1))^2)
cat("\n=== Validación retrospectiva -- California M3->M5 (modelo de conteo) ===\n")
cat("Brier Score climatología:", round(brier_clima_c, 4), "\n")
cat("Brier Score modelo: ", round(brier_modelo_c, 4), "\n")
cat("Brier Skill Score: ", round(1 - brier_modelo_c / brier_clima_c, 4), "\n")
## Diagnóstico de precisión: valores sin redondear y conteo de casos negativos en el conjunto de prueba de Chile M3->M5 ----
base_conteo <- chileM3_5_no_traslapadas[!is.na(chileM3_5_no_traslapadas$Y_previo), ]
base_conteo$X_prof_media[is.na(base_conteo$X_prof_media)] <- 0
base_conteo$X_mag_max[is.na(base_conteo$X_mag_max)] <- 0
base_conteo <- base_conteo[order(base_conteo$inicio_obs), ]
n_train <- floor(0.8 * nrow(base_conteo))
train_c <- base_conteo[1:n_train, ]
test_c <- base_conteo[(n_train + 1):nrow(base_conteo), ]
modelo_final_chileM3M5 <- glm(Y_conteo_objetivo ~ X_conteo_menores + X_mag_max, data = train_c,
family = poisson(link = "log"))
prob_clima_c <- mean(train_c$Y_conteo_objetivo >= 1)
prob_modelo_c <- probabilidad_alerta(modelo_final_chileM3M5, test_c)
cat("\n=== Diagnóstico de precisión: Chile M3->M5 (conteo) ===\n")
cat("N ventanas de prueba:", nrow(test_c), "\n")
cat("N ventanas de prueba con Y_conteo_objetivo == 0 (casos 'tranquilos'):",
sum(test_c$Y_conteo_objetivo == 0), "\n")
cat("\nBrier Score climatología (sin redondear):",
format(mean((prob_clima_c - (test_c$Y_conteo_objetivo >= 1))^2), scientific = FALSE, digits = 8), "\n")
cat("Brier Score modelo (sin redondear): ",
format(mean((prob_modelo_c - (test_c$Y_conteo_objetivo >= 1))^2), scientific = FALSE, digits = 8), "\n")
indices_tranquilos <- which(test_c$Y_conteo_objetivo == 0)
if (length(indices_tranquilos) > 0) {
cat("\nProbabilidad pronosticada por el modelo en las ventanas tranquilas:\n")
print(round(prob_modelo_c[indices_tranquilos], 4))
cat("(Comparar con la probabilidad climatológica:", round(prob_clima_c, 4), ")\n")
}