Colombia se encuentra en un régimen tectónico complejo, ya que está cerca de la unión de tres placas tectónicas: la Sudamericana, la de Nazca y la del Caribe (Cooper et al., 1005; Cortés y Colletta, 2005; Villagómez et al., 2011; Escobar-Rey et al., 2024). Como consecuencia, se generan diferentes esfuerzos que dan lugar a movimientos telúricos, como terremotos o sismos (Vergara et al., 1996; Guzman et al., 1998; Restrepo-Pace, 1999; Taboada et al., 2000).
Históricamente, Colombia se ha visto afectada por eventos sísmicos significativos. Por ejemplo, el 18 de octubre de 1992 ocurrió un sismo de magnitud 6.6 Mw y otro de 7.1 Mw en la ciudad de Bogotá. El 25 de enero de 1999, un sismo de magnitud 6.1 Mw afectó el Eje Cafetero. Asimismo, en el departamento del Huila se registraron sismos en 1967 y 2016 con magnitudes de 7.0 Mw y 5.3 Mw, respectivamente (Ramírez, 1975; Espinosa, 2012; Servicio Geológico Colombiano, s.f.).
En esta investigación, se analizará la variabilidad de la sismicidad en Colombia mediante pruebas estadísticas que permitan identificar diferencias significativas en la magnitud y la profundidad de los sismos según su distribución regional y temporal. Para ello, se plantean las siguientes hipótesis: 1) Determinar si hay diferencias significativas en la magnitud promedio entre las cinco regiones de Colombia. 2) Comparar la magnitud promedio entre las regiones Andina y Pacífica. 3) Evaluar si existe una relación lineal entre la profundidad y la magnitud de los sismos. 4) Determinar se hay diferencias significativas en la media de la magnitud entre diferentes años.
Los datos de la investigación fueron tomados del catálogo de sismicidad del Servicio Geológico Colombiano, el cual tiene dos catálogos. El primero corresponde al período comprendido entre el 1 de junio de 1993 y el 28 de febrero de 2018; este incluye todos los eventos localizados en el territorio nacional, registrados por una red de estaciones sísmicas y adquiridos a través del software Earthworm, luego procesados con el software SEISAN (ervicio Geológico Colombiano, s.f).
El segundo catálogo abarca desde el 1 de marzo de 2018 hasta la fecha; los eventos reportados en este período han sido adquiridos y procesados utilizando el software SeisComP3, desarrollado por GFZ (GeoForschungsZentrum Potsdam) y GEMPA. Actualmente, la red del SGC cuenta con 159 estaciones sísmicas distribuidas a lo largo del país, encargadas de registrar los datos. Todos los eventos reportados en ambos catálogos han sido revisados manualmente por analistas de sismología del SGC (servicio Geológico Colombiano, s.f).
Para el análisis realizado, se utilizó el segundo catálogo, comprendido entre el 1 de marzo de 2018 y el 31 de diciembre de 2024
Mag: Magnitud del sismo. Es una medida de la energía liberada durante un sismo. Se expresa en una escala logarítmica, como la escala de Richter.
Región: Región geográfica de Colombia donde se encuentra el departamento. Las regiones pueden incluir Andina, Caribe, Pacífica, etc.
Año: Año en el que ocurrió el sismo.
## Prof Mag Región año
## Min. : 0.00 Min. :0.100 Length:160224 Min. :2018
## 1st Qu.: 25.29 1st Qu.:1.400 Class :character 1st Qu.:2019
## Median :120.65 Median :1.700 Mode :character Median :2021
## Mean : 92.02 Mean :1.719 Mean :2021
## 3rd Qu.:141.70 3rd Qu.:1.900 3rd Qu.:2023
## Max. :319.00 Max. :6.100 Max. :2024
# Resumen agrupado por regiones
datos_sismico %>%
group_by(Región) %>%
summarise(n = length(Mag),
media = mean(Mag),
ds = sd(Mag),
mediana = median(Mag),
minimo = min(Mag),
maximo = max(Mag)) %>%
mutate(proporcion = n / sum(n)) %>%
kable(digits = 2, col.names = c("Región",
"Cantidad", "Media", "Desv. Est.", "Mediana",
"Mínimo", "Máximo", "Proporción"))
Región | Cantidad | Media | Desv. Est. | Mediana | Mínimo | Máximo | Proporción |
---|---|---|---|---|---|---|---|
Amazonia | 45 | 2.12 | 0.37 | 2.1 | 1.3 | 3.2 | 0.00 |
Andina | 126062 | 1.75 | 0.49 | 1.7 | 0.1 | 5.6 | 0.79 |
Caribe | 10204 | 1.60 | 0.46 | 1.6 | 0.2 | 5.8 | 0.06 |
Orinoquía | 9713 | 1.50 | 0.53 | 1.4 | 0.1 | 6.1 | 0.06 |
Pacífica | 14200 | 1.70 | 0.53 | 1.6 | 0.1 | 6.1 | 0.09 |
datos_sismico%>%
ggplot(aes(x=Región, y=Mag))+
geom_boxplot()+
stat_summary(fun = mean, geom = "point", color = "red")+
theme_bw()
### Resumen agrupado por años
datos_sismico %>%
group_by(año) %>%
summarise(n = length(Mag),
media = mean(Mag),
ds = sd(Mag),
mediana = median(Mag),
minimo = min(Mag),
maximo = max(Mag)) %>%
mutate(proporcion = n / sum(n)) %>%
kable(digits = 2, col.names = c("Año",
"Cantidad", "Media", "Desv. Est.", "Mediana",
"Mínimo", "Máximo", "Proporción"))
Año | Cantidad | Media | Desv. Est. | Mediana | Mínimo | Máximo | Proporción |
---|---|---|---|---|---|---|---|
2018 | 16093 | 1.71 | 0.52 | 1.7 | 0.1 | 5.6 | 0.10 |
2019 | 24188 | 1.65 | 0.51 | 1.6 | 0.1 | 6.1 | 0.15 |
2020 | 24644 | 1.68 | 0.50 | 1.6 | 0.2 | 5.8 | 0.15 |
2021 | 21700 | 1.74 | 0.49 | 1.7 | 0.2 | 5.1 | 0.14 |
2022 | 21019 | 1.75 | 0.49 | 1.7 | 0.1 | 5.5 | 0.13 |
2023 | 24604 | 1.73 | 0.51 | 1.7 | 0.1 | 6.1 | 0.15 |
2024 | 27976 | 1.78 | 0.46 | 1.7 | 0.2 | 5.6 | 0.17 |
datos_sismico%>%
ggplot(aes(x=año, y=Mag))+
geom_boxplot()+
stat_summary(fun = mean, geom = "point", color = "red")+
theme_bw()
ggplot(datos_sismico, aes(x = Prof)) +
geom_histogram(binwidth = 10, fill = "skyblue", color = "black", alpha = 0.7) +
theme_minimal() +
labs(title = "Distribución de Prof", x = "Profundidad - km", y = "Frecuencia") +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"), # Centrar y aumentar el título
axis.title.x = element_text(size = 14, face = "bold"), # Aumentar y negrita en el título del eje x
axis.title.y = element_text(size = 14, face = "bold"), # Aumentar y negrita en el título del eje y
axis.text.x = element_text(size = 12), # Aumentar el tamaño de las etiquetas del eje x
axis.text.y = element_text(size = 12), # Aumentar el tamaño de las etiquetas del eje y
panel.grid = element_blank() # Quitar fondo cuadriculado
) +
geom_vline(aes(xintercept = mean(Prof)), color = "red", linetype = "dashed", size = 1) + # Línea de la media
annotate("text", x = mean(datos_sismico$Prof), y = 3, label = paste("Media:", round(mean(datos_sismico$Prof), 2)), color = "red", vjust = -1)
ggplot(resumen, aes(x = Región, y = media, fill = Región)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_errorbar(aes(ymin = media - se, ymax = media + se), width = 0.2, position = position_dodge(0.9)) +
theme_minimal() +
labs(title = "Medias de Magnitudes por Región con Errores Estándar", x = "Región", y = "Magnitud Media") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(resumen, aes(x = año, y = media, fill = año)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_errorbar(aes(ymin = media - se, ymax = media + se), width = 0.2, position = position_dodge(0.9)) +
theme_minimal() +
labs(title = "Medias de Magnitudes por año con Errores Estándar", x = "Año", y = "Magnitud Media") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Prueba 1: Determinar si hay diferencias significativas en la magnitud promedio entre las cinco regiones de Colombia. Ho: Miu_Amazonia = Miu_Andina = Miu_Caribe = Miu_Orinoquia = Miu_Pacifica H1: Miu_i /= Miu_j
# ANOVA
datos_sismico$Región <- as.factor(datos_sismico$Región) # Se convierte a factor
str(datos_sismico$Región)
## Factor w/ 5 levels "Amazonia","Andina",..: 2 2 2 2 2 2 2 3 2 2 ...
anova_result <- aov(Mag ~ Región, data = datos_sismico)
residuos <- residuals(anova_result)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## Región 4 716 179.10 738.3 <2e-16 ***
## Residuals 160219 38868 0.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Supuestos:
# 1. Normalidad
# Prueba de Kolmogorov-Smirnof sobre los residuos
ks.test(residuos, "pnorm", mean = mean(residuos), sd = sd(residuos))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.10659, p-value < 2.2e-16
## alternative hypothesis: two-sided
p-value < 0.05 No normal, dado que hay ties en los datos, se usa la prueba lillierfors.
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.10659, p-value < 2.2e-16
p-value < 0.05 No normal
Gráfico Q-Q
p-value < 2.2e-16 se rechaza Ho (varianzas no homogéneas)
Valores de autocorrelación bajos y dentro de los intervalos de
confianza, los datos se consideran independientes en el tiempo.
Dado que los residuos no se distribuyen normalmente y las varianzas no son homogéneas, pero sí hay independencia entre las variables se usa una prueba no paramétrica Kruskal Wallis Boxplot
ggplot(datos_sismico, aes(x = Región, y = Mag, fill = Región)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Distribución de Magnitudes por Región",
x = "Región",
y = "Magnitud")
Las medianas no coinciden entre regiones y los rangos intercuartílicos
son diferentes entre regiones, para confirmar si hay diferencias
significativas entre las medianas se procede con la prueba
##
## Kruskal-Wallis rank sum test
##
## data: Mag by Región
## Kruskal-Wallis chi-squared = 4359.5, df = 4, p-value < 2.2e-16
p-value 2.26 e-16, se rechaza Ho, existen diferencias significativas entre las magnitudes sísmicas de al menos una las cinco regiones de Colombia. Prueba post Hoc
## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Bonferroni method.
## Comparison Z P.unadj P.adj
## 1 Amazonia - Andina 5.866692 4.445743e-09 4.445743e-08
## 2 Amazonia - Caribe 8.099778 5.505952e-16 5.505952e-15
## 3 Andina - Caribe 32.586409 6.389598e-233 6.389598e-232
## 4 Amazonia - Orinoquía 9.946359 2.615784e-23 2.615784e-22
## 5 Andina - Orinoquía 58.064341 0.000000e+00 0.000000e+00
## 6 Caribe - Orinoquía 19.472743 1.869871e-84 1.869871e-83
## 7 Amazonia - Pacífica 6.985313 2.842217e-12 2.842217e-11
## 8 Andina - Pacífica 19.006964 1.493489e-80 1.493489e-79
## 9 Caribe - Pacífica -12.879339 5.883702e-38 5.883702e-37
## 10 Orinoquía - Pacífica -33.658332 2.354961e-248 2.354961e-247
Todos los p-value < 0.05, lo que indica que todas las regiones tienen magnitudes significativamente distintas entre sí
Prueba 2: Comparar la magnitud promedio entre las regiones Andina y Pacífica. Ho: MiuAndina - MiuPacífica es igual a cero H1: MiuAndina - MiuPacífica es diferente de cero
En primer lugar,se filtran los datos solo por región Andina y Paçífica
Supuestos: Normalidad:Se rechaza Ho, no normales.Por ende, no se puede hacer uso de una prueba T Student para la diferencia de medias
Homogeneidad de Varianzas: se rechaza Ho, las varianzas no son homogéneas, así que refuerza el ítem anterior.
Se procede a aplicar una prueba no paramétrica, para este caso, se usa la prueba no paramétrica Mann - Whitney - Wilcoxon
Esta prueba ordena todos los valores de menor a mayor y luego asigna rangos a cada valor, suma estos rangos y finalmente calcula el estadístico W
datos_sismico_prueba_2 <- datos_sismico_prueba_2 %>%
mutate(Rango = rank(Mag))
head(datos_sismico_prueba_2)
Se suman los rangos por región
suma_rangos <- datos_sismico_prueba_2 %>%
group_by(Región) %>%
summarise(Suma_Rangos = sum(Rango))
print(suma_rangos)
## # A tibble: 2 × 2
## Región Suma_Rangos
## <fct> <dbl>
## 1 Andina 8927360461
## 2 Pacífica 909423992
Prueba Mann-Whitney-Wilcoxon
prueba_wilcoxon <- wilcox.test(Mag ~ Región, data = datos_sismico_prueba_2, exact = FALSE)
prueba_wilcoxon
##
## Wilcoxon rank sum test with continuity correction
##
## data: Mag by Región
## W = 981483508, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
p-value < 0.05, se rechaza la hipótesis nula, esto significa que la magnitud promedio entre las regiones Andina y pacífica es significativamente diferente.
En conclusión, aunque la diferencia en la mediana es pequeña (0.1), la prueba de Wilcoxon determina que es estadísticamente significativa.
Prueba 3: Evaluar si existe una relación lineal entre la profundidad (Prof) y la magnitud (Mag) de los sismos. Ho: B1 = 0 H1: B1 /= 0
datos_sismico <- na.omit(datos_sismico)
#Se estandarizan la profundidad y la magnitud en primer lugar:
datos_sismico$Prof_scaled <- scale(datos_sismico$Prof)
datos_sismico$Mag_scaled <- scale(datos_sismico$Mag)
Gráfico de dispersión
ggplot(datos_sismico, aes(x = Prof_scaled, y = Mag_scaled)) +
geom_point(color = "blue", alpha = 0.5) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Relación entre Profundidad y Magnitud (Escaladas)",
x = "Profundidad (escalada)",
y = "Magnitud (escalada)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
Modelo de regresión lineal
##
## Call:
## lm(formula = Mag_scaled ~ Prof_scaled, data = datos_sismico)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3618 -0.5677 -0.1460 0.4201 9.3009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.760e-15 2.341e-03 0.0 1
## Prof_scaled 3.495e-01 2.341e-03 149.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9369 on 160222 degrees of freedom
## Multiple R-squared: 0.1221, Adjusted R-squared: 0.1221
## F-statistic: 2.229e+04 on 1 and 160222 DF, p-value: < 2.2e-16
Supuesto de normalidad del modelo:
los residuos no son normales Dado que hay ties en el ks.test, se aplica la lillie.test
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.098233, p-value < 2.2e-16
p-value < 0.05, se rechaza Ho, se descarta normalidad.
Ya que los residuos del modelo lineal no son normales, se hace uso de pruebas no paramétricas para confirmar si las variables se encuentran correlacionadas Ho : p = 0 H1: p /= 0
##
## Spearman's rank correlation rho
##
## data: datos_sismico$Prof and datos_sismico$Mag
## S = 4.0745e+14, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.405656
rho es 0.406 lo que indica una correlación moderadamente positiva y p-value es menor a 0,05, se rechaza Ho, lo que indica que existe una relación significativa entre la profundidad y la magnitud.
Prueba 4: Determinar se hay diferencias
significativas en la media de la magnitud entre diferentes los distintos
años.
Ho: Miu_2018 = Miu_2019 = Miu_2020 = Miu_2021 = Miu_2022 = Miu_2023 =
Miu_2024 H1: Miu_i /= Miu_j
## Factor w/ 7 levels "2018","2019",..: 1 1 1 1 1 1 1 1 1 1 ...
anova_result <- aov(Mag ~ año, data = datos_sismico)
residuos <- residuals(anova_result)
summary(anova_result)
## Df Sum Sq Mean Sq F value Pr(>F)
## año 6 276 46.05 187.7 <2e-16 ***
## Residuals 160217 39308 0.25
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Supuestos:
# 1. Normalidad
# Prueba de Kolmogorov-Smirnof sobre los residuos
ks.test(residuos, "pnorm", mean = mean(residuos), sd = sd(residuos))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: residuos
## D = 0.084558, p-value < 2.2e-16
## alternative hypothesis: two-sided
p-value < 0.05 No normal, dado que hay ties en los datos, se usa la prueba lillierfors.
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.084558, p-value < 2.2e-16
p-value < 0.05 No normal Gráfico Q-Q
Homogeneidad Dado que los datos no son normales, se usa la Prueba de
Levene para la homocedeasticidad de varianzas
p-value < 2.2e-16 se rechaza Ho (varianzas no homogéneas)
Valores de autocorrelación bajos y dentro de los intervalos de
confianza, los datos se consideran independientes en el tiempo.
Dado que los residuos no se distribuyen normalmente y las varianzas no son homogéneas, pero sí hay independencia entre las variables se usa una prueba no paramétrica Kruskal Wallis Boxplot
ggplot(datos_sismico, aes(x = año, y = Mag, fill = año)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Distribución de Magnitudes por Año",
x = "Año",
y = "Magnitud")
Las medianas no coinciden entre los años y los rangos intercuartílicos son diferentes, para confirmar si hay diferencias significativas entre las medianas se procede con la prueba Kruskal Wallis
##
## Kruskal-Wallis rank sum test
##
## data: Mag by año
## Kruskal-Wallis chi-squared = 1791.2, df = 6, p-value < 2.2e-16
p-value: < 2.2e-16 indica que la probabilidad de observar una diferencia tan extrema como la observada, bajo la hipótesis nula de que todas las medianas son iguales, es extremadamente baja. Dado que el valor p es mucho menor que 0.05, rechazamos la hipótesis nula. Prueba post Hoc
## Dunn (1964) Kruskal-Wallis multiple comparison
## p-values adjusted with the Bonferroni method.
## Comparison Z P.unadj P.adj
## 1 2018 - 2019 14.911195 2.787191e-50 5.853102e-49
## 2 2018 - 2020 7.459657 8.674797e-14 1.821707e-12
## 3 2019 - 2020 -8.405971 4.243358e-17 8.911052e-16
## 4 2018 - 2021 -4.621708 3.805929e-06 7.992452e-05
## 5 2019 - 2021 -21.364841 2.838156e-101 5.960128e-100
## 6 2020 - 2021 -13.286087 2.787820e-40 5.854423e-39
## 7 2018 - 2022 -9.027079 1.763152e-19 3.702619e-18
## 8 2019 - 2022 -26.113248 2.578522e-150 5.414897e-149
## 9 2020 - 2022 -18.122968 2.099567e-73 4.409092e-72
## 10 2021 - 2022 -4.802199 1.569326e-06 3.295585e-05
## 11 2018 - 2023 -4.852507 1.219104e-06 2.560119e-05
## 12 2019 - 2023 -22.185415 4.750282e-109 9.975592e-108
## 13 2020 - 2023 -13.847581 1.315680e-43 2.762927e-42
## 14 2021 - 2023 -0.119845 9.046059e-01 1.000000e+00
## 15 2022 - 2023 4.829197 1.370847e-06 2.878778e-05
## 16 2018 - 2024 -17.929513 6.938505e-72 1.457086e-70
## 17 2019 - 2024 -37.480037 1.947831e-307 4.090445e-306
## 18 2020 - 2024 -28.958616 2.186015e-184 4.590631e-183
## 19 2021 - 2024 -14.294726 2.360326e-46 4.956685e-45
## 20 2022 - 2024 -9.074665 1.140277e-19 2.394581e-18
## 21 2023 - 2024 -14.667196 1.045777e-48 2.196131e-47
Todos los p-value < 0.05, lo que indica que todos los años tienen magnitudes significativamente distintas entre sí
En primer lugar, al analizar las diferencias en la magnitud promedio de los sismos entre las cinco regiones de Colombia, se encontró que existen diferencias significativas, como lo confirmó la prueba de Kruskal-Wallis (p < 2.2 × 10−16). Sin embargo, al verificar los supuestos, se observó que los residuos no seguían una distribución normal y las varianzas no eran homogéneas, lo que llevó a la necesidad de utilizar métodos no paramétricos. A pesar de esto, la independencia de los datos se mantuvo, lo que permitió realizar un análisis robusto.
Posteriormente, al comparar específicamente las regiones Andina y Pacífica, se aplicó la prueba de Mann-Whitney-Wilcoxon, la cual arrojó un valor p < 2.2 × 10−16, confirmando que existe una diferencia significativa en las medianas de las magnitudes entre estas dos regiones. Aunque la diferencia en las medianas fue pequeña (0.1), esta resultó estadísticamente significativa, lo que sugiere que factores regionales podrían estar influyendo en la magnitud de los sismos.
Por otro lado, al evaluar la relación entre la profundidad y la magnitud de los sismos, se encontró que existe una correlación positiva y significativa, como lo indicó la prueba de Spearman (ρ = 0.3495, p < 2.2 × 10−16). Sin embargo, el coeficiente de determinación (R2 = 0.1221) sugiere que solo el 12.21% de la variabilidad en la magnitud puede explicarse por la profundidad, lo que indica que otros factores no considerados en este análisis podrían estar influyendo.
Finalmente, al comparar la magnitud promedio de los sismos entre diferentes años, la prueba de Kruskal-Wallis reveló diferencias significativas (p < 2.2 × 10−16). La prueba post-hoc de Dunn confirmó que la mayoría de los años presentan diferencias significativas en las medianas de las magnitudes, con excepción de la comparación entre 2021 y 2023, que no mostró diferencias significativas (p = 1.000). Esto sugiere que, aunque la magnitud de los sismos varía significativamente a lo largo del tiempo, existen períodos en los que esta magnitud se mantiene relativamente estable.
En resumen, los resultados indican que tanto la región como el año tienen un impacto significativo en la magnitud de los sismos, mientras que la profundidad también muestra una correlación, aunque moderada. Estos hallazgos resaltan la importancia de considerar múltiples factores al analizar la actividad sísmica y sugieren la necesidad de investigaciones adicionales para comprender mejor las dinámicas detrás de estos fenómenos.