Decarga y lectura de datos

# Se descarga el archivo csv
if(file.exists(here::here("./data/class/raw/sismos1980_2021.csv")) == FALSE){
url_data <- "http://www2.ssn.unam.mx:8080/catalogo/csv/20211007163631YPQUH.csv"
download.file(url_data, here::here("./data/class/raw/sismos1980_2021.csv")) # Se debe modificar la ruta 
}

  # Se cargan los datos y se limpian los nombres de las variables
sismos_raw <- read_csv(here::here("./data/class/raw/sismos1980_2021.csv"), skip = 4) %>%
  janitor::clean_names() 
corazon <- read_csv(here::here("./data/class/raw/corazon.csv")) %>%
  janitor::clean_names()

Pequeño webscrapping

Se extrae una tabla de wikipedia que contiene los nombres y claves de cada estado de la república mexicana.

url_claves <- "https://es.wikipedia.org/wiki/ISO_3166-2:MX"
  
r <- GET(url_claves)
  
doc <- readHTMLTable(
  doc = content(r, "text"))
  
claves <- doc[1][[1]] %>%
  select(nombre_estado = V1,
         clave = V5) %>%
  filter(nombre_estado != "Estado" & !is.na(clave))

rm(doc, r, url_claves)

Limpieza de los datos

La limpieza consistió en convertir las variables de fecha a tipo Date y añadir variables específicas de año, mes, hora. Además se extrajo la clave de cada Estado. En este último paso se encontraron tres observaciones que no estaban bien codificadas y pertenecían al Estado de Nuevo Léon.

sismos_clean <- sismos_raw %>%
  filter(!is.na(hora)) %>% # las ultimas 6 filas contienen notas
  mutate_at(vars(fecha, fecha_utc), as.Date, format = "%Y-%m-%d") %>%
  mutate(year = year(fecha),
         mes = month(fecha), # Se agregan variables de mes, año
         hour =  hour(hora), # Se agrega variable de hora
         estado = str_extract(referencia_de_localizacion, "[A-z]*$"), # identificador de estado
         estado = case_when(estado == "N" ~ "NL", # en 3 observaciones NL aparecía como N
                            TRUE ~ estado))

Ejercicio 1

1 Total de sismos detectados por mes.

Se puede apreciar que en los meses de enero, febrero, y septiembre, la cantidad de sismos supera los 20,000 dentro del período de 1980 a 2021.

sismos_clean %>%
  group_by(mes) %>%
  summarise(mes = mes,
            sismos = n()) %>%
  distinct()
Sismos detectados por mes
Mes Total de sismos
1 21195
2 20644
3 19155
4 17959
5 17807
6 18129
7 17633
8 17827
9 21497
10 18254
11 18094
12 18570
Período de 1980 a 2021.

2 El total de sismos detectados por cada estado de la república.

En este caso, primero se cruzan los datos del web scraping para agregar una variable que contenga los nombres de los estados.

 # Homogeneizar claves y nombres
nombre_estado <- claves$nombre_estado
estado <- sort(unique(sismos_clean$estado))
cbind(nombre_estado, estado)
##       nombre_estado         estado
##  [1,] "Aguascalientes"      "AGS" 
##  [2,] "Baja California"     "BC"  
##  [3,] "Baja California Sur" "BCS" 
##  [4,] "Campeche"            "CAMP"
##  [5,] "Chiapas"             "CDMX"
##  [6,] "Chihuahua"           "CHIH"
##  [7,] "Ciudad de México"    "CHIS"
##  [8,] "Coahuila"            "COAH"
##  [9,] "Colima"              "COL" 
## [10,] "Durango"             "DGO" 
## [11,] "Guanajuato"          "GRO" 
## [12,] "Guerrero"            "GTO" 
## [13,] "Hidalgo"             "HGO" 
## [14,] "Jalisco"             "JAL" 
## [15,] "México"              "MEX" 
## [16,] "Michoacán"           "MICH"
## [17,] "Morelos"             "MOR" 
## [18,] "Nayarit"             "NAY" 
## [19,] "Nuevo León"          "NL"  
## [20,] "Oaxaca"              "OAX" 
## [21,] "Puebla"              "PUE" 
## [22,] "Querétaro"           "QR"  
## [23,] "Quintana Roo"        "QRO" 
## [24,] "San Luis Potosí"     "SIN" 
## [25,] "Sinaloa"             "SLP" 
## [26,] "Sonora"              "SON" 
## [27,] "Tabasco"             "TAB" 
## [28,] "Tamaulipas"          "TAMS"
## [29,] "Tlaxcala"            "TLAX"
## [30,] "Veracruz"            "VER" 
## [31,] "Yucatán"             "YUC" 
## [32,] "Zacatecas"           "ZAC"
  # No concuerda CHIS - CDMX, GRO - GTO, QR - QRO, SIN - SLP
  claves_mal <- c("CHIS", "CDMX", "GRO", "GTO", "QR", "QRO", "SIN", "SLP")
  index_mal <- which(estado %in% claves_mal)
  
  # Hay que reemplazar los datos
estado <- replace(estado, index_mal, estado[c(7, 5, 11, 12, 23, 22, 25, 24)])
nombre_clave <- as.data.frame(cbind(nombre_estado, estado))

  # Se agregan los nombres completos de los estados al df
sismos_clean <- merge(sismos_clean, nombre_clave, by = "estado")

rm(claves_mal, estado, index_mal, nombre_estado)

Ahora con la variable de nombres agregada, es posible generar una tabla más entendible.

Chiapas y Oaxaca son las entidades que lideran esta estadística.

sismos_clean %>%
  group_by(nombre_estado) %>%
  summarise(nombre_estado = nombre_estado,
            total_sismos = n()) %>%
  distinct()
Sismos detectados por cada estado
Estado Total sismos
Aguascalientes 42
Baja California 8170
Baja California Sur 2846
Campeche 62
Chiapas 42913
Chihuahua 400
Ciudad de México 285
Coahuila 109
Colima 5288
Durango 61
Guanajuato 39249
Guerrero 82
Hidalgo 548
Jalisco 6832
México 583
Michoacán 13979
Morelos 233
Nayarit 263
Nuevo León 423
Oaxaca 94592
Puebla 1130
Querétaro 15
Quintana Roo 86
San Luis Potosí 189
Sinaloa 920
Sonora 1309
Tabasco 530
Tamaulipas 115
Tlaxcala 120
Veracruz 5159
Yucatán 3
Zacatecas 228
Período de 1980 a 2021.

3 Gráfco circular de frecuencia de los sismos de acuerdo a la hora del día (0-23 ó 1-24) cuando acontenció.

Como podemos observar, la mayoría de los sismos dentro del período de 1980 a 2021 han ocurrido durante la noche o la madrugada.

animate(sismos_clean %>%
          group_by(hour) %>%
          summarise(hour = hour,
                    total_sismos = n()) %>%
          distinct() %>%
          ggplot(aes(hour, total_sismos, color = total_sismos, fill = total_sismos)) +
          coord_polar(theta = "x", start = -.13) +
          geom_bar(stat = "identity", width = .9) +
          geom_hline(yintercept = seq(0, 11000, by = 2000), color = "grey80", size = 0.3) +
          scale_x_continuous(breaks = 0:24, expand = c(.002,0)) +
          scale_color_gradient(low = "#F0CB35", high = "#dc2430", name = "Número de sismos\nen el período:\n{frame_time}") +
          scale_fill_gradient(low = "#F0CB35", high = "#dc2430", name = "Número de sismos\nen el período:\n{frame_time}") +
          labs(x = "", y = "Total de sismos", title = "Total de sismos por hora en México",
               subtitle = "Enero de 1980 a Octubre de 2021",
               caption = "Fuente: Elaboración propia con datos\nobtenidos del Servicio Sismológico Nacional.") +
          theme_bw() + 
          theme(text = element_text(colour = "black"),
                axis.text = element_text(colour = "black", size = 12),
                axis.text.y = element_text(size = 12),
                axis.title.y = element_text(size = 14),
                plot.title = element_text(face = "bold", colour = "black", size = 18), 
                plot.subtitle = element_text(size = 16),
                plot.caption = element_text(size = 12),
                legend.title = element_text(size = 12)) +
          transition_time(cumsum(total_sismos)) +
  shadow_mark(past = T, future = F, alpha = 0.5), nframes = 120, height = 700, width = 700, type = "cairo")

4. Considerando solamente aquellos sismos cuya magnitud fue mayor o igual a 5 ¿cuáles son los 5 estados con mayor cantidad de sismos?

sismos_clean %>%
  filter(magnitud >= 5 & magnitud != "no calculable") %>%
  group_by(nombre_estado) %>%
  summarise(nombre_estado = nombre_estado,
            total_sismos = n()) %>%
  distinct() %>%
  ungroup() %>%
  slice_max(order_by = total_sismos, n = 5)
Los cinco estados con mayor cantidad de sismos
nombre_estado total_sismos
Chiapas 499
Oaxaca 272
Guanajuato 185
Baja California Sur 115
Jalisco 105
Sismos cuya magnitud fue mayor o igual a 5.

5. Histograma de la magnitud de los sísmo detectados en el estado de Chiapas.

La mayoría de los sismos se concentran entre los 3 y 4 grados de magnitud.

plot(sismos_clean %>%
  filter(estado == "CHIS" & magnitud != "no calculable") %>%
  select(magnitud, nombre_estado) %>%
  ggplot(aes(as.numeric(magnitud), fill = magnitud, color = magnitud)) +
  geom_histogram(bins = 100, alpha = 0.5) +
  labs(x = "Magnitud",
       y = "Número de sismos",
       title = "Sismos en el Estado de Chiapas por magnitud",
       subtitle = "De Enero 1980 a Septiembre de 2021",
       caption = "Fuente: Elaboración propia con datos obtenidos del Servicio Sismológico Nacional.") +
  scale_y_continuous(breaks = seq(0, 6000, 500)) +
  scale_x_discrete(limits = as.character(c(0:10))) +
  theme_minimal() +
  theme(text = element_text(colour = "black"),
        axis.text = element_text(colour = "black", size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        plot.title = element_text(face = "bold", colour = "black", size = 14), 
        plot.subtitle = element_text(size = 12),
        plot.caption = element_text(size = 10),
        legend.position = "none"))
Magnitud discreta

Magnitud discreta

Lo anterior se puede apreciar mejor en el siguiente histograma clasificado en intervalos más concretos.

plot(sismos_clean %>%
  filter(estado == "CHIS" & magnitud != "no calculable") %>%
  select(magnitud, nombre_estado) %>%
  mutate(mag_fac = factor(case_when(magnitud >= 1 & magnitud < 2.0 ~ "Menor a 2",
                                    magnitud >= 2 & magnitud < 3.1 ~ "2 a 3",
                                    magnitud >= 3.1 & magnitud < 4.1 ~ "3.1 a 4",
                                    magnitud >= 4.1 & magnitud < 5.1 ~ "4.1 a 5",
                                    magnitud >= 5.1 & magnitud < 6.1 ~ "5.1 a 6",
                                    magnitud >= 6.1 & magnitud < 7.1 ~ "6.1 a 7",
                                    magnitud >= 7.1 & magnitud < 8.1 ~ "7.1 a 8",
                                    magnitud >= 8.1 & magnitud < 9.1 ~ "8.1 a 9",
                                    magnitud >= 9.1 & magnitud < 10.0 ~ "9.1 a 10"), 
                          levels = c("Menor a 2", 
                                     "2 a 3", 
                                     "3.1 a 4", 
                                     "4.1 a 5", 
                                     "5.1 a 6", 
                                     "6.1 a 7", 
                                     "7.1 a 8", 
                                     "8.1 a 9", 
                                     "9.1 a 10"))) %>%
  ggplot(aes(mag_fac, fill = mag_fac, color = mag_fac)) +
  geom_histogram(stat = "count", alpha = 0.5) +
  scale_y_continuous(breaks = seq(0, 30000, 3000)) +
  scale_fill_discrete() +
  labs(x = "Magnitud",
       y = "Número de sismos",
       title = "Sismos en el Estado de Chiapas por magnitud",
       subtitle = "De Enero 1980 a Septiembre de 2021",
       caption = "Fuente: Elaboración propia con datos obtenidos del Servicio Sismológico Nacional.") +
  theme_minimal() +
  theme(text = element_text(colour = "black"),
        axis.text = element_text(colour = "black", size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        plot.title = element_text(face = "bold", colour = "black", size = 14), 
        plot.subtitle = element_text(size = 12),
        plot.caption = element_text(size = 10),
        legend.position = "none") )
Magnitud por categorías

Magnitud por categorías

Ejercicio 2

Obtener el vector de estimaciones puntuales para los parámetros \(b\) cuando \(X\) = matriz de variables \(biking\ y\ smoking\) (498 x 2) y, \(y\) = \(heart\ disease\) (498 x 1).

Utilizando la fórmula para obtener el vector de estimaciones puntuales:

\[b = (X^TX)^{-1}X^Ty\]

y_heart_disease <- as.matrix(corazon$heart_disease) # y
biking_smoking <- as.matrix(corazon[c("biking", "smoking")]) # X
t_biking_smoking <- t(biking_smoking) # matriz traspuesta X^t

# vector de estimaciones puntuales
b <- (solve(t_biking_smoking%*%biking_smoking))%*%(t_biking_smoking%*%y_heart_disease)

Los valores del vector cobran sentido, ya que sugieren que por cada unidad de ciclismo el riesgo de padecer una enfermedad del corazón se reduce en 0.038 unidades, mientras que cada unidad de hábito de fumar lo incrementa en 0.62 unidades.

Vector de estimaciones puntuales
Parámetros b
biking -0.0382738
smoking 0.6230371