AEST-1erEntregaAlg-Alg2
Alith Castillo Mejia
2023-03-02
Introducción
El Servicio Geologico de Estados Unidos (USGS) nos brinda infromación para comprender todo lo que sucede en el planeta tierra, desde los desastres naturales hasta como mejorar y proteger nuestra calidad de vida. El USGS vigila los terremotos globales, informando y evaluando el impacto y los riesgos que estos traen para la vida humana, el USGS lleva a cabo constantes investigaciones sobre la causas y los efectos de los terremotos en el planeta Tierra, toda la informacion la podemos conseguir en su pagina oficial.
Analisis explorarorio de datos (EDA)
Mientras se extraían los de datos se encontró que no todos los terremotos son naturales, algunos de ellos son provocados por el hombre. También se analizó que entre febrero y el 25 de marzo de 2019 hubo más de 8500 terremotos en todo el mundo, tambien nos dimos cuenta que los terremotos no son terremotos naturales, sino explosiones de canteras, explosiones químicas, terremotos de hielo, etc.
El primer paso es importar las librerias que contienen los datos al sistema de R:
library(tidyverse) ## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(gghighlight)
library(leaflet)
library(IRdisplay)
options(scipen = 999)
options(warn = -1)
library(lubridate)##
## Attaching package: 'lubridate'
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(viridis)## Loading required package: viridisLitelibrary(data.table)##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
##
## The following objects are masked from 'package:dplyr':
##
## between, first, last
##
## The following object is masked from 'package:purrr':
##
## transposelibrary(caTools)
library(kernlab)##
## Attaching package: 'kernlab'
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## alphalibrary(caret)## Loading required package: lattice
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## liftlibrary(ROSE)## Loaded ROSE 0.0-4library(dplyr)Tenemos que el conjunto de datos hay 22 variables que serán descritas mas adelante y 8161 datos registrados.
earthquake <- read_csv("all_month.csv")## Rows: 8161 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): magType, net, id, place, type, status, locationSource, magSource
## dbl (12): latitude, longitude, depth, mag, nst, gap, dmin, rms, horizontalE...
## dttm (2): time, updated
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.glimpse(earthquake) ## Rows: 8,161
## Columns: 22
## $ time <dttm> 2019-03-11 04:30:18, 2019-03-11 04:24:26, 2019-03-11 …
## $ latitude <dbl> 33.52100, 61.06940, 38.77333, 65.01810, 37.29130, 35.1…
## $ longitude <dbl> -116.7945, -151.8329, -122.7435, -148.7292, -117.5088,…
## $ depth <dbl> 2.82, 105.60, 1.16, 8.60, 6.90, 5.00, 6.20, 48.70, 2.5…
## $ mag <dbl> 0.44, 1.10, 0.73, 1.60, 0.70, 1.80, 1.70, 1.90, 0.56, …
## $ magType <chr> "ml", "ml", "md", "ml", "ml", "mb_lg", "ml", "ml", "md…
## $ nst <dbl> 18, NA, 5, NA, 8, NA, 30, NA, 10, NA, NA, NA, NA, NA, …
## $ gap <dbl> 53.00, NA, 213.00, NA, 176.36, 64.00, 93.14, NA, 77.00…
## $ dmin <dbl> 0.009898, NA, 0.011220, NA, 0.375000, 0.219000, 0.1370…
## $ rms <dbl> 0.1400, 0.2500, 0.0100, 0.5100, 0.2800, 0.5100, 0.3300…
## $ net <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …
## $ id <chr> "ci38263479", "ak01937u587x", "nc73150441", "ak01937u2…
## $ updated <dttm> 2019-03-11 04:33:55, 2019-03-11 04:27:53, 2019-03-11 …
## $ place <chr> "11km NE of Aguanga, CA", "51km NW of Nikiski, Alaska"…
## $ type <chr> "earthquake", "earthquake", "earthquake", "earthquake"…
## $ horizontalError <dbl> 0.25, NA, 1.55, NA, NA, 1.10, NA, NA, 0.42, 5.80, NA, …
## $ depthError <dbl> 1.02, 1.50, 0.70, 0.20, 24.50, 1.70, 0.70, 1.50, 1.31,…
## $ magError <dbl> 0.129, NA, 0.230, NA, NA, 0.255, NA, NA, NA, 0.080, NA…
## $ magNst <dbl> 11, NA, 2, NA, NA, 4, NA, NA, 1, 15, NA, NA, 11, 60, 1…
## $ status <chr> "automatic", "automatic", "automatic", "automatic", "a…
## $ locationSource <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …
## $ magSource <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …head(earthquake)## # A tibble: 6 × 22
## time latitude longi…¹ depth mag magType nst gap dmin
## <dttm> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 2019-03-11 04:30:18 33.5 -117. 2.82 0.44 ml 18 53 0.00990
## 2 2019-03-11 04:24:26 61.1 -152. 106. 1.1 ml NA NA NA
## 3 2019-03-11 04:22:42 38.8 -123. 1.16 0.73 md 5 213 0.0112
## 4 2019-03-11 04:13:42 65.0 -149. 8.6 1.6 ml NA NA NA
## 5 2019-03-11 04:10:13 37.3 -118. 6.9 0.7 ml 8 176. 0.375
## 6 2019-03-11 04:04:12 35.1 -97.6 5 1.8 mb_lg NA 64 0.219
## # … with 13 more variables: rms <dbl>, net <chr>, id <chr>, updated <dttm>,
## # place <chr>, type <chr>, horizontalError <dbl>, depthError <dbl>,
## # magError <dbl>, magNst <dbl>, status <chr>, locationSource <chr>,
## # magSource <chr>, and abbreviated variable name ¹longitudeDefinición de las variables:
Depth:Profundidad a la que comienza la ruptura del terremoto, y esta puede variar dependiendo del geoide WGS84, al nivel del mar o a la elevación promedio de las estaciones sismicas.
DepthError: Incertidumbre de la profundida a la que ocurrió el terremoto.
Dmin- Distancia horizontal desde el epicentro hasta la estación más cercana esta distancia es medida en grados en grados en donde su equivalencia es 1 grado = 111.2 kilómetros; entre mas pequeño sea este numero mas acertada está la profundidad calculada.
Gap: Mayor separación azimutal de las estaciones cercanas, entre menor sea la separación mas confiables es la posiciones horizontal calculada del terremoto. Nota: las separaciones mayores a 180° suelen tener grandes grandes incertidumbres de localización y profundidad.
HorizontalError: Incertidumbre de la localización del tereremoto, y se expresa en kilómetros.
Id: “Clave” unica que idenifica el terremoto, esta puede cambiar con el pasar del tiempo.
Latitude: latitud en grados, esta tiene valores negativos para latitudes meridionales.
locationSource: La red que creó la ubicación del terremoto.
longitude Longitud en grados, esta tiene valores negativos para longitudes occidentales.
Mag: maginitud del terremoto.
magError: Incertidumbre de la magnitud del terremoto.
MagNst: Número total de estaciones sísmicas utilizadas para calcular la magnitud de este terremoto.
MagSource: Red que originó la magnitud notificada.
MagType El método o algoritmo utilizado para calcular la magnitud.
Net- El ID de un contribuidor de datos.
Nst- Número de estaciones utilizadas para determinar la localización del terremoto.
Place: Descripción textual de la región geográficacercana al suceso. NOta: Si no hay ninguna ciudad cercana en un radio de 300 kilómetros , o si la base de datos de ciudades cercanas no está disponible por alguna razón, se utiliza el esquema de regionalización sísmica y geográfica de Flinn-Engdahl (F-E).
Status: Indica si el evento ha sido revisado por un humano, el estado puede ser automático o revisado. Nota: Los sucesos automáticos son publicados directamente por sistemas de procesamiento automático y no han sido verificados ni alterados por un humano. Los eventos revisados han sido examinados por un humano. El nivel de revisión puede variar desde una rápida comprobación de validez hasta un cuidadoso reanálisis del suceso.
Time: Hora en la que se produjo el terremoto.
type- Tipo de evento sísmico.
updated- Hora de la última actualización.
Analisis exploratorio:
options(repr.plot.width=11, repr.plot.height=8)
ggplot(earthquake, aes(time, mag)) +
geom_line(color = 'steelblue')+ theme_bw()+
ggtitle("Magnitud del Terremoto")+
xlab('Fecha')+
ylab('Magnitud')+
theme(plot.title = element_text(hjust = 0.5))
Se encontró que hay 7 terremotos que tienen una magnitud mayor a 6 en la escala de Richter.
1. Magnitud: 7.5 Fecha: 2019-02-22
2. Magnitud: 7 Fecha: 2019-03-01
3. Magnitud: 6.4 Fecha: 2019-03-06
4. Magnitud: 6.4 Fecha: 2019-02-17
5. Magnitud: 6.3 Fecha: 2019-03
6. Magnitud: 6.2 Fecha: 2019-03-1
7. Magnitud: 6.1 Fecha: 2019-03-10
Extración de caracteristicas:
Extraemos las siguientes caracteristicas:
- Ubicación del terremoto
- Hora del día
En el siguiente grafico de linea se muestra el numero de terremotos que ocurre en una hora determinada del día:
#Localización
earthquake$location <- sub('.*,\\s*','', earthquake$place)
#Hora del Día (Horas)'
earthquake$hour <- ymd_hms(earthquake$time)
earthquake$hour <- hour(earthquake$hour)
#Visualización del numero de terremotos en una hora especifica del Día
earthquake %>%
filter(!is.na(mag))%>%
group_by(hour)%>%
summarise(count = length(id),max_magnitude = max(mag))%>%
ggplot(aes(hour,count, color = max_magnitude))+geom_line()+
scale_color_viridis(direction = -1)+
scale_x_continuous(breaks=seq(0,24,1))+
xlab("Hora del Día")+
ylab("Numero de terremotos")+
ggtitle("Número de terremotos según la hora del día")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
Podemos ver como la mayor cantidad de terremotos se da en la media
noche, es decir, entre las 12 y la 1 de la mañana; a partir de ahí vemos
como los picos mas altos en sentido de cantidad se dan a las 5 de la
mañana, las 3, las 8 y las 11 de la noche respectivamente. Las horas en
las que se presentan menor numero de terremotos son las 12 del dia, las
5 de la tarde y las 9 de noche. Una particularidad, que la mayor
magnitud de toda la cantidad de terremotos se da entre las 8 y 9 de la
mañana, y tiene una cantidad aproximandamente entre los 347 y 170
registros por dia.
En el siguiente mapa vemos donde se distribuyen los terremotos y sus
magnitudes. 
Vemos como en la mayoria de las partes donde ocurren mas terremotos es
la costa Oeste de Estados Unidos y Alaska y en el cinturon del fuego del
pacifico. Para el caso de Estados Unidos y Alaska la siguiente grafica
nos muestra que el 32% de os terremotos se dan en california y en 31% se
localizan en Alaska, la magnitud promedio de estas dos ubicaciones está
entre 1 y 2. Para el caso del cinturon del fuego del pacifico es donde
tenemos actividad sismica mas dispersa en el sentido de la ubicación,
registrando magnitudes mayores a 4 en la escala de Rither.
earthquake %>%
group_by(location) %>%
filter(!(is.na(mag)))%>%
summarise(Number_of_quakes = length(location),
Average_Magnitude = mean(mag))%>%
mutate(Percent = round(prop.table(Number_of_quakes)*100,2))%>%
arrange(desc(Number_of_quakes))%>% top_n(25, Number_of_quakes)## # A tibble: 25 × 4
## location Number_of_quakes Average_Magnitude Percent
## <chr> <int> <dbl> <dbl>
## 1 CA 2643 1.02 32.4
## 2 Alaska 2551 1.72 31.3
## 3 Nevada 480 0.776 5.88
## 4 Utah 480 1.39 5.88
## 5 Puerto Rico 269 2.37 3.3
## 6 Hawaii 241 1.83 2.95
## 7 Montana 174 1.27 2.13
## 8 Washington 114 1.06 1.4
## 9 Indonesia 102 4.60 1.25
## 10 California 86 0.980 1.05
## # … with 15 more rowsEn el siguiente grafico de dispersión vemos como es la distribución de los terremotos dependiendo de la profundidad a la que ocurre.
ggplot(earthquake, aes(mag,depth,color = hour))+
geom_jitter(alpha = 0.5)+
gghighlight(max(depth>200)| max(mag>4))+
scale_x_continuous(breaks=seq(0,9,1))+
scale_color_viridis(direction = -1)+
theme_bw()+
xlab('Magnitude')+
ylab('Depth')+
ggtitle('Magnitude Vs. Depth in Km')+
theme(plot.title = element_text(hjust = 0.5))
Vemos como los terremotos de magnitud enre 4 y 5 ocurren a una
profundidad mayor a los 200Km, y los terremotos de magnitus superiores a
5 son relativamente dispersas. Vemos tambien como las magnitudes menores
a 4 se dan en una profundidad menor a los 200Km. Hay que tener en cuenta
que la profundida es independiente a la magnitud del terremoto, la
magnitud está ligada a la cantidad de energia que libera.
earthquake %>%
group_by(type)%>%
summarise(count = length(type))%>%
mutate(Percent = prop.table(count*100))%>%
ggplot(aes(type,Percent, fill = type))+
geom_col()+theme_bw()+
xlab('Type of Seismic Activity')+
ylab('Percent')+
ggtitle('Percent of Seismic Activity')+
theme(plot.title = element_text(hjust = 0.5))
Podemos ver que en el total de los datos hay distintas categorias de
actividades que generan terremotos; en este caso el 98% de los registros
se tienen que son de terremotos naturales.
En el sigiente mapa vemos como es la distribución de las actividades
que producen terremotos en el mundo: 
Podemor ver como los terremotos naturales son los que tienen mas
presencia, y esta tiene mayor registro en la costa Oeste de los Estado
Unidos en donde tambien se tiene datos de explosiones, y Alaska, en esta
ultima se tiene registro de criosismos,o terremotos de hielo.
Que variables determinan a un sismo y a un terremoto?
Los “sismos” y los “terremotos” se refieren los dos a movimientos o vibraciones de la tierra, pero se utilizan de manera diferente según el contexto geológico o geográfico.
Un sismo, es un movimiento de menor intensidad que se produce en la corteza terrestre. Los sismos pueden ser causados por diversas fuerzas geológicas, como la actividad volcánica, la tectónica de placas, la intrusión de magma o la falla de rocas.
Por otro lado, un terremoto es un movimiento sísmico de mayor intensidad que puede tener consecuencias devastadoras en las áreas afectadas. Un terremoto ocurre cuando hay una liberación repentina de energía acumulada en las fallas geológicas subterráneas. Esta liberación de energía causa fuertes movimientos de tierra que pueden provocar deslizamientos de tierra, tsunamis, inundaciones y daños estructurales en edificios y demas.
Profundidad: Donde se da la “ruptura, dezplazamineto o choque”.
ggplot(earthquake, aes('',depth, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 100))+theme_bw()+
xlab('Type of Seismic Activity')+
ggtitle('Depth of Seismic Activity')
En el anterior grafico pudimos ver como las actividades sismicas por lo
general se dan en la superficie, a diferencia de los terremotos ocurren
entre los 10 y 20 km y alacanzan una profundidad de varios kilometros
mas.
Magnitud:
ggplot(earthquake, aes('',mag, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 8))+theme_bw()+
xlab('Type of Seismic Activity')+
ggtitle('Magnitude of Seismic Activity')
Un terremoto por lo general tiene una magnitud entre 0 y 8, pero vemos
como los sismos debidos a otras actividades antropicas rodan
aproximandamente entre el 0 y 3 en la escala de Rither
Hora del Día:
ggplot(earthquake, aes('',hour, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 24))+theme_bw()+
xlab('Tipo de actividad sismica')+
ggtitle('Hora del evento sismico')
Por lo general todas estos eventos suceden entre las 5 de la mañana y
las 12 de la noche.
#Predicción En este apartado se intenta predecir el tipo de actividad sísmica utilizando algunas características como la profundidad, la magnitud, la latitud y la longitud.
el primer paso es clasificar todas las actividades sísmicas que no son terremotos naturales en una categoría, ya que hemos visto que la mayoría de los terremotos tienen más o menos el mismo rango de profundidad, magnitud y ubicación.
#Escogemos las columnas que vamos a necesitar, en este caso son: latitud, longitud, profundidad, magnitud, tipo de actividad y hora
quake <- earthquake[,c(2,3,4,5,24,15)]
#Caracterizamos el tipo de evento sismico
quake$type <- ifelse(quake$type == 'earthquake', 'Earthquake', 'Otherquake')
#convertimos la columna "tipo" en un factor
quake$type <- factor(quake$type)
#Se comprueba y se eliminan las filas que no contengan NA
colSums(is.na(quake))## latitude longitude depth mag hour type
## 0 0 0 4 0 0quake <- quake[!is.na(quake$mag),]#Visualizamos una tabla con nuestras variables y datos
head(quake)## # A tibble: 6 × 6
## latitude longitude depth mag hour type
## <dbl> <dbl> <dbl> <dbl> <int> <fct>
## 1 33.5 -117. 2.82 0.44 4 Earthquake
## 2 61.1 -152. 106. 1.1 4 Earthquake
## 3 38.8 -123. 1.16 0.73 4 Earthquake
## 4 65.0 -149. 8.6 1.6 4 Earthquake
## 5 37.3 -118. 6.9 0.7 4 Earthquake
## 6 35.1 -97.6 5 1.8 4 EarthquakeSe utilizaron maquinas de vectores para hacer estas predicciones porque el desequilibro que hay entre los datos es que el 98% de los datos son registrados como terremotos, y el otro 2% de los datos son de las otras actividades sismicas.
EL modelo se iniciará despresiando el desequilibrio, y con la librería ROSE instalada inicialmente se tratará de arreglar el desequilibrio
table(quake$type)##
## Earthquake Otherquake
## 8005 152set.seed(1234)
indices = sample.split(quake$type, SplitRatio = 0.7)
train = quake[indices,]
test = quake[!(indices),]
head(train)## # A tibble: 6 × 6
## latitude longitude depth mag hour type
## <dbl> <dbl> <dbl> <dbl> <int> <fct>
## 1 33.5 -117. 2.82 0.44 4 Earthquake
## 2 61.1 -152. 106. 1.1 4 Earthquake
## 3 38.8 -123. 1.16 0.73 4 Earthquake
## 4 65.0 -149. 8.6 1.6 4 Earthquake
## 5 35.1 -97.6 5 1.8 4 Earthquake
## 6 37.3 -117. 6.2 1.7 3 Earthquaketable(train$type)##
## Earthquake Otherquake
## 5604 106table(test$type)##
## Earthquake Otherquake
## 2401 46set.seed(1234)
#Usando RBF Kernel
Model_RBF <- ksvm(type~ ., data = train, scale = FALSE, kernel = "rbfdot")
Eval_RBF <- predict(Model_RBF, test[,-6])
#Matrix de confusión RBF Kernel
confusionMatrix(Eval_RBF,test$type)## Confusion Matrix and Statistics
##
## Reference
## Prediction Earthquake Otherquake
## Earthquake 2401 46
## Otherquake 0 0
##
## Accuracy : 0.9812
## 95% CI : (0.975, 0.9862)
## No Information Rate : 0.9812
## P-Value [Acc > NIR] : 0.5391
##
## Kappa : 0
##
## Mcnemar's Test P-Value : 0.00000000003247
##
## Sensitivity : 1.0000
## Specificity : 0.0000
## Pos Pred Value : 0.9812
## Neg Pred Value : NaN
## Prevalence : 0.9812
## Detection Rate : 0.9812
## Detection Prevalence : 1.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : Earthquake
## en la predicción anterior se ve el desequilibrio en los datos. Los terremotos se predicen correctamente, lo que hace que la sensibilidad sea alta, sin embargo las otras actividades sismicas tambien se predicen correctamente pero estas llevan una especifidad de 0.
Tratamiento del desequilibrio:
Se utilizará la libreria ROSE para tratar el desequilibrio utilizando el muestreo.
set.seed(1234)
over <- ovun.sample(type~., data = train, method="over", N= 11640)$dataset.seed(1234)
#Using RBF Kernel
Model_RBF <- ksvm(type~ ., data = over, scale = FALSE, kernel = "rbfdot")
Eval_RBF <- predict(Model_RBF, test[,-6])
#confusion matrix - RBF Kernel
confusionMatrix(Eval_RBF,test$type)## Confusion Matrix and Statistics
##
## Reference
## Prediction Earthquake Otherquake
## Earthquake 2285 4
## Otherquake 116 42
##
## Accuracy : 0.951
## 95% CI : (0.9416, 0.9592)
## No Information Rate : 0.9812
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.3941
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.9517
## Specificity : 0.9130
## Pos Pred Value : 0.9983
## Neg Pred Value : 0.2658
## Prevalence : 0.9812
## Detection Rate : 0.9338
## Detection Prevalence : 0.9354
## Balanced Accuracy : 0.9324
##
## 'Positive' Class : Earthquake
## Finalmente podemos ver el aumento significativo de la especifidad, 42 de 46 eventos sismicos se predicen correctamente, y tenemos una especifidad del 91%.