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 II:
############### Script Asesoria II ################
# Librerias ----
library(readr)
library(dplyr)
library(stringr)
library(corrplot)
library(ggplot2)
library(vcd)
library(MASS)
library(lmtest)
library(car)
library(stats)
library(lubridate)
library(tidyr)
library(pscl)
# Obtención 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)
}
Base_CA <- descargar_usgs_anual(
minlat = 32,
maxlat = 42,
minlon = -125,
maxlon = -114,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 5
)
Base_CHILE <- descargar_usgs_anual(
minlat = -46,
maxlat = -17,
minlon = -76,
maxlon = -66,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 5
)
# GUARDADO DE LA BASE CRUDA
getwd()
write.csv(chile,"chile.csv", row.names = FALSE)
write.csv(california, "california.csv", row.names = FALSE)
# verificación de datos faltantes
variables_clave <- c(
"id", "time", "latitude", "longitude", "place",
"depth", "mag", "magType", "type"
)
eliminadas_california <- sum(!complete.cases(california[, variables_clave]))
eliminadas_chile <- sum(!complete.cases(chile[, variables_clave]))
eliminadas_california #Sin obs con datos faltantes
eliminadas_chile #Sin obs con datos faltantes
# duplicados
duplicados_exactos_california <- california[duplicated(california), ]
nrow(duplicados_exactos_california)
duplicados_exactos_chile <- chile[duplicated(chile), ]
nrow(duplicados_exactos_chile)
# Transformacion Escalas magnitud ----
table(chile$magType)
table(california$magType)
california %>%
count(magType) %>%
mutate(porcentaje = round(100 * n / sum(n), 2))
chile %>%
count(magType) %>%
mutate(porcentaje = round(100 * n / sum(n), 2))
## Homogenizacion Chile ----
chile <- chile %>%
mutate(
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_
)
)
## Homogenizacion California -- NO REALIZADA -- Filtracion de datos solo a Mw
california <- california[
tolower(california$magType) %in% c("mw","mwr"),
]
table(california$magType)
# Categorización Magnitud y Profundidad ----
## Magnitud ----
chile <- chile %>% mutate(mag_class=ifelse(Mw_hom < 4.0,"Menor",
ifelse(Mw_hom < 5.0, "Ligera",
ifelse(Mw_hom < 6.0,"Moderada",
ifelse(Mw_hom < 7.0, "Fuerte",
ifelse(Mw_hom < 8.0, "Mayor","Gran Terremoto"))))))
california <- california %>% mutate(mag_class=ifelse(mag < 4.0,"Menor",
ifelse(mag < 5.0, "Ligera",
ifelse(mag < 6.0,"Moderada",
ifelse(mag < 7.0, "Fuerte",
ifelse(mag < 8.0, "Mayor","Gran Terremoto"))))))
table(chile$mag_class)
table(california$mag_class)
## Profundidad ----
chile <- chile %>% mutate(depth_class=ifelse(depth <= 20, "Superficial Superior",
ifelse(depth <= 70, "Superficial Inferior",
ifelse(depth <= 300, "Intermedio", "Profundo"))))
california <- california %>% mutate(depth_class=ifelse(depth <= 20, "Superficial Superior",
ifelse(depth <= 70, "Superficial Inferior",
ifelse(depth <= 300, "Intermedio", "Profundo"))))
table(chile$depth_class)
table(california$depth_class)
# Tests de medidas de asociación ----
chile$magnitud <- chile$Mw_hom
california$magnitud <- california$mag
## Definición de zonas y variables a testear ----
zonas <- list(Chile = chile, California = california)
variables <- c("magnitud", "depth", "latitude", "longitude")
pares <- combn(variables, 2) # 6 combinaciones de a pares
n_pares <- ncol(pares)
resultados <- data.frame(
zona = character(), var1 = character(), var2 = character(),
metodo = character(), n = integer(),
estadistico = numeric(), valor_asociacion = numeric(),
p_valor = numeric(), stringsAsFactors = FALSE
)
for (nombre_zona in names(zonas)) {
df <- zonas[[nombre_zona]]
for (j in 1:n_pares) {
v1 <- pares[1, j]
v2 <- pares[2, j]
x <- df[[v1]]
y <- df[[v2]]
n <- length(x)
# --- Spearman ---
# exact = FALSE: para aproximación asintótica con corrección por empates
test_spearman <- cor.test(x, y, method = "spearman", exact = FALSE)
resultados <- rbind(resultados, data.frame(
zona = nombre_zona, var1 = v1, var2 = v2, metodo = "Spearman", n = n,
estadistico = as.numeric(test_spearman$statistic),
valor_asociacion = as.numeric(test_spearman$estimate),
p_valor = test_spearman$p.value
))
# --- Kendall (tau-b) ---
# exact = FALSE: misma razón
test_kendall <- cor.test(x, y, method = "kendall", exact = FALSE)
resultados <- rbind(resultados, data.frame(
zona = nombre_zona, var1 = v1, var2 = v2, metodo = "Kendall", n = n,
estadistico = as.numeric(test_kendall$statistic),
valor_asociacion = as.numeric(test_kendall$estimate),
p_valor = test_kendall$p.value
))
}
}
resultados <- resultados[order(resultados$zona, resultados$var1, resultados$var2), ]
print(resultados, digits = 4, row.names = FALSE)
## Corrección por multiplicidad ----
resultados$p_fdr <- p.adjust(resultados$p_valor, method = "BH")
alfa <- 0.05
resultados$sig_cruda <- resultados$p_valor < alfa
resultados$sig_fdr <- resultados$p_fdr < alfa
print(resultados[, c("zona", "var1", "var2", "metodo", "valor_asociacion",
"p_valor", "sig_cruda",
"p_fdr", "sig_fdr")], digits = 3, row.names = FALSE)
### Matrices de Correlacion ----
resultados$etiqueta <- paste0(
sprintf("%.2f", resultados$valor_asociacion),
"\n(p=", format.pval(resultados$p_valor, digits = 2, eps = 0.001),
ifelse(resultados$sig_fdr, "*", ""), ")"
)
filas_orden <- c("magnitud", "depth", "latitude")
columnas_orden <- c("depth", "latitude", "longitude")
resultados$var1 <- factor(resultados$var1, levels = filas_orden)
resultados$var2 <- factor(resultados$var2, levels = columnas_orden)
resultados$metodo <- factor(resultados$metodo, levels = c("Spearman", "Kendall"))
for (nombre_zona in c("Chile", "California")) {
datos_zona <- resultados[resultados$zona == nombre_zona, ]
grafico <- ggplot(datos_zona, aes(x = var2, y = var1, fill = valor_asociacion)) +
geom_tile(color = "white") +
geom_text(aes(label = etiqueta, color = abs(valor_asociacion) > 0.5),
size = 3, show.legend = FALSE) +
scale_color_manual(values = c("black", "white")) +
facet_wrap(~ metodo, nrow = 1) +
scale_fill_gradient2(low = "firebrick", mid = "white", high = "steelblue",
midpoint = 0, limits = c(-1, 1), name = "Valor de\nasociación") +
labs(title = paste("Matrices de correlación —", nombre_zona),
subtitle = "Valores de coeficientes junto con sus p-valores sin corrección — * significativo tras FDR",
x = NULL, y = NULL) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(face = "bold"))
print(grafico)
ggsave(paste0("matriz_correlacion_", tolower(nombre_zona), ".png"),
plot = grafico, width = 7, height = 4, dpi = 300)
}
# Chi cuadrado y V de Cramer ----
# Agrupacion Magnitud Chile (3 niveles)
chile <- chile %>%
mutate(mag_class2 = factor(
case_when(
Mw_hom < 6.0 ~ "Moderada",
Mw_hom < 7.0 ~ "Fuerte",
TRUE ~ "Mayor o superior"
),
levels = c("Moderada", "Fuerte", "Mayor o superior")
))
# Magnitud California (2 niveles)
california <- california %>%
mutate(mag_class2 = factor(
case_when(
mag < 6.0 ~ "Moderada",
TRUE ~ "Fuerte o superior"
),
levels = c("Moderada", "Fuerte o superior")
))
# Profundidad a Factor Chile (3 niveles)
chile <- chile %>%
mutate(depth_factor = factor(
case_when(
depth <= 20 ~ "Superficial Superior",
depth <= 70 ~ "Superficial Inferior",
TRUE ~ "Intermedio"
),
levels = c("Superficial Superior", "Superficial Inferior", "Intermedio")
))
# Profundidad a Factor California (2 niveles)
california <- california %>%
mutate(depth_factor = factor(
case_when(
depth <= 20 ~ "Superficial Superior",
TRUE ~ "Superficial Inferior"
),
levels = c("Superficial Superior", "Superficial Inferior")
))
## Caso: Profundidad vs Magnitud - Chile ----
tabla1<-table("Profundidad" = chile$depth_factor, "Magnitud"= chile$mag_class2)
chisq.test(tabla1, simulate.p.value = TRUE)
# Test exacto de Fisher Generalizado
fisher_chile <- tryCatch({
fisher.test(tabla1)
}, error = function(e) {
# Si la memoria colapsa, R entra aquí automáticamente y usa Monte Carlo
cat("\n[Aviso]: N muy grande. Calculando p-valor de Fisher mediante Simulación de Monte Carlo (10.000 réplicas)...\n")
fisher.test(tabla_chile, simulate.p.value = TRUE, B = 10000)
})
print(fisher_chile)
# V de Cramer
assocstats(tabla1)$cramer
## Caso: Profundidad vs Magnitud - California ----
tabla2<-table("Profundidad" = california$depth_factor, "Magnitud"= california$mag_class2)
chisq.test(tabla2, simulate.p.value = TRUE)
fisher.test(tabla2)
assocstats(tabla2)$cramer
## Caso: Zona vs Magnitud ----
chile <- chile %>%
mutate(mag_class_conjunta = factor(
ifelse(Mw_hom < 6.0, "Moderada", "Fuerte o superior"),
levels = c("Moderada", "Fuerte o superior")
))
california <- california %>%
mutate(mag_class_conjunta = factor(
ifelse(mag < 6.0, "Moderada", "Fuerte o superior"),
levels = c("Moderada", "Fuerte o superior")
))
conjunto <- bind_rows(
chile %>% mutate(zona = "Chile") %>%
dplyr::select(zona, mag_class_conjunta, depth_factor),
california %>% mutate(zona = "California") %>%
dplyr::select(zona, mag_class_conjunta, depth_factor)
) %>%
mutate(zona = factor(zona, levels = c("Chile", "California")))
tabla3<-table("Zona" = conjunto$zona, "Magnitud" = conjunto$mag_class_conjunta)
chisq.test(tabla3, simulate.p.value = TRUE)
assocstats(tabla3)$cramer
### Caso: Zona VS Profundidad ----
tabla4<-table("Zona" = conjunto$zona, "Profundidad" = conjunto$depth_factor)
chisq.test(tabla4, simulate.p.value = TRUE)
assocstats(tabla4)$cramer
# Pruebas No Parametricas ----
## Profundidad ----
# Mann-Whitney (Wilcoxon Rank Sum Test)
wilcox.test(chile$depth, california$depth, exact = FALSE)
# Kolmogorov-Smirnov
ks.test(chile$depth, california$depth)
## Magnitud ----
# Mann-Whitney (Wilcoxon Rank Sum Test)
wilcox.test(chile$Mw_hom, california$mag, exact = FALSE)
# Kolmogorov-Smirnov
ks.test(chile$Mw_hom, california$mag)
# Tests de dependencia temporal (réplicas) ----
## Preparación: fecha y tiempo ----
chile$time_fecha <- as.POSIXct(chile$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
california$time_fecha <- as.POSIXct(california$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
chile$anio_mes <- format(chile$time_fecha, "%Y-%m")
california$anio_mes <- format(california$time_fecha, "%Y-%m")
zonas_tiempo <- list(Chile = chile, California = california)
## Indice de dispersion de Fisher sobre conteos mensuales ----
for (nombre_zona in names(zonas_tiempo)) {
df <- zonas_tiempo[[nombre_zona]]
# Se completan los meses sin eventos (conteo = 0), para no subestimar
# la varianza real de la serie
fecha_min <- as.Date(format(min(df$time_fecha, na.rm = TRUE), "%Y-%m-01"))
fecha_max <- as.Date(format(max(df$time_fecha, na.rm = TRUE), "%Y-%m-01"))
todos_los_meses <- format(seq(fecha_min, fecha_max, by = "month"), "%Y-%m")
conteo_mensual <- table(df$anio_mes)
conteo_completo <- setNames(rep(0, length(todos_los_meses)), todos_los_meses)
conteo_completo[names(conteo_mensual)] <- conteo_mensual
T_meses <- length(conteo_completo)
media_n <- mean(conteo_completo)
var_n <- var(conteo_completo)
D <- (T_meses - 1) * var_n / media_n
p_D <- pchisq(D, df = T_meses - 1, lower.tail = FALSE)
cat("\n===", nombre_zona, "- Índice de dispersión (conteos mensuales, M>=5.0) ===\n")
cat("Media:", round(media_n, 3), "| Varianza:", round(var_n, 3),
"| Razón Var/Media:", round(var_n / media_n, 3), "\n")
cat("D =", round(D, 2), "| df =", T_meses - 1,
"| p-valor =", format.pval(p_D, digits = 3, eps = 0.001), "\n")
}
## Tiempos entre eventos: Cv y KS vs Exponencial ----
for (nombre_zona in names(zonas_tiempo)) {
df <- zonas_tiempo[[nombre_zona]]
df <- df[order(df$time_fecha), ]
tiempos_espera <- as.numeric(diff(df$time_fecha), units = "days")
tiempos_espera <- tiempos_espera[tiempos_espera > 0]
media_tiempo <- mean(tiempos_espera)
sd_tiempo <- sd(tiempos_espera)
CV <- sd_tiempo / media_tiempo
# Test de Kolmogorov-Smirnov contra una Exponencial con la misma media
# observada (tasa = 1/media)
test_ks <- ks.test(tiempos_espera, "pexp", rate = 1 / media_tiempo)
cat("\n===", nombre_zona, "- Tiempos entre eventos (M>=5.0, n =", nrow(df), ") ===\n")
cat("Media:", round(media_tiempo, 3), "días | Desv. Est.:", round(sd_tiempo, 3), "\n")
cat("Coeficiente de variación (CV):", round(CV, 3),
"(CV=1 esperado bajo independencia/Poisson homogéneo)\n")
cat("KS D =", round(test_ks$statistic, 4),
"| p-valor =", format.pval(test_ks$p.value, digits = 3, eps = 0.001), "\n")
}
# Declustering: algoritmo de Gardner-Knopoff (1974) ----
## Chile ----
if (!"magnitud" %in% names(chile)) chile$magnitud <- chile$Mw_hom
if (!"time_fecha" %in% names(chile)) chile$time_fecha <- as.POSIXct(chile$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
orden_mag <- order(chile$magnitud, decreasing = TRUE)
df_gk <- chile[orden_mag, ]
n <- nrow(df_gk)
lat <- df_gk$latitude
lon <- df_gk$longitude
tiempo <- df_gk$time_fecha
mag <- df_gk$magnitud
es_dependiente <- rep(FALSE, n)
cluster_id <- df_gk$id
for (i in 1:n) {
if (es_dependiente[i]) next
M_i <- mag[i]
L_i <- 10^(0.1238 * M_i + 0.983)
T_i <- if (M_i >= 6.5) 10^(0.032 * M_i + 2.7389) else 10^(0.5409 * M_i - 0.547)
# Distancia de Haversine (km) del evento i a todos los demás eventos
phi1 <- lat[i] * pi / 180
phi2 <- lat * pi / 180
dphi <- (lat - lat[i]) * pi / 180
dlambda <- (lon - lon[i]) * pi / 180
a <- sin(dphi / 2)^2 + cos(phi1) * cos(phi2) * sin(dlambda / 2)^2
distancia <- 6371 * 2 * atan2(sqrt(a), sqrt(1 - a))
# Diferencia temporal en días (positiva = ocurre después del evento i)
dif_tiempo <- as.numeric(difftime(tiempo, tiempo[i], units = "days"))
en_ventana <- (distancia <= L_i) &
(dif_tiempo > 0) & (dif_tiempo <= T_i) &
(mag <= M_i) &
(!es_dependiente)
en_ventana[i] <- FALSE # excluir el propio evento
if (any(en_ventana)) {
es_dependiente[en_ventana] <- TRUE
cluster_id[en_ventana] <- cluster_id[i]
}
}
df_gk$es_dependiente <- es_dependiente
df_gk$cluster_id <- cluster_id
chile <- df_gk[order(df_gk$time_fecha), ]
#Eventos totales:
nrow(chile)
#Eventos independientes:
sum(!chile$es_dependiente)
#Réplicas removidas:
sum(chile$es_dependiente)
#N° de secuencias con al menos una réplica:
length(unique(chile$cluster_id[chile$es_dependiente]))
evento_max_chile <- chile[which.max(chile$magnitud), ]
#Evento de mayor magnitud en Chile:
evento_max_chile$id
format(evento_max_chile$time_fecha, "%Y-%m-%d")
evento_max_chile$magnitud
#Réplicas atribuidas a este evento:
sum(chile$cluster_id == evento_max_chile$id & chile$es_dependiente)
## California ----
if (!"magnitud" %in% names(california)) california$magnitud <- california$mag
if (!"time_fecha" %in% names(california)) california$time_fecha <- as.POSIXct(california$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
orden_mag <- order(california$magnitud, decreasing = TRUE)
df_gk <- california[orden_mag, ]
n <- nrow(df_gk)
lat <- df_gk$latitude
lon <- df_gk$longitude
tiempo <- df_gk$time_fecha
mag <- df_gk$magnitud
es_dependiente <- rep(FALSE, n)
cluster_id <- df_gk$id
for (i in 1:n) {
if (es_dependiente[i]) next
M_i <- mag[i]
L_i <- 10^(0.1238 * M_i + 0.983)
T_i <- if (M_i >= 6.5) 10^(0.032 * M_i + 2.7389) else 10^(0.5409 * M_i - 0.547)
phi1 <- lat[i] * pi / 180
phi2 <- lat * pi / 180
dphi <- (lat - lat[i]) * pi / 180
dlambda <- (lon - lon[i]) * pi / 180
a <- sin(dphi / 2)^2 + cos(phi1) * cos(phi2) * sin(dlambda / 2)^2
distancia <- 6371 * 2 * atan2(sqrt(a), sqrt(1 - a))
dif_tiempo <- as.numeric(difftime(tiempo, tiempo[i], units = "days"))
en_ventana <- (distancia <= L_i) &
(dif_tiempo > 0) & (dif_tiempo <= T_i) &
(mag <= M_i) &
(!es_dependiente)
en_ventana[i] <- FALSE
if (any(en_ventana)) {
es_dependiente[en_ventana] <- TRUE
cluster_id[en_ventana] <- cluster_id[i]
}
}
df_gk$es_dependiente <- es_dependiente
df_gk$cluster_id <- cluster_id
california <- df_gk[order(df_gk$time_fecha), ]
#Eventos totales:
nrow(california)
#Eventos independientes:
sum(!california$es_dependiente)
#Réplicas removidas:
sum(california$es_dependiente)
#N° de secuencias con al menos una réplica:
length(unique(california$cluster_id[california$es_dependiente]))
evento_max_california <- california[which.max(california$magnitud), ]
#Evento de mayor magnitud en California:
evento_max_california$id
format(evento_max_california$time_fecha, "%Y-%m-%d")
evento_max_california$magnitud
#Réplicas atribuidas a este evento:
sum(california$cluster_id == evento_max_california$id & california$es_dependiente)
# Tests de dependencia temporal (catalogo desclusterizado) ----
# Bases con solo sismos principales
chile_declust <- chile[!chile$es_dependiente, ]
calif_declust <- california[!california$es_dependiente, ]
# Fecha y tiempo
chile_declust$time_fecha <- as.POSIXct(chile_declust$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
chile_declust$anio_mes <- format(chile_declust$time_fecha, "%Y-%m")
calif_declust$time_fecha <- as.POSIXct(calif_declust$time, format = "%Y-%m-%dT%H:%M:%OS", tz = "UTC")
calif_declust$anio_mes <- format(calif_declust$time_fecha, "%Y-%m")
chile_declust <- chile_declust[order(chile_declust$time_fecha), ]
calif_declust <- calif_declust[order(calif_declust$time_fecha), ]
zonas_declust <- list(Chile_Declust = chile_declust, California_Declust = calif_declust)
## Test Indice de dispersion sobre conteos mensuales (Solo Principales) ----
for (nombre_zona in names(zonas_declust)) {
df <- zonas_declust[[nombre_zona]]
fecha_min <- as.Date(format(min(df$time_fecha, na.rm = TRUE), "%Y-%m-01"))
fecha_max <- as.Date(format(max(df$time_fecha, na.rm = TRUE), "%Y-%m-01"))
todos_los_meses <- format(seq(fecha_min, fecha_max, by = "month"), "%Y-%m")
conteo_mensual <- table(df$anio_mes)
conteo_completo <- setNames(rep(0, length(todos_los_meses)), todos_los_meses)
conteo_completo[names(conteo_mensual)] <- conteo_mensual
T_meses <- length(conteo_completo)
media_n <- mean(conteo_completo)
var_n <- var(conteo_completo)
D <- (T_meses - 1) * var_n / media_n
p_D <- pchisq(D, df = T_meses - 1, lower.tail = FALSE)
cat("\n===", nombre_zona, "- Dispersión Mensual (Eventos Independientes) ===\n")
cat("N =", nrow(df), "sismos principales\n")
cat("Media:", round(media_n, 3), "| Varianza:", round(var_n, 3),
"| Razón Var/Media:", round(var_n / media_n, 3), "\n")
cat("D =", round(D, 2), "| df =", T_meses - 1,
"| p-valor =", format.pval(p_D, digits = 3, eps = 0.001), "\n")
}
## Test Tiempos de Espera y CV ----
for (nombre_zona in names(zonas_declust)) {
df <- zonas_declust[[nombre_zona]]
tiempos_espera <- as.numeric(diff(df$time_fecha), units = "days")
tiempos_espera <- tiempos_espera[tiempos_espera > 0]
media_tiempo <- mean(tiempos_espera)
sd_tiempo <- sd(tiempos_espera)
CV <- sd_tiempo / media_tiempo
test_ks <- ks.test(tiempos_espera, "pexp", rate = 1 / media_tiempo)
cat("\n===", nombre_zona, "- Tiempos entre eventos (Eventos Independientes) ===\n")
cat("Coeficiente de variación (CV):", round(CV, 3), "(Se espera ~1.0 bajo Poisson)\n")
cat("KS D =", round(test_ks$statistic, 4),
"| p-valor =", format.pval(test_ks$p.value, digits = 3, eps = 0.001), "\n")
}
# Modelo de conteo: Productividad de réplicas por sismo principal ----
## Numero de replicas por sismo principal Chile ----
mainshocks_chile <- chile[!chile$es_dependiente, ]
conteo_replicas_chile <- table(chile$cluster_id[chile$es_dependiente])
mainshocks_chile$N_replicas <- 0
posicion <- match(mainshocks_chile$id, names(conteo_replicas_chile))
mainshocks_chile$N_replicas[!is.na(posicion)] <-
as.numeric(conteo_replicas_chile[posicion[!is.na(posicion)]])
# Chile: distribución de numero de réplicas por sismo principal
summary(mainshocks_chile$N_replicas)
# Sismos principales sin réplicas
sum(mainshocks_chile$N_replicas == 0)
nrow(mainshocks_chile)
### Modelo Poisson ---
modelo_prod_chile_pois <- glm(N_replicas ~ magnitud, data = mainshocks_chile,
family = poisson(link = "log"))
cat("\n=== Chile: modelo Poisson (N_replicas ~ Magnitud) ===\n")
summary(modelo_prod_chile_pois)
# Sobre-dispersión (estadístico de Pearson / df)
dispersion_chile <- sum(residuals(modelo_prod_chile_pois, type = "pearson")^2) /
modelo_prod_chile_pois$df.residual
# Estadístico de dispersión (Pearson/df) (>>1 indica sobre-dispersión; se prefiere Binomial Negativa)
round(dispersion_chile, 3)
### Modelo Binomial Negativa
---
modelo_prod_chile_nb <- glm.nb(N_replicas ~ magnitud, data = mainshocks_chile)
summary(modelo_prod_chile_nb)
# Comparación de ajuste (AIC): Poisson - Binomial Negativa
round(AIC(modelo_prod_chile_pois), 1)
round(AIC(modelo_prod_chile_nb), 1)
# Interpretación: razón de tasas de incidencia (IRR = exp(beta))
coefs_prod_chile <- summary(modelo_prod_chile_nb)$coefficients
IRR_chile <- exp(coefs_prod_chile[, "Estimate"])
ic_inf_chile <- exp(coefs_prod_chile[, "Estimate"] - qnorm(0.975) * coefs_prod_chile[, "Std. Error"])
ic_sup_chile <- exp(coefs_prod_chile[, "Estimate"] + qnorm(0.975) * coefs_prod_chile[, "Std. Error"])
tabla_irr_chile <- data.frame(
coeficiente = rownames(coefs_prod_chile),
IRR = round(IRR_chile, 3),
IC_95_inf = round(ic_inf_chile, 3),
IC_95_sup = round(ic_sup_chile, 3),
p_valor = format.pval(coefs_prod_chile[, "Pr(>|z|)"], digits = 3, eps = 0.001)
)
#Razón de tasas (IRR) e IC 95%, modelo Binomial Negativa
print(tabla_irr_chile, row.names = FALSE)
## Numero de replicas por sismo principal California ----
mainshocks_california <- california[!california$es_dependiente, ]
conteo_replicas_california <- table(california$cluster_id[california$es_dependiente])
mainshocks_california$N_replicas <- 0
posicion <- match(mainshocks_california$id, names(conteo_replicas_california))
mainshocks_california$N_replicas[!is.na(posicion)] <- as.numeric(conteo_replicas_california[posicion[!is.na(posicion)]])
# California: distribución de numero de réplicas por sismo principal
summary(mainshocks_california$N_replicas)
# Sismos principales sin réplicas
sum(mainshocks_california$N_replicas == 0)
nrow(mainshocks_california)
### Modelo Poisson ---
modelo_prod_california_pois <- glm(N_replicas ~ magnitud, data = mainshocks_california,
family = poisson(link = "log"))
print(summary(modelo_prod_california_pois)$coefficients)
# Chequeo de sobre-dispersión
dispersion_california <- sum(residuals(modelo_prod_california_pois, type = "pearson")^2) /
modelo_prod_california_pois$df.residual
# Estadístico de dispersión (Pearson/df) (>>1 indica sobre-dispersión; se prefiere Binomial Negativa)
round(dispersion_california, 3)
### Modelo Binomial Negativa ---
modelo_prod_california_nb <- tryCatch(
glm.nb(N_replicas ~ magnitud, data = mainshocks_california),
warning = function(w) {
# Advertencia al ajustar Binomial Negativa en California
conditionMessage(w)
suppressWarnings(glm.nb(N_replicas ~ magnitud, data = mainshocks_california))
}
)
summary(modelo_prod_california_nb)
# Theta (parámetro de dispersión)
round(modelo_prod_california_nb$theta, 3)
# Comparación de ajuste (AIC): Poisson vs Binomial Negativa
round(AIC(modelo_prod_california_pois), 1)
round(AIC(modelo_prod_california_nb), 1)
coefs_prod_california <- summary(modelo_prod_california_nb)$coefficients
IRR_california <- exp(coefs_prod_california[, "Estimate"])
ic_inf_california <- exp(coefs_prod_california[, "Estimate"] - qnorm(0.975) * coefs_prod_california[, "Std. Error"])
ic_sup_california <- exp(coefs_prod_california[, "Estimate"] + qnorm(0.975) * coefs_prod_california[, "Std. Error"])
tabla_irr_california <- data.frame(
coeficiente = rownames(coefs_prod_california),
IRR = round(IRR_california, 3),
IC_95_inf = round(ic_inf_california, 3),
IC_95_sup = round(ic_sup_california, 3),
p_valor = format.pval(coefs_prod_california[, "Pr(>|z|)"], digits = 3, eps = 0.001)
)
# Razón de tasas (IRR) e IC 95%, modelo Binomial Negativa
print(tabla_irr_california, row.names = FALSE)
# Resumen: productividad de réplicas por unidad de magnitud
resumen_productividad <- data.frame(
zona = c("Chile", "California"),
IRR_mag = c(tabla_irr_chile$IRR[tabla_irr_chile$coeficiente == "magnitud"],
tabla_irr_california$IRR[tabla_irr_california$coeficiente == "magnitud"]),
IC_95_inf = c(tabla_irr_chile$IC_95_inf[tabla_irr_chile$coeficiente == "magnitud"],
tabla_irr_california$IC_95_inf[tabla_irr_california$coeficiente == "magnitud"]),
IC_95_sup = c(tabla_irr_chile$IC_95_sup[tabla_irr_chile$coeficiente == "magnitud"],
tabla_irr_california$IC_95_sup[tabla_irr_california$coeficiente == "magnitud"]),
p_valor = c(tabla_irr_chile$p_valor[tabla_irr_chile$coeficiente == "magnitud"],
tabla_irr_california$p_valor[tabla_irr_california$coeficiente == "magnitud"])
)
print(resumen_productividad, row.names = FALSE)
## Chile: diagnóstico de influencia y sensibilidad ----
### Distancia de Cook sobre el modelo Binomial Negativa ---
cooks_chile <- cooks.distance(modelo_prod_chile_nb)
umbral_cook <- 4 / nrow(mainshocks_chile) # regla práctica habitual: 4/n
influyentes <- which(cooks_chile > umbral_cook)
#Umbral (4/n)
round(umbral_cook, 5)
# N° de observaciones influyentes
length(influyentes)
nrow(mainshocks_chile)
tabla_influyentes <- data.frame(
id = mainshocks_chile$id[influyentes],
fecha = format(mainshocks_chile$time_fecha[influyentes], "%Y-%m-%d"),
magnitud = mainshocks_chile$magnitud[influyentes],
N_replicas = mainshocks_chile$N_replicas[influyentes],
cooks_dist = round(cooks_chile[influyentes], 4)
)
tabla_influyentes <- tabla_influyentes[order(-tabla_influyentes$cooks_dist), ]
# Mainshocks más influyentes (ordenados por distancia de Cook)
print(head(tabla_influyentes, 10), row.names = FALSE)
### Sensibilidad 1: reajuste excluyendo las observaciones influyentes (según distancia de Cook) ----
chile_sin_influyentes <- mainshocks_chile[-influyentes, ]
modelo_chile_nb_sin_influyentes <- glm.nb(N_replicas ~ magnitud,
data = chile_sin_influyentes)
coefs_sin_influyentes <- summary(modelo_chile_nb_sin_influyentes)$coefficients
IRR_sin_influyentes <- exp(coefs_sin_influyentes["magnitud", "Estimate"])
ic_inf_sin_influyentes <- exp(coefs_sin_influyentes["magnitud", "Estimate"] -
qnorm(0.975) * coefs_sin_influyentes["magnitud", "Std. Error"])
ic_sup_sin_influyentes <- exp(coefs_sin_influyentes["magnitud", "Estimate"] +
qnorm(0.975) * coefs_sin_influyentes["magnitud", "Std. Error"])
# Chile: IRR de magnitud, excluyendo observaciones influyentes (Cook)
cat("N =", nrow(chile_sin_influyentes), "(se excluyeron", length(influyentes), "observaciones)\n")
cat("IRR =", round(IRR_sin_influyentes, 3),
"| IC 95% = [", round(ic_inf_sin_influyentes, 3), ",", round(ic_sup_sin_influyentes, 3), "]\n")
### Sensibilidad 2: reajuste excluyendo eventos M>=8.0 ----
chile_sin_M8 <- mainshocks_chile[mainshocks_chile$magnitud < 8.0, ]
# Eventos M>=8.0 en Chile (excluidos en esta sensibilidad)
print(mainshocks_chile[mainshocks_chile$magnitud >= 8.0,
c("id", "magnitud", "N_replicas")])
modelo_chile_nb_sin_M8 <- glm.nb(N_replicas ~ magnitud, data = chile_sin_M8)
coefs_sin_M8 <- summary(modelo_chile_nb_sin_M8)$coefficients
IRR_sin_M8 <- exp(coefs_sin_M8["magnitud", "Estimate"])
ic_inf_sin_M8 <- exp(coefs_sin_M8["magnitud", "Estimate"] -
qnorm(0.975) * coefs_sin_M8["magnitud", "Std. Error"])
ic_sup_sin_M8 <- exp(coefs_sin_M8["magnitud", "Estimate"] +
qnorm(0.975) * coefs_sin_M8["magnitud", "Std. Error"])
cat("\n=== Chile: IRR de magnitud, excluyendo M>=8.0 ===\n")
cat("N =", nrow(chile_sin_M8), "(se excluyeron",
nrow(mainshocks_chile) - nrow(chile_sin_M8), "observaciones)\n")
cat("IRR =", round(IRR_sin_M8, 3),
"| IC 95% = [", round(ic_inf_sin_M8, 3), ",", round(ic_sup_sin_M8, 3), "]\n")
### Comparación final: modelo completo vs. las dos sensibilidades ----
# Chile: comparación del IRR de magnitud ante distintos criterios de exclusión
tabla_sensibilidad_chile <- data.frame(
criterio = c("Catálogo completo", "Sin obs. influyentes (Cook)", "Sin M>=8.0"),
N = c(nrow(mainshocks_chile), nrow(chile_sin_influyentes), nrow(chile_sin_M8)),
IRR = round(c(tabla_irr_chile$IRR[tabla_irr_chile$coeficiente == "magnitud"],
IRR_sin_influyentes, IRR_sin_M8), 3),
IC_95_inf = round(c(tabla_irr_chile$IC_95_inf[tabla_irr_chile$coeficiente == "magnitud"],
ic_inf_sin_influyentes, ic_inf_sin_M8), 3),
IC_95_sup = round(c(tabla_irr_chile$IC_95_sup[tabla_irr_chile$coeficiente == "magnitud"],
ic_sup_sin_influyentes, ic_sup_sin_M8), 3)
)
print(tabla_sensibilidad_chile, row.names = FALSE)
### Diagnóstico: Verificacion replicas Maule ----
evento_maule <- mainshocks_chile[which.max(mainshocks_chile$magnitud), ]
# IDs a verificar
ids_sospechosos <- c("usp000h7zt", "usp000h7ru", "usp000h7uh", "usp000h7rk",
"usp000h7sy", "usp000h7vb", "usp000h80w")
eventos_sospechosos <- mainshocks_chile[mainshocks_chile$id %in% ids_sospechosos, ]
phi1 <- evento_maule$latitude * pi / 180
phi2 <- eventos_sospechosos$latitude * pi / 180
dphi <- (eventos_sospechosos$latitude - evento_maule$latitude) * pi / 180
dlambda <- (eventos_sospechosos$longitude - evento_maule$longitude) * pi / 180
a <- sin(dphi / 2)^2 + cos(phi1) * cos(phi2) * sin(dlambda / 2)^2
distancia_al_maule <- 6371 * 2 * atan2(sqrt(a), sqrt(1 - a))
L_maule <- 10^(0.1238 * evento_maule$magnitud + 0.983)
tabla_diagnostico <- data.frame(
id = eventos_sospechosos$id,
magnitud = eventos_sospechosos$magnitud,
N_replicas_propio = eventos_sospechosos$N_replicas,
distancia_km = round(distancia_al_maule, 1),
ventana_L_maule_km = round(L_maule, 1),
fuera_de_ventana = distancia_al_maule > L_maule
)
print(tabla_diagnostico, row.names = FALSE)
# Corrección: fusión de sub-clusters fragmentados con la secuencia del Maule ----
id_maule <- evento_maule$id
ids_impostores <- c("usp000h7rk", "usp000h7ru", "usp000h7sy", "usp000h7uh",
"usp000h7vb", "usp000h7zt", "usp000h80w")
# Numero de replicas del Maule antes de corregir
sum(chile$cluster_id == id_maule & chile$es_dependiente)
chile$cluster_id[chile$cluster_id %in% ids_impostores] <- id_maule
chile$es_dependiente[chile$id %in% ids_impostores] <- TRUE
#N réplicas del Maule después de corregir
sum(chile$cluster_id == id_maule & chile$es_dependiente)
### Reconstrucción de numero de replicas con el catálogo corregido ----
mainshocks_chile_corregido <- chile[!chile$es_dependiente, ]
conteo_replicas_chile_corr <- table(chile$cluster_id[chile$es_dependiente])
mainshocks_chile_corregido$N_replicas <- 0
posicion <- match(mainshocks_chile_corregido$id, names(conteo_replicas_chile_corr))
mainshocks_chile_corregido$N_replicas[!is.na(posicion)] <-
as.numeric(conteo_replicas_chile_corr[posicion[!is.na(posicion)]])
# Chile (corregido): distribución de numero de réplicas por sismo principal
summary(mainshocks_chile_corregido$N_replicas)
cat("N mainshocks -- original:", nrow(mainshocks_chile),
"| corregido:", nrow(mainshocks_chile_corregido), "\n")
cat("N réplicas del Maule -- original:",
mainshocks_chile$N_replicas[mainshocks_chile$id == id_maule],
"| corregido:",
mainshocks_chile_corregido$N_replicas[mainshocks_chile_corregido$id == id_maule], "\n")
### Reajuste del modelo con el catálogo corregido ----
modelo_prod_chile_pois_corr <- glm(N_replicas ~ magnitud, data = mainshocks_chile_corregido,
family = poisson(link = "log"))
dispersion_chile_corr <- sum(residuals(modelo_prod_chile_pois_corr, type = "pearson")^2) /
modelo_prod_chile_pois_corr$df.residual
# Chile (corregido) - dispersión Poisson (Pearson/df)
round(dispersion_chile_corr, 3)
modelo_prod_chile_nb_corr <- glm.nb(N_replicas ~ magnitud, data = mainshocks_chile_corregido)
summary(modelo_prod_chile_nb_corr)
cat("AIC -- Poisson:", round(AIC(modelo_prod_chile_pois_corr), 1),
"| Binomial Negativa:", round(AIC(modelo_prod_chile_nb_corr), 1), "\n")
coefs_corr <- summary(modelo_prod_chile_nb_corr)$coefficients
IRR_corr <- exp(coefs_corr["magnitud", "Estimate"])
ic_inf_corr <- exp(coefs_corr["magnitud", "Estimate"] - qnorm(0.975) * coefs_corr["magnitud", "Std. Error"])
ic_sup_corr <- exp(coefs_corr["magnitud", "Estimate"] + qnorm(0.975) * coefs_corr["magnitud", "Std. Error"])
cat("\nChile (corregido) -- IRR magnitud:", round(IRR_corr, 3),
"| IC 95% = [", round(ic_inf_corr, 3), ",", round(ic_sup_corr, 3), "]\n")
### Diagnóstico de influencia sobre el modelo corregido ----
cooks_chile_corr <- cooks.distance(modelo_prod_chile_nb_corr)
umbral_cook_corr <- 4 / nrow(mainshocks_chile_corregido)
influyentes_corr <- which(cooks_chile_corr > umbral_cook_corr)
# Chile (corregido) - N° de observaciones influyentes (Cook)
length(influyentes_corr)
nrow(mainshocks_chile_corregido)
chile_corr_sin_influyentes <- mainshocks_chile_corregido[-influyentes_corr, ]
modelo_chile_nb_corr_sin_influyentes <- glm.nb(N_replicas ~ magnitud,
data = chile_corr_sin_influyentes)
coefs_sin_inf_corr <- summary(modelo_chile_nb_corr_sin_influyentes)$coefficients
IRR_sin_inf_corr <- exp(coefs_sin_inf_corr["magnitud", "Estimate"])
ic_inf_sin_inf_corr <- exp(coefs_sin_inf_corr["magnitud", "Estimate"] -
qnorm(0.975) * coefs_sin_inf_corr["magnitud", "Std. Error"])
ic_sup_sin_inf_corr <- exp(coefs_sin_inf_corr["magnitud", "Estimate"] +
qnorm(0.975) * coefs_sin_inf_corr["magnitud", "Std. Error"])
cat("Chile (corregido, sin influyentes) -- IRR:", round(IRR_sin_inf_corr, 3),
"| IC 95% = [", round(ic_inf_sin_inf_corr, 3), ",", round(ic_sup_sin_inf_corr, 3), "]\n")
### Comparación final: original (fragmentado) vs. corregido (fusionado) ----
# Comparación: antes vs. después de corregir la fragmentación del Maule
tabla_comparacion_correccion <- data.frame(
escenario = c("Original (fragmentado)",
"Corregido (fusionado)",
"Corregido, sin influyentes (Cook)"),
N_mainshocks = c(nrow(mainshocks_chile),
nrow(mainshocks_chile_corregido),
nrow(chile_corr_sin_influyentes)),
N_replicas_maule = c(mainshocks_chile$N_replicas[mainshocks_chile$id == id_maule],
mainshocks_chile_corregido$N_replicas[mainshocks_chile_corregido$id == id_maule],
NA),
IRR_magnitud = round(c(tabla_irr_chile$IRR[tabla_irr_chile$coeficiente == "magnitud"],
IRR_corr, IRR_sin_inf_corr), 3),
IC_95_inf = round(c(tabla_irr_chile$IC_95_inf[tabla_irr_chile$coeficiente == "magnitud"],
ic_inf_corr, ic_inf_sin_inf_corr), 3),
IC_95_sup = round(c(tabla_irr_chile$IC_95_sup[tabla_irr_chile$coeficiente == "magnitud"],
ic_sup_corr, ic_sup_sin_inf_corr), 3)
)
print(tabla_comparacion_correccion, row.names = FALSE)
# Regresion Logistica ----
chile$Y <- ifelse(chile$magnitud >= 6.0, 1, 0)
california$Y <- ifelse(california$magnitud >= 6.0, 1, 0)
# Convertir fecha
datos_Y_completo$time_fecha <- as.POSIXct(
datos_Y_completo$time,
tz = "UTC",
tryFormats = c(
"%Y-%m-%dT%H:%M:%OSZ",
"%Y-%m-%dT%H:%M:%OS",
"%Y-%m-%d %H:%M:%OS"
)
)
# Eliminar fechas que no pudieron convertirse
datos_Y_completo <- datos_Y_completo[
!is.na(datos_Y_completo$time_fecha),
]
# Ordenar cronológicamente
datos_Y_completo <- datos_Y_completo[
order(datos_Y_completo$time_fecha),
]
# Revisar distribución de Y
table(datos_Y_completo$Y)
prop.table(table(datos_Y_completo$Y))
acf(
datos_Y_completo$Y,
lag.max = 30,
main = "ACF de la variable Y ordenada temporalmente"
)
Box.test(
datos_Y_completo$Y,
lag = 30,
type = "Ljung-Box"
)
## Sismos de fondo
# Mantener solo eventos clasificados como independientes
datos_Y <- chile[
!is.na(chile$es_dependiente) &
chile$es_dependiente == FALSE &
!is.na(chile$Y) &
!is.na(chile$time),
]
# Convertir fecha
datos_Y$time_fecha <- as.POSIXct(
datos_Y$time,
tz = "UTC",
tryFormats = c(
"%Y-%m-%dT%H:%M:%OSZ",
"%Y-%m-%dT%H:%M:%OS",
"%Y-%m-%d %H:%M:%OS"
)
)
# Eliminar fechas que no pudieron convertirse
datos_Y <- datos_Y[!is.na(datos_Y$time_fecha), ]
# Ordenar cronológicamente
datos_Y <- datos_Y[order(datos_Y$time_fecha), ]
# Revisar variable respuesta
table(datos_Y$Y)
prop.table(table(datos_Y$Y))
acf(
datos_Y$Y,
lag.max = 30,
main = "ACF de la variable Y ordenada temporalmente"
)
Box.test(
datos_Y$Y,
lag = 30,
type = "Ljung-Box"
)
## Regresión Logistica Chile ----
datos_modelo_chile <- chile[
!is.na(chile$es_dependiente) &
chile$es_dependiente == FALSE &
complete.cases(
chile[, c(
"Y",
"depth_class",
"latitude",
"longitude",
"time"
)]
),
]
# Convertir la fecha
datos_modelo_chile$time_fecha <- as.POSIXct(
datos_modelo_chile$time,
tz = "UTC",
tryFormats = c(
"%Y-%m-%dT%H:%M:%OSZ",
"%Y-%m-%dT%H:%M:%OS",
"%Y-%m-%d %H:%M:%OS"
)
)
# Eliminar fechas que no pudieron convertirse
datos_modelo_chile <- datos_modelo_chile[
!is.na(datos_modelo_chile$time_fecha),
]
# Ordenar cronológicamente
datos_modelo_chile <- datos_modelo_chile[
order(datos_modelo_chile$time_fecha),
]
# Asegurar que la profundidad sea categórica
datos_modelo_chile$depth_class <- factor(
datos_modelo_chile$depth_class
)
# Centrar latitud y longitud
datos_modelo_chile$latitude_c <-
datos_modelo_chile$latitude -
mean(datos_modelo_chile$latitude)
datos_modelo_chile$longitude_c <-
datos_modelo_chile$longitude -
mean(datos_modelo_chile$longitude)
# Revisar la variable respuesta
# Frecuencias de Y
print(table(datos_modelo_chile$Y))
# Proporciones de Y
print(prop.table(table(datos_modelo_chile$Y)))
# Número de observaciones utilizadas
print(nrow(datos_modelo_chile))
# Modelo logístico completo
modelo_chile <- glm(
Y ~ depth_class+ latitude_c + longitude_c ,
data = datos_modelo_chile,
family = binomial(link = "logit")
)
summary(modelo_chile)
# Modelo nulo
modelo_nulo_chile <- glm(
Y ~ 1,
data = datos_modelo_chile,
family = binomial(link = "logit")
)
# Significancia global
# Comparación entre modelo nulo y modelo completo
anova(
modelo_nulo_chile,
modelo_chile,
test = "Chisq"
)
# Aporte individual de cada variable
drop1(
modelo_chile,
test = "Chisq"
)
modelo_chile<- step(
modelo_chile,
direction = "backward"
)
summary(modelo_chile)
anova(
modelo_nulo_chile,
modelo_chile,
test = "Chisq"
)
# Multicolinealidad
# Factores de inflación de varianza
vif_chile <- car::vif(modelo_chile)
print(vif_chile)
# Independencia temporal de los residuos
residuos_chile <- residuals(
modelo_chile,
type = "pearson"
)
# Como la base ya está ordenada por fecha,
# los residuos también quedan en orden cronológico
lag_chile <- min(
30,
length(residuos_chile) - 1
)
acf(
residuos_chile,
lag.max = lag_chile,
main = "ACF de residuos de Pearson - Chile"
)
prueba_temporal_chile <- Box.test(
residuos_chile,
lag = lag_chile,
type = "Ljung-Box"
)
cat("\nPrueba de Ljung-Box:\n")
print(prueba_temporal_chile)
# Odds ratios
coeficientes_chile <- coef(modelo_chile)
errores_chile <- sqrt(
diag(vcov(modelo_chile))
)
odds_ratio_chile <- data.frame(
OR = exp(coeficientes_chile),
LI_95 = exp(
coeficientes_chile -
1.96 * errores_chile
),
LS_95 = exp(
coeficientes_chile +
1.96 * errores_chile
)
)
# Odds ratios e intervalos de confianza
print(
round(
odds_ratio_chile,
4
)
)
# Revisar convergencia
print(modelo_chile$converged)
# Rango de probabilidades estimadas
print(range(fitted(modelo_chile)))
### Pseudo R2
R2_McFadden <- 1 -
as.numeric(logLik(modelo_chile)) /
as.numeric(logLik(modelo_nulo_chile))
R2_McFadden
pR2(modelo_chile)
## Regresion Logistica California ----
datos_modelo_california <- california[
!is.na(california$es_dependiente) &
california$es_dependiente == FALSE &
complete.cases(
california[, c(
"Y",
"depth_class",
"latitude",
"longitude",
"time"
)]
),
]
# Convertir la fecha
datos_modelo_california$time_fecha <- as.POSIXct(
datos_modelo_california$time,
tz = "UTC",
tryFormats = c(
"%Y-%m-%dT%H:%M:%OSZ",
"%Y-%m-%dT%H:%M:%OS",
"%Y-%m-%d %H:%M:%OS"
)
)
# Eliminar fechas que no pudieron convertirse
datos_modelo_california <- datos_modelo_california[
!is.na(datos_modelo_california$time_fecha),
]
# Ordenar cronológicamente
datos_modelo_california <- datos_modelo_california[
order(datos_modelo_california$time_fecha),
]
# Asegurar que la profundidad sea categórica
datos_modelo_california$depth_class <- factor(
datos_modelo_california$depth_class
)
# Centrar latitud y longitud
datos_modelo_california$latitude_c <-
datos_modelo_california$latitude -
mean(datos_modelo_california$latitude)
datos_modelo_california$longitude_c <-
datos_modelo_california$longitude -
mean(datos_modelo_california$longitude)
# Revisar la variable respuesta
# Frecuencias de Y
print(table(datos_modelo_california$Y))
# Proporciones de Y
print(prop.table(table(datos_modelo_california$Y)))
# Número de observaciones utilizadas
print(nrow(datos_modelo_california))
# Modelo logístico completo
modelo_california <- glm(
Y ~ depth_class + latitude_c + longitude_c,
data = datos_modelo_california,
family = binomial(link = "logit")
)
summary(modelo_california)
# Modelo nulo
modelo_nulo_california <- glm(
Y ~ 1,
data = datos_modelo_california,
family = binomial(link = "logit")
)
# Significancia global
anova(
modelo_nulo_california,
modelo_california,
test = "Chisq"
)
# Aporte individual de cada variable
drop1(
modelo_california,
test = "Chisq"
)
# Multicolinealidad
vif_california <- car::vif(
modelo_california
)
print(vif_california)
# Independencia temporal de los residuos
residuos_california <- residuals(
modelo_california,
type = "pearson"
)
lag_california <- min(
30,
length(residuos_california) - 1
)
acf(
residuos_california,
lag.max = lag_california,
main = "ACF de residuos de Pearson - California"
)
prueba_temporal_california <- Box.test(
residuos_california,
lag = lag_california,
type = "Ljung-Box"
)
# Prueba de Ljung-Box
print(prueba_temporal_california)
# Odds ratios
coeficientes_california <- coef(
modelo_california
)
errores_california <- sqrt(
diag(vcov(modelo_california))
)
odds_ratio_california <- data.frame(
OR = exp(coeficientes_california),
LI_95 = exp(
coeficientes_california -
1.96 * errores_california
),
LS_95 = exp(
coeficientes_california +
1.96 * errores_california
)
)
# Odds ratios e intervalos de confianza
print(
round(
odds_ratio_california,
4
)
)
# Revisar convergencia
print(modelo_california$converged)
# Rango de probabilidades estimadas
print(range(fitted(modelo_california)))
# Ley Gutenber-Ritcher ----
calcular_GR <- function(
datos,
columna_mag,
region,
solo_fondo = FALSE,
Mc = 5,
deltaM = 0.001
) {
magnitudes <- datos[[columna_mag]]
if (solo_fondo) {
magnitudes <- magnitudes[
!is.na(datos$es_dependiente) &
datos$es_dependiente == FALSE
]
catalogo <- "Sin replicas"
} else {
catalogo <- "Con replicas"
}
magnitudes <- magnitudes[
!is.na(magnitudes) &
magnitudes >= Mc
]
N <- length(magnitudes)
if (N < 2) {
stop("No hay suficientes eventos con magnitud mayor o igual a Mc.")
}
media_M <- mean(magnitudes)
# Estimación de b por máxima verosimilitud
b <- log10(exp(1)) /
(
media_M -
(Mc - deltaM / 2)
)
# Cálculo de a a partir de b y N
a <- log10(N) + b * Mc
data.frame(
region = region,
catalogo = catalogo,
variable_magnitud = columna_mag,
Mc = Mc,
deltaM = deltaM,
N = N,
media_M = media_M,
a = a,
b = b,
metodo_b = "Maxima verosimilitud",
ecuacion = paste0(
"log10(N) = ",
round(a, 4),
" - ",
round(b, 4),
" M"
)
)
}
GR_chile_con_replicas <- calcular_GR(
datos = chile,
columna_mag = "Mw_hom",
region = "Chile",
solo_fondo = FALSE
)
GR_chile_sin_replicas <- calcular_GR(
datos = chile,
columna_mag = "Mw_hom",
region = "Chile",
solo_fondo = TRUE
)
GR_california_con_replicas <- calcular_GR(
datos = california,
columna_mag = "mag",
region = "California",
solo_fondo = FALSE
)
GR_california_sin_replicas <- calcular_GR(
datos = california,
columna_mag = "mag",
region = "California",
solo_fondo = TRUE
)
resultados_GR <- rbind(
GR_chile_con_replicas,
GR_chile_sin_replicas,
GR_california_con_replicas,
GR_california_sin_replicas
)
resultados_GR
# Dependencia Temporal Mensual ----
# conteos mensuales
construir_conteo_mensual_total <- function(df,
anio_inicio = 2000,
anio_fin = 2025) {
df$year_month <- format(df$time, "%Y-%m")
todos_los_meses <- format(
seq(as.Date(paste0(anio_inicio, "-01-01")),
as.Date(paste0(anio_fin, "-12-01")),
by = "month"),
"%Y-%m"
)
conteo <- df %>%
count(year_month) %>%
right_join(data.frame(year_month = todos_los_meses),
by = "year_month") %>%
mutate(n = ifelse(is.na(n), 0, n)) %>%
arrange(year_month)
conteo
}
#construir todos los conteos
conteo_mensual_total_chile <- construir_conteo_mensual_total(chile)
conteo_mensual_total_california <- construir_conteo_mensual_total(california)
conteo_mensual_total_chile$pais <- "Chile"
conteo_mensual_total_california$pais <- "California"
print(conteo_mensual_total_chile, n = Inf)
print(conteo_mensual_total_california, n = Inf)
# Serie mensual Chile
ts_mensual_total_chile <- ts(
conteo_mensual_total_chile$n,
start = c(2000, 1),
frequency = 12
)
# Serie mensual California
ts_mensual_total_california <- ts(
conteo_mensual_total_california$n,
start = c(2000, 1),
frequency = 12
)
#graficar Chile
conteo_mensual_total_chile <- conteo_mensual_total_chile %>%
mutate(year_month = as.Date(paste0(year_month, "-01")))
ggplot(conteo_mensual_total_chile,
aes(x = year_month, y = n)) +
geom_line(color = "#DC143C", linewidth = 0.7) +
labs(
title = "Conteo mensual de sismos M ≥ 5",
subtitle = "Chile",
x = "Año",
y = "Número de sismos"
) +
theme_minimal(base_size = 12)
#graficar California
conteo_mensual_total_california <- conteo_mensual_total_california %>%
mutate(year_month = as.Date(paste0(year_month, "-01")))
ggplot(conteo_mensual_total_california,
aes(x = year_month, y = n)) +
geom_line(linewidth = 0.7, color = "#1f78b4") +
scale_x_date(
date_breaks = "5 years",
date_labels = "%Y"
) +
labs(
title = "Conteo mensual de sismos M ≥ 5",
subtitle = "California",
x = "Año",
y = "Número de sismos"
) +
theme_minimal(base_size = 12)
#ACF Chile
acf(ts_mensual_total_chile,
main = "ACF conteos mensuales de sismos M≥5 - Chile",
xlab = "Rezago (meses)")
# PACF Chile
pacf(
ts_mensual_total_chile,
main = "PACF Conteos mensuales de sismos M≥5 - Chile",
xlab = "Rezago (meses)"
)
#ACF California
acf(ts_mensual_total_california,
main = "ACF conteos mensuales de sismos M≥5 - California",
xlab = "Rezago (meses)")
# PACF California
pacf(
ts_mensual_total_california,
main = "PACF Conteos mensuales de sismos M≥5 - California",
xlab = "Rezago (meses)"
)
#Box-Ljung test Chile
Box.test(
ts_mensual_total_chile,
lag =24,
type = "Ljung-Box"
)
#Box-Ljung test California
Box.test(
ts_mensual_total_california,
lag =24,
type = "Ljung-Box"
)