1 . CARGA DE LIBRERÍAS Y DATOS

En esta sección inicial, preparamos el entorno de trabajo cargando los paquetes estadísticos y de formato, y procedemos a importar y limpiar nuestra base de datos para aislar la variable cuantitativa discreta de estudio: Número de erupciones

# 1.1 Carga de librerías necesarias
library(tidyverse) 
library(gt)        
library(e1071)     

# 1.2 Importar el dataset 
datos_volcanes <- read.csv("global_volcano_eruption_intelligence (1).csv", sep = ";")

# 1.3 Limpieza a prueba de balas
datos_limpios_pluma <- datos_volcanes

# Convertimos a número. 
datos_limpios_pluma$pluma <- suppressWarnings(as.numeric(datos_limpios_pluma$volcano_eruption_count_in_dataset))

# Filtramos quitando todos los NA 
datos_limpios_pluma <- datos_limpios_pluma[!is.na(datos_limpios_pluma$pluma), ]

# Aislamos la variable para los cálculos matemáticos
var_discreta <- datos_limpios_pluma$pluma

2 . DEFINICIÓN DE INTERVALOS REPRESENTATIVOS

En esta sección calculamos los parámetros estadísticos base utilizando la Regla de Sturges. Posteriormente, determinamos el rango y la amplitud. Debido a la asimetría de los datos, ajustamos el número de intervalos a \(K = 9\) para obtener una distribución representativa con numeros enteros.

# 2.1 Parámetros estadísticos base
n_var <- length(var_discreta)
min_var <- min(var_discreta)
max_var <- max(var_discreta)

# Rango
R_var <- max_var - min_var

# 2.2 Regla de Sturges (Teórica)
k_sturges <- round(1 + 3.322 * log10(n_var))
amplitud_sturges <- R_var / k_sturges

# 2.3 Ajuste de Intervalos Representativos
k_red <- 9  
amplitud_red <- R_var / k_red

# 2.4 Cálculo de límites matemáticos exactos (Usando K=6)
limites_red <- seq(min_var, max_var, length.out = k_red + 1)

# 2.5 Redondeo exclusivo para las etiquetas visuales de la tabla
limites_rd_red <- round(limites_red, 2)
Intervalo_txt_red <- paste0("[", limites_rd_red[-length(limites_rd_red)], " - ", limites_rd_red[-1], ")")
Intervalo_txt_red[k_red] <- paste0("[", limites_rd_red[k_red], " - ", limites_rd_red[length(limites_rd_red)], "]")

3 . CÁLCULO DE FRECUENCIAS Y TABLAS DE DISTRIBUCIÓN

En esta fase procedemos a realizar el conteo de los datos. Para evidenciar el comportamiento de la distribución, primero se calculan las frecuencias y utilizando la Regla de Sturges original, lo cual nos permitirá identificar la presencia de clases nulas. Posteriormente, se presenta la tabla optimizada.

3.1 Distribución original de datos

# 2.1 Construcción de la tabla de frecuencias puntuales
tdf_discreta <- data.frame(n_eru = var_discreta) %>%
  count(n_eru, name = "ni") %>%
  arrange(n_eru) %>%
  mutate(
    hi = ni / sum(ni),
    Ni_asc = cumsum(ni),
    Hi_asc = cumsum(hi),
    # Frecuencia descendente
    Ni_desc = sum(ni) - lag(cumsum(ni), default = 0),
    Hi_desc = sum(hi) - lag(cumsum(hi), default = 0)
  )

# 2.2 Añadir fila de totales
tdf_final_disc <- tdf_discreta %>%
  add_row(
    n_eru = NA, 
    ni = sum(tdf_discreta$ni), 
    hi = 1.0, 
    Ni_asc = sum(tdf_discreta$ni), 
    Hi_asc = 1.0, 
    Ni_desc = sum(tdf_discreta$ni), 
    Hi_desc = 1.0
  )

# 2.3 Visualización con formato GT (Estilo APA)
tdf_final_disc %>%
  gt() %>%
  cols_label(
    n_eru = html("N° Erupciones<br>(X)"),
    ni = html("n<sub>i</sub>"),
    hi = html("h<sub>i</sub>"),
    Ni_asc = html("N<sub>i</sub> &uarr;"),
    Hi_asc = html("H<sub>i</sub> &uarr;"),
    Ni_desc = html("N<sub>i</sub> &darr;"),
    Hi_desc = html("H<sub>i</sub> &darr;")
  ) %>%
  tab_header(
    title = md("**Tabla N° 1**"),
    subtitle = "Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico"
  ) %>%
  fmt_number(columns = c(hi, Hi_asc, Hi_desc), decimals = 4) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.font.names = "Times New Roman",
    table_body.border.bottom.color = "black"
  )

N° Erupciones (X)	n_i	h_i	N_i ↑	H_i ↑	N_i ↓	H_i ↓
Tabla N° 1
Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico
1	129	0.1437	129	0.1437	898	1.0000
2	96	0.1069	225	0.2506	769	0.8563
3	72	0.0802	297	0.3307	673	0.7494
4	88	0.0980	385	0.4287	601	0.6693
5	70	0.0780	455	0.5067	513	0.5713
6	54	0.0601	509	0.5668	443	0.4933
7	14	0.0156	523	0.5824	389	0.4332
8	32	0.0356	555	0.6180	375	0.4176
9	36	0.0401	591	0.6581	343	0.3820
10	30	0.0334	621	0.6915	307	0.3419
11	33	0.0367	654	0.7283	277	0.3085
15	45	0.0501	699	0.7784	244	0.2717
16	48	0.0535	747	0.8318	199	0.2216
18	54	0.0601	801	0.8920	151	0.1682
19	19	0.0212	820	0.9131	97	0.1080
25	50	0.0557	870	0.9688	78	0.0869
28	28	0.0312	898	1.0000	28	0.0312
NA	898	1.0000	898	1.0000	898	1.0000

3.2 Distribución Original con Regla de Sturges

# 1. Límites para Sturges
limites_st <- seq(min_var, max_var, length.out = k_sturges + 1)
limites_rd_st <- round(limites_st, 0)

Intervalo_txt_st <- paste0("[", limites_rd_st[-length(limites_rd_st)], " - ", limites_rd_st[-1], ")")
Intervalo_txt_st[k_sturges] <- paste0("[", limites_rd_st[k_sturges], " - ", limites_rd_st[length(limites_rd_st)], "]")

# 2. Conteo de frecuencias (ni) y relativas (hi)
ni_st <- numeric(k_sturges)
for (i in 1:k_sturges) {
  if (i < k_sturges) {
    ni_st[i] <- sum(var_discreta >= limites_st[i] & var_discreta < limites_st[i+1])
  } else {
    ni_st[i] <- sum(var_discreta >= limites_st[i] & var_discreta <= max_var)
  }
}

hi_st <- round((ni_st / sum(ni_st)) * 100, 2)
hi_st[k_sturges] <- 100 - sum(hi_st[1:(k_sturges-1)])

# 3. Frecuencias acumuladas
Ni_asc_st <- cumsum(ni_st)
Ni_dsc_st <- rev(cumsum(rev(ni_st)))
Hi_asc_st <- cumsum(hi_st)
Hi_dsc_st <- rev(cumsum(rev(hi_st)))

# 4. Construcción y presentación de la Tabla
TDF_sturges <- data.frame(Intervalo = as.character(Intervalo_txt_st), ni = ni_st, hi = hi_st, Ni_dsc = Ni_dsc_st, Hi_dsc = Hi_dsc_st, Ni_asc = Ni_asc_st, Hi_asc = Hi_asc_st)
TDF_final_sturges <- rbind(TDF_sturges, data.frame(Intervalo = "Total", ni = sum(ni_st), hi = 100, Ni_dsc = NA, Hi_dsc = NA, Ni_asc = NA, Hi_asc = NA))

TDF_final_sturges %>%
  mutate(across(c(ni, hi, Ni_dsc, Hi_dsc, Ni_asc, Hi_asc), as.numeric)) %>%
  gt() %>%
  cols_label(
    Intervalo = html("Intervalo<br>de clase (km)"),
    ni = html("Frecuencia<br>absoluta<br>n<sub>i</sub>"),
    hi = html("Frecuencia<br>relativa<br>h<sub>i%</sub>"),
    Ni_dsc = html("N<sub>i</sub> &darr;"), Hi_dsc = html("H<sub>i%</sub> &darr;"),
    Ni_asc = html("N<sub>i</sub> &uarr;"), Hi_asc = html("H<sub>i%</sub> &uarr;")
  ) %>%
  tab_header(title = md("**Tabla N° 2**"), subtitle = paste0("Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico. Según Sturges (K = ", k_sturges, ")")) %>%
  tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
  cols_align(align = "center", columns = everything()) %>%
  sub_missing(missing_text = "") %>%
  fmt_number(columns = c(hi, Hi_dsc, Hi_asc), decimals = 2) %>%
  tab_options(table.font.names = "Times New Roman", column_labels.border.top.color = "black", column_labels.border.top.width = px(2), column_labels.border.bottom.color = "black", column_labels.border.bottom.width = px(1), table_body.border.bottom.color = "black", table_body.border.bottom.width = px(2))

Intervalo de clase (km)	Frecuencia absoluta n_i	Frecuencia relativa h_i%	N_i ↓	H_i% ↓	N_i ↑	H_i% ↑
Tabla N° 2
Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico. Según Sturges (K = 11)
[1 - 3)	297	33.07	898	100.00	297	33.07
[3 - 6)	158	17.59	601	66.93	455	50.66
[6 - 8)	100	11.14	443	49.34	555	61.80
[8 - 11)	66	7.35	343	38.20	621	69.15
[11 - 13)	33	3.67	277	30.85	654	72.82
[13 - 16)	45	5.01	244	27.18	699	77.83
[16 - 18)	102	11.36	199	22.17	801	89.19
[18 - 21)	19	2.12	97	10.81	820	91.31
[21 - 23)	0	0.00	78	8.69	820	91.31
[23 - 26)	50	5.57	78	8.69	870	96.88
[26 - 28]	28	3.12	28	3.12	898	100.00
Total	898	100.00
Elaborado por: Grupo 2 - Carrera de Geología.

3.3 Distribución Ajustada (Intervalos Representativos)

Como punto de partida teórico, se utilizó la Regla de Sturges, el método estadístico estándar para determinar el número de clases (\(K\)). Para una muestra de \(N = 898\) registros, esta regla sugiere matemáticamente la creación de aproximadamente 11 intervalos. Sin embargo, al aplicar este valor a la naturaleza de nuestra variable (cuyos valores oscilan entre 1 y 28 km), el cálculo de la amplitud generó fracciones decimales. Al redondear estos límites para trabajar con números enteros, se produjo una inconsistencia visual y matemática en la amplitud de las clases (algunos intervalos quedaron con un ancho de 2 km y otros de 3 km. Para solucionar la distorsión generada por el redondeo de Sturges y optimizar la presentación de la información, se aplicó un criterio técnico ajustando el número de clases a \(K = 9\). Dado que el rango exacto de los datos es de 27 (28 - 1), dividir este rango entre 9 intervalos produce una amplitud constante y exacta de 3 km para cada clase.

# 1. Conteo de frecuencias usando los límites reducidos (K=9)
ni_red <- numeric(k_red)
for (i in 1:k_red) {
  if (i < k_red) {
    ni_red[i] <- sum(var_discreta >= limites_red[i] & var_discreta < limites_red[i+1])
  } else {
    ni_red[i] <- sum(var_discreta >= limites_red[i] & var_discreta <= max_var)
  }
}

hi_red <- round((ni_red / sum(ni_red)) * 100, 2)
hi_red[k_red] <- 100 - sum(hi_red[1:(k_red-1)]) 

# 2. Frecuencias acumuladas
Ni_asc_red <- cumsum(ni_red)
Ni_dsc_red <- rev(cumsum(rev(ni_red)))
Hi_asc_red <- cumsum(hi_red)
Hi_dsc_red <- rev(cumsum(rev(hi_red)))

# 3. Construcción y presentación de la Tabla Ajustada
TDF_reducida <- data.frame(Intervalo = as.character(Intervalo_txt_red), ni = ni_red, hi = hi_red, Ni_dsc = Ni_dsc_red, Hi_dsc = Hi_dsc_red, Ni_asc = Ni_asc_red, Hi_asc = Hi_asc_red)
TDF_final_reducida <- rbind(TDF_reducida, data.frame(Intervalo = "Total", ni = sum(ni_red), hi = 100, Ni_dsc = NA, Hi_dsc = NA, Ni_asc = NA, Hi_asc = NA))

TDF_final_reducida %>%
  mutate(across(c(ni, hi, Ni_dsc, Hi_dsc, Ni_asc, Hi_asc), as.numeric)) %>%
  gt() %>%
  cols_label(
    Intervalo = html("Intervalo<br>de clase (km)"),
    ni = html("Frecuencia<br>absoluta<br>n<sub>i</sub>"),
    hi = html("Frecuencia<br>relativa<br>h<sub>i%</sub>"),
    Ni_dsc = html("N<sub>i</sub> &darr;"), Hi_dsc = html("H<sub>i%</sub> &darr;"),
    Ni_asc = html("N<sub>i</sub> &uarr;"), Hi_asc = html("H<sub>i%</sub> &uarr;")
  ) %>%
  tab_header(title = md("**Tabla N° 3**"), subtitle = "Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico. (K = 9)") %>%
  tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
  cols_align(align = "center", columns = everything()) %>%
  sub_missing(missing_text = "") %>%
  fmt_number(columns = c(hi, Hi_dsc, Hi_asc), decimals = 2) %>%
  tab_options(table.font.names = "Times New Roman", column_labels.border.top.color = "black", column_labels.border.top.width = px(2), column_labels.border.bottom.color = "black", column_labels.border.bottom.width = px(1), table_body.border.bottom.color = "black", table_body.border.bottom.width = px(2))

Intervalo de clase (km)	Frecuencia absoluta n_i	Frecuencia relativa h_i%	N_i ↓	H_i% ↓	N_i ↑	H_i% ↑
Tabla N° 3
Cantidad de eurpciones volcánicas registradas a nivel global durante el registro histórico. (K = 9)
[1 - 4)	297	33.07	898	100.00	297	33.07
[4 - 7)	212	23.61	601	66.93	509	56.68
[7 - 10)	82	9.13	389	43.32	591	65.81
[10 - 13)	63	7.02	307	34.19	654	72.83
[13 - 16)	45	5.01	244	27.17	699	77.84
[16 - 19)	102	11.36	199	22.16	801	89.20
[19 - 22)	19	2.12	97	10.80	820	91.32
[22 - 25)	0	0.00	78	8.68	820	91.32
[25 - 28]	78	8.68	78	8.68	898	100.00
Total	898	100.00
Elaborado por: Grupo 2 - Carrera de Geología.

4 . REPRESENTACIÓN GRÁFICA: HISTOGRAMAS

En esta sección se presentan las distribuciones de frecuencias de forma individual y posteriormente una comparativa integrada.

4.1 Frecuencia Absoluta Local (\(n_i\))

par(mar=c(5, 5, 9, 2))
ni_graf <- ni_red
max_ni <- max(ni_graf)

bp_ni <- barplot(ni_graf,
                 col = "#FF0000",
                 border = "black",
                 space = 0,
                 las = 1,
                 ylim = c(0, max_ni * 1.15),
                 yaxt = "n",
                 main = "",
                 xlab = "Número de erupciones",
                 ylab = "Cantidad (ni)")
mtext("Gráfica N°1.Frecuencia absoluta del número de erupciones de volcanes activos a nivel mundial (local)", 
      side = 3, line = 3, cex = 0.8, font = 2)
ticks_y <- round(seq(0, max_ni, length.out = 5), 0)
axis(side = 2, at = ticks_y, labels = ticks_y, las = 1)
axis(side = 1, at = 0:length(ni_graf), labels = limites_rd_red, cex.axis = 0.8)

par(mar=c(5, 4, 4, 2) + 0.1)

4.2 Frecuencia Absoluta Global (\(n_i\))

par(mar=c(5, 5, 8, 2))
n_total_global <- 898 # Dato global de tu script original

bp_ni_glob <- barplot(ni_graf,
                      col = "#7B3F00",
                      border = "black",
                      space = 0,
                      las = 1,
                      ylim = c(0, n_total_global * 1.05),
                      yaxt = "n",
                      main = "",
                      xlab = "Número de erupciones",
                      ylab = "Cantidad (ni)")
mtext("Gráfica N°2.Frecuencia absoluta del número de erupciones de volcanes activos a nivel mundial (global)", 
      side = 3, line = 3, cex = 0.8, font = 2)
ticks_y_global <- c(0, 200, 400, 600, 800, n_total_global)
axis(side = 2, at = ticks_y_global, labels = ticks_y_global, las = 1)
axis(side = 1, at = 0:length(ni_graf), labels = limites_rd_red, cex.axis = 0.8)

# Línea de referencia global
abline(h = n_total_global, col = "blue", lty = 2, lwd = 1.5)

4.3 Frecuencia Relativa Local (\(h_i\%\))

par(mar=c(5, 5, 9, 2))
hi_graf <- hi_red
max_hi <- max(hi_graf)

bp_hi <- barplot(hi_graf,
                 col = "red",
                 border = "black",
                 space = 0,
                 las = 1,
                 ylim = c(0, max_hi * 1.15),
                 yaxt = "n",
                 main = "",
                 xlab = "Número de erupciones",
                 ylab = "Porcentaje (%)")
mtext("Gráfica N°3.Frecuencia relativa del número de erupciones de volcanes activos a nivel mundial (local)", 
      side = 3, line = 3, cex = 0.8, font = 2)
ticks_y_hi <- seq(0, max_hi * 1.15, length.out = 5)
axis(side = 2, at = ticks_y_hi, labels = round(ticks_y_hi, 2), las = 1)
axis(side = 1, at = 0:length(hi_graf), labels = limites_rd_red, cex.axis = 0.8)

4.4 Frecuencia Relativa Global (\(h_i\%\))

par(mar=c(5, 5, 9, 2))

bp_hi_glob <- barplot(hi_graf,
                      col = "#7B3F00",
                      border = "black",
                      space = 0,
                      las = 1,
                      ylim = c(0, 100),
                      yaxt = "n",
                      main = "",
                      xlab = "Número de erupciones",
                      ylab = "Porcentaje (%)")
mtext("Gráfica N°4.Frecuencia relativa del número de erupciones de volcanes activos a nivel mundial (global)", 
      side = 3, line = 3, cex = 0.7, font = 2)

ticks_hi_global <- seq(0, 100, by = 20)
axis(side = 2, at = ticks_hi_global, labels = paste0(ticks_hi_global, "%"), las = 1)
axis(side = 1, at = 0:length(hi_graf), labels = limites_rd_red, cex.axis = 0.8)

# Línea de referencia al 100%
abline(h = 100, col = "blue", lty = 2, lwd = 1.5)

# Restaurar la configuración gráfica a un solo panel para los siguientes gráficos (Caja y Ojivas)
par(mfrow = c(1, 1))

4.5 Panel Comparativo Integrado

# Configuramos el lienzo para 2 filas y 2 columnas
par(mfrow = c(2, 2), mar = c(5, 5, 6, 2))

# 1. Frecuencia Absoluta Local
bp_ni <- barplot(ni_red, space = 0, col = "#FF0000", border = "black", las = 1,
                 ylim = c(0, max(ni_red) * 1.2), yaxt = "n",
                 main = "Frecuencia Absoluta (Local)", xlab = "Altura (km)", ylab = "Cantidad (ni)", cex.names = 0.5, las = 2)
axis(side = 2, at = round(seq(0, max(ni_red), length.out = 5), 0), las = 1)
axis(side = 1, at = 0:length(ni_red), labels = limites_rd_red, cex.axis = 0.6)


# 2. Frecuencia Absoluta Global
barplot(ni_graf, space = 0, col = "#7B3F00", border = "black", las = 1,
        ylim = c(0, n_total_global * 1.05), yaxt = "n",
        main = "Frecuencia Absoluta (Global)", xlab = "Altura (km)", ylab = "Cantidad (ni)", cex.names = 0.5, las = 2)
axis(side = 2, at = c(0, 200, 400, 600, 800, n_total_global), las = 1)
axis(side = 1, at = 0:length(ni_graf), labels = limites_rd_red, cex.axis = 0.6)
abline(h = n_total_global, col = "blue", lty = 2, lwd = 1)

# 3. Frecuencia Relativa Local
bp_hi <- barplot(hi_red, space = 0, col = "red", border = "black", las = 1,
                 ylim = c(0, max(hi_red) * 1.2), yaxt = "n",
                 main = "Frecuencia Relativa (Local)", xlab = "Altura (km)", ylab = "Porcentaje (%)", cex.names = 0.5, las = 2)
axis(side = 2, at = round(seq(0, max(hi_red), length.out = 5), 2), las = 1)
axis(side = 1, at = 0:length(hi_red), labels = limites_rd_red, cex.axis = 0.6)


# 4. Frecuencia Relativa Global
barplot(hi_graf, space = 0, col = "#7B3F00", border = "black", las = 1,
        ylim = c(0, 110), yaxt = "n",
        main = "Frecuencia Relativa (Global)", xlab = "Altura (km)", ylab = "Porcentaje (%)", cex.names = 0.5, las = 2)
axis(side = 2, at = seq(0, 100, by = 25), labels = paste0(seq(0, 100, by = 25), "%"), las = 1)
axis(side = 1, at = 0:length(hi_graf), labels = limites_rd_red, cex.axis = 0.6)
abline(h = 100, col = "blue", lty = 2, lwd = 1)

# Restaurar panel
par(mfrow = c(1, 1))

5 . DIAGRAMA DE CAJA (BOXPLOT)

El diagrama de caja nos permite visualizar la dispersión central de los datos, identificando la mediana, los cuartiles y la magnitud de los valores atípicos (outliers) .

# Ajustamos márgenes para asegurar que el título largo se visualice correctamente
par(mar = c(5, 5, 6, 2))

# Generamos el Boxplot
boxplot(var_discreta, 
        col = "lightblue", 
        border = "black",
        horizontal = TRUE, 
        xlab = "Altura de la pluma (km)",
        main = "", 
        notch = TRUE) 

# Título personalizado y centrado
title(main = "Gráfica N°5 Distribución del nnúmero de erupciones en volcanes activos con detección de valores atípicos", 
      line = 3, cex.main = 0.9)

# Agregamos la leyenda explicativa
legend("topright", 
       legend = c("Caja: 50% central (Q1 - Q3)", 
                  "Línea gruesa: Mediana", 
                  "Círculos: Valores Atípicos (Outliers)"), 
       fill = c("lightblue", NA, NA),
       border = c("black", NA, NA),
       lty = c(NA, 1, NA),
       lwd = c(NA, 2, NA),
       pch = c(NA, NA, 1),
       bty = "n", 
       cex = 0.7)

6 . GRÁFICA DE OJIVAS

En esta sección construimos las ojivas para observar el comportamiento acumulativo del número de erupciones tanto en sentido ascendente como descendente.

6.1 Ojiva de frecuencia absoluta

# Configuración de márgenes para dar espacio a la leyenda a la derecha
par(mar = c(5, 5, 4, 10), xpd = TRUE)

# Coordenadas exactas usando los límites y frecuencias de la tabla reducida
x_asc <- limites_red[-1]                # Ls (Límites superiores)
x_desc <- limites_red[-length(limites_red)] # Li (Límites inferiores)
y_asc <- Ni_asc_red                     # Acumulada ascendente
y_desc <- Ni_dsc_red                    # Acumulada descendente

# 1. Dibujar la Ascendente 
plot(x_asc, y_asc,
     type = "b", 
     main = "",
     xlab = "Número de erupciones",
     ylab = "Frecuencia absoluta acumulada (ni)",
     col = "blue",
     
     pch = 19, 
     xlim = c(min(limites_red), max(limites_red)), 
     ylim = c(0, sum(ni_red)),
     bty = "l")

# 2. Agregar la Descendente 
lines(x_desc, y_desc, col = "red", type = "b", pch = 19)

# Cuadrícula para facilitar la lectura del punto de cruce
grid()

# Título ajustado con mtext
mtext("Gráfica N°6 Ojivas Ascendentes y Descendentes de la\nDistribución de la Altura de la Pluma Volcánica", 
      side = 3, 
      line = 2, 
      adj = 0.5, 
      cex = 0.9, 
      font = 2)

# Leyenda lateral
legend("right", 
       legend = c("Ascendente", "Descendente"), 
       col = c("red", "blue"), 
       lty = 1, 
       pch = 19, 
       cex = 0.8, 
       inset = c(-0.35, 0), # Ajuste para que quede fuera del área de trazado
       bty = "n")

6.2 Ojivas Relativas Acumuladas (\(H_i\%\))

# Configuración de márgenes para la leyenda externa
par(mar = c(5, 5, 4, 10), xpd = TRUE)

# Coordenadas usando los límites y las frecuencias relativas acumuladas
x_asc <- limites_red
y_asc <- c(0, Hi_asc_red)      # Empezamos en 0 para la ascendente
x_desc <- limites_red
y_desc <- c(Hi_dsc_red, 0)     # Terminamos en 0 para la descendente

# 1. Dibujar la Ascendente Relativa
plot(x_asc, y_asc,
     type = "b", 
     main = "",
     xlab = "Número de erupciones",
     ylab = "Frecuencia relativa acumulada (%)",
     col = "blue",
     pch = 19, 
     xlim = c(min(limites_red), max(limites_red)), 
     ylim = c(0, 100),
     bty = "l")

# 2. Agregar la Descendente Relativa
lines(x_asc, y_desc, col = "red", type = "b", pch = 19)

# Cuadrícula para localizar la mediana (punto de cruce)
grid()

# Título ajustado con mtext
mtext("Gráfica N°7: Ojivas de Frecuencia Relativa Acumulada\nde la Altura de la Pluma Volcánica", 
      side = 3, 
      line = 2, 
      adj = 0.5, 
      cex = 0.9, 
      font = 2)

# Leyenda lateral
legend("right", 
       legend = c("Ascendente (%)", "Descendente (%)"), 
       col = c("blue", "red"), 
       lty = 1, 
       pch = 19, 
       cex = 0.8, 
       inset = c(-0.35, 0), 
       bty = "n")

# Restaurar márgenes
par(mar = c(5, 4, 4, 2) + 0.1, xpd = FALSE)

7 . INDICADORES ESTADÍSTICOS DE LA ALTURA DE LA PLUMA

A continuación, se presentan los indicadores estadísticos calculados para el número de erupciones, permitiendo describir el comportamiento central y la variabilidad de la muestra.

7.1 Indicadores de Tendencia Central

# 7.1 Cálculo de indicadores
media_val <- mean(var_discreta)
mediana_val <- median(var_discreta)
# Función para la moda en datos continuos (intervalos)
get_mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}
moda_val <- get_mode(var_discreta)
min_val <- min(var_discreta)
max_val <- max(var_discreta)
rango_val <- max_val - min_val

# 7.2 Creación del Data Frame
TDF_indicadores <- data.frame(
  Indicador = c("Mínimo", "Media", "Mediana", "Moda", "Máximo", "Rango"),
  Valor = c(min_val, media_val, mediana_val, moda_val, max_val, rango_val),
  Unidad = rep("km", 6),
  Interpretación = c(
    "Umbral base; todos los volcanes registran al menos un evento eruptivo.",
    "Promedio afectado por valores extremos, denotando un fuerte sesgo positivo.",
    "Punto central: el 50% de los volcanes registran 5 o menos erupciones.",
    "El nivel de actividad eruptiva más común y recurrente en el dataset.",
    "Máxima recurrencia registrada por un único sistema volcánico hiperactivo.",
    "Amplitud total que evidencia la alta heterogeneidad del dinamismo global."
  )
)

# 7.3 Generación de la Tabla con formato gt
TDF_indicadores %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 4**"),
    subtitle = "Indicadores de tendencia central de la altura de la pluma volcánica"
  ) %>%
  cols_label(
    Indicador = "Indicador",
    Valor = "Valor",
    Unidad = "Unidad",
    Interpretación = "Interpretación"
  ) %>%
  fmt_number(columns = Valor, decimals = 2) %>%
  tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
  tab_options(table.font.names = "Times New Roman")

Indicador	Valor	Unidad	Interpretación
Tabla N° 4
Indicadores de tendencia central de la altura de la pluma volcánica
Mínimo	1.00	km	Umbral base; todos los volcanes registran al menos un evento eruptivo.
Media	8.59	km	Promedio afectado por valores extremos, denotando un fuerte sesgo positivo.
Mediana	5.00	km	Punto central: el 50% de los volcanes registran 5 o menos erupciones.
Moda	1.00	km	El nivel de actividad eruptiva más común y recurrente en el dataset.
Máximo	28.00	km	Máxima recurrencia registrada por un único sistema volcánico hiperactivo.
Rango	27.00	km	Amplitud total que evidencia la alta heterogeneidad del dinamismo global.
Elaborado por: Grupo 2 - Carrera de Geología.

7.2 Indicadores de Dispersión

# 1. Cálculos de los indicadores de dispersión
varianza_val <- var(var_discreta)
desviacion_val <- sd(var_discreta)
cv_val <- (desviacion_val / media_val) * 100

# 2. Creación del Data Frame con unidades y análisis vulcanológico
TDF_dispersion <- data.frame(
  Indicador = c("Varianza", "Desviación estándar", "Coeficiente de variación"),
  Valor = c(varianza_val, desviacion_val, cv_val),
  Unidad = c("erupciones²", "erupciones", "%"),
  Interpretación = c(
    "Magnifica la dispersión; un valor alto refleja el peso de los volcanes con erupciones extremas.",
    "Muestra la fluctuación típica y la marcada diferencia en el dinamismo entre conductos volcánicos.",
    "Porcentaje elevado que define una muestra altamente heterogénea y sin un comportamiento estándar."
  )
) # <--- Si falta este paréntesis final, ocurre el error que te salió

# 3. Presentación formal con la librería gt
TDF_dispersion %>%
  mutate(Valor = ifelse(Indicador == "Coeficiente de variación", 
                        paste0(round(Valor, 2), "%"), 
                        as.character(round(Valor, 4)))) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 5**"), 
    subtitle = "Indicadores de dispersión del número de erupciones en volcanes activos"
  ) %>%
  cols_label(
    Indicador = "Indicador",
    Valor = "Valor",
    Unidad = "Unidad",
    Interpretación = "Interpretación vulcanológica"
  ) %>%
  tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
  tab_options(table.font.names = "Times New Roman",
              column_labels.border.top.color = "black", column_labels.border.top.width = px(2),
              column_labels.border.bottom.color = "black", column_labels.border.bottom.width = px(1),
              table_body.border.bottom.color = "black", table_body.border.bottom.width = px(2))

Indicador	Valor	Unidad	Interpretación vulcanológica
Tabla N° 5
Indicadores de dispersión del número de erupciones en volcanes activos
Varianza	58.6524	erupciones²	Magnifica la dispersión; un valor alto refleja el peso de los volcanes con erupciones extremas.
Desviación estándar	7.6585	erupciones	Muestra la fluctuación típica y la marcada diferencia en el dinamismo entre conductos volcánicos.
Coeficiente de variación	89.15%	%	Porcentaje elevado que define una muestra altamente heterogénea y sin un comportamiento estándar.
Elaborado por: Grupo 2 - Carrera de Geología.

7.3 Indicadores de Posición

# 1. Cálculos de posición basados en la variable discreta,
q1 <- quantile(var_discreta, 0.25)
q2 <- median(var_discreta)
q3 <- quantile(var_discreta, 0.75)
iqr <- q3 - q1
lim_inf <- q1 - 1.5 * iqr
lim_sup <- q3 + 1.5 * iqr

# 2. Detección de outliers (Uso correcto de var_discreta),
outliers <- var_discreta[var_discreta < lim_inf | var_discreta > lim_sup]
num_outliers <- length(outliers)
porc_outliers <- (num_outliers / length(var_discreta)) * 100

# 3. Construcción de la tabla con unidades e interpretaciones corregidas ,
TDF_posicion <- data.frame(
  Indicador = c("Cuartil 1 (Q1)", "Cuartil 2 (Mediana)", "Cuartil 3 (Q3)", "Rango Intercuartílico (IQR)", "Límite Inferior Outliers", "Límite Superior Outliers", "Número de Outliers"),
  Valor = c(q1, q2, q3, iqr, lim_inf, lim_sup, num_outliers),
  Unidad = c(rep("erupciones", 6), "volcanes"),
  Interpretación = c(
    "El 25% de los volcanes registran un número de erupciones menor o igual a este valor.",
    "El 50% de los volcanes presentan una frecuencia acumulada menor o igual a este nivel.",
    "El 75% de los sistemas volcánicos se sitúan por debajo de este umbral de erupciones.",
    "Amplitud o dispersión del 50% central de los volcanes evaluados.",
    "Umbral estadístico mínimo para identificar anomalías inferiores.",
    "Umbral crítico a partir del cual el número de erupciones se considera anómalo.",
    paste0("Cantidad y porcentaje de volcanes clasificados como sistemas hiperactivos (", round(porc_outliers, 2), "%).")
  )
)

# 4. Tabla final con formato gt
TDF_posicion %>%
  mutate(Valor = ifelse(Indicador == "Número de Outliers", as.character(round(Valor, 0)), as.character(round(Valor, 4)))) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 6**"), 
    subtitle = "Indicadores de posición y detección de valores atípicos en el número de erupciones"
  ) %>%
  cols_label(
    Indicador = "Indicador",
    Valor = "Valor",
    Unidad = "Unidad",
    Interpretación = "Interpretación vulcanológica"
  ) %>%
  tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
  tab_options(table.font.names = "Times New Roman",
              column_labels.border.top.color = "black", column_labels.border.top.width = px(2),
              column_labels.border.bottom.color = "black", column_labels.border.bottom.width = px(1),
              table_body.border.bottom.color = "black", table_body.border.bottom.width = px(2))

Indicador	Valor	Unidad	Interpretación vulcanológica
Tabla N° 6
Indicadores de posición y detección de valores atípicos en el número de erupciones
Cuartil 1 (Q1)	2.25	erupciones	El 25% de los volcanes registran un número de erupciones menor o igual a este valor.
Cuartil 2 (Mediana)	5	erupciones	El 50% de los volcanes presentan una frecuencia acumulada menor o igual a este nivel.
Cuartil 3 (Q3)	15	erupciones	El 75% de los sistemas volcánicos se sitúan por debajo de este umbral de erupciones.
Rango Intercuartílico (IQR)	12.75	erupciones	Amplitud o dispersión del 50% central de los volcanes evaluados.
Límite Inferior Outliers	-16.875	erupciones	Umbral estadístico mínimo para identificar anomalías inferiores.
Límite Superior Outliers	34.125	erupciones	Umbral crítico a partir del cual el número de erupciones se considera anómalo.
Número de Outliers	0	volcanes	Cantidad y porcentaje de volcanes clasificados como sistemas hiperactivos (0%).
Elaborado por: Grupo 2 - Carrera de Geología.

7.4 Indicadores de Forma

# Aseguramos el uso de la librería 'moments'
library(moments)

# 1. Cálculos basados en la variable discreta (número de erupciones)
asimetria <- skewness(var_discreta)
curtosis <- kurtosis(var_discreta) - 3 

# 2. Interpretación lógica y vulcanológica dinámica
if(abs(asimetria) < 0.5) {
    interp_asimetria <- "Distribución aproximadamente simétrica."
} else if(asimetria > 0) {
    interp_asimetria <- "Asimetría positiva (sesgo a la derecha): Predominan los volcanes con pocas erupciones, con una cola extendida de sistemas extremadamente activos."
} else {
    interp_asimetria <- "Asimetría negativa (sesgo a la izquierda)."
}

if(abs(curtosis) < 0.5) {
    interp_curtosis <- "Distribución mesocúrtica (similar a la normal)."
} else if(curtosis > 0) {
    interp_curtosis <- "Distribución leptocúrtica (muy picuda): Alta concentración de volcanes en un rango de actividad bajo, pero con presencia notable de valores extremos (outliers)."
} else {
    interp_curtosis <- "Distribución platicúrtica (más aplanada)."
}

# 3. Construcción del Data Frame
forma_df <- data.frame(
    Indicador = c("Coeficiente de Asimetría (Fisher)", "Interpretación de la Asimetría", 
                  "Coeficiente de Curtosis (Exceso)", "Interpretación de la Curtosis"),
    Valor = c(as.character(round(asimetria, 4)), interp_asimetria, as.character(round(curtosis, 4)), interp_curtosis),
    Formula = c("g₁ = E[(X-μ)³]/σ³", "|g₁|<0.5: Simétrica; g₁>0: Positiva; g₁<0: Negativa", 
                "g₂ = E[(X-μ)⁴]/σ⁴ - 3", "|g₂|<0.5: Mesocúrtica; g₂>0: Leptocúrtica; g₂<0: Platicúrtica")
)

# 4. Generación de la Tabla N° 10 con formato gt
forma_df %>%
    gt() %>%
    cols_label(
        Indicador = "Indicador",
        Valor = "Valor / Análisis",
        Formula = "Fórmula / Criterio"
    ) %>%
    cols_align(align = "center", columns = c(Indicador, Formula)) %>%
    cols_align(align = "left", columns = Valor) %>% # Se alinea a la izquierda para facilitar la lectura del texto largo
    tab_header(
        title = md("**Tabla N° 7**"),
        subtitle = "Indicadores de forma de la distribución del número de erupciones"
    ) %>%
    tab_source_note(source_note = "Elaborado por: Grupo 2 - Carrera de Geología.") %>%
    tab_options(
        table.font.names = "Times New Roman",
        column_labels.border.top.color = "black", 
        column_labels.border.top.width = px(2),
        column_labels.border.bottom.color = "black", 
        column_labels.border.bottom.width = px(1),
        table_body.border.bottom.color = "black", 
        table_body.border.bottom.width = px(2)
    )

Indicador	Valor / Análisis	Fórmula / Criterio
Tabla N° 7
Indicadores de forma de la distribución del número de erupciones
Coeficiente de Asimetría (Fisher)	1.0301	g₁ = E[(X-μ)³]/σ³
Interpretación de la Asimetría	Asimetría positiva (sesgo a la derecha): Predominan los volcanes con pocas erupciones, con una cola extendida de sistemas extremadamente activos.	\|g₁\|<0.5: Simétrica; g₁>0: Positiva; g₁<0: Negativa
Coeficiente de Curtosis (Exceso)	-0.015	g₂ = E[(X-μ)⁴]/σ⁴ - 3
Interpretación de la Curtosis	Distribución mesocúrtica (similar a la normal).	\|g₂\|<0.5: Mesocúrtica; g₂>0: Leptocúrtica; g₂<0: Platicúrtica
Elaborado por: Grupo 2 - Carrera de Geología.

Análisis Estadístico Descriptivo Variable Cuantitativa Discreta

Grupo 2 - Geología

05/06/2026