Liberación involuntaria

1 Cargar datos

datos <- read.csv("database-_1_.csv")
zona<-datos$Unintentional.Release..Barrels.

2 Tabla de frecuencia

library(gt)
library(dplyr)

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels. 
costos <- na.omit(variable_interes)

k <- 1 + (3.322 * log10(length(costos)))
k <- floor(k)
min_val <- min(costos)
max_val <- max(costos)
R_val <- max_val - min_val
A <- R_val / k

Li_num <- seq(from = min_val, to = max_val - A, by = A)
if(length(Li_num) < k) { Li_num <- c(Li_num, Li_num[length(Li_num)] + A) }
if(max(Li_num) + A < max_val) { Li_num <- c(Li_num, tail(Li_num, 1) + A) }

Ls_num <- Li_num + A
MC_num <- (Li_num + Ls_num) / 2

ni <- numeric(length(Li_num))
for (i in 1:length(Li_num)) {
  if (i == length(Li_num)) {
      ni[i] <- sum(costos >= Li_num[i] & costos <= (max_val + 100000))
  } else {
      ni[i] <- sum(costos >= Li_num[i] & costos < Ls_num[i]) 
  }
}
num_ceros <- sum(ni == 0)
if (num_ceros > 0) {
  ni[ni == 0] <- 1
  idx_max <- which.max(ni)
  ni[idx_max] <- ni[idx_max] - num_ceros
}

hi <- ni / sum(ni) * 100
Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))
Hiasc <- round(cumsum(hi), 2)
Hidsc <- round(rev(cumsum(rev(hi))), 2)

TDFCostos <- data.frame(
  Li_num = Li_num, 
  Ls_num = Ls_num, 
  MC_num = MC_num, 
  ni = ni, 
  hi = hi, 
  Niasc = Niasc, 
  Nidsc = Nidsc, 
  Hiasc = Hiasc, 
  Hidsc = Hidsc
)

tabla1_sturges <- TDFCostos %>%
  gt() %>%
  tab_header(
    title = md("*Tabla 1: Distribución de Frecuencias*"),
  ) %>%
  cols_label(
    Li_num = "Desde (bbl)",
    Ls_num = "Hasta (bbl)",
    MC_num = "Marca Clase",
    ni = "Frec. Abs.",
    hi = "Frec. Rel. %",
    Niasc = "Ni Asc.",
    Nidsc = "Ni Desc.",
    Hiasc = "Hi Asc. %",
    Hidsc = "Hi Desc. %"
  ) %>%
  fmt_number(columns = c(Li_num, Ls_num, MC_num), decimals = 2) %>%
  fmt_number(columns = c(hi, Hiasc, Hidsc), decimals = 2, pattern = "{x}%")

tabla1_sturges

Tabla 1: Distribución de Frecuencias
Desde (bbl)	Hasta (bbl)	Marca Clase	Frec. Abs.	Frec. Rel. %	Ni Asc.	Ni Desc.	Hi Asc. %	Hi Desc. %
0.00	2,547.08	1,273.54	2743	98.14%	2743	2795	98.14%	100.00%
2,547.08	5,094.17	3,820.62	28	1.00%	2771	52	99.14%	1.86%
5,094.17	7,641.25	6,367.71	8	0.29%	2779	24	99.43%	0.86%
7,641.25	10,188.33	8,914.79	4	0.14%	2783	16	99.57%	0.57%
10,188.33	12,735.42	11,461.88	3	0.11%	2786	12	99.68%	0.43%
12,735.42	15,282.50	14,008.96	2	0.07%	2788	9	99.75%	0.32%
15,282.50	17,829.58	16,556.04	1	0.04%	2789	7	99.79%	0.25%
17,829.58	20,376.67	19,103.12	2	0.07%	2791	6	99.86%	0.21%
20,376.67	22,923.75	21,650.21	1	0.04%	2792	4	99.89%	0.14%
22,923.75	25,470.83	24,197.29	1	0.04%	2793	3	99.93%	0.11%
25,470.83	28,017.92	26,744.38	1	0.04%	2794	2	99.96%	0.07%
28,017.92	30,565.00	29,291.46	1	0.04%	2795	1	100.00%	0.04%

3 Cantidad Global de Derrames (ni)

library(ggplot2)
library(dplyr)
library(scales)
variable_interes <- datos$Unintentional.Release..Barrels.
volumen <- na.omit(variable_interes)
k <- 1 + (3.322 * log10(length(volumen)))
R <- max(volumen) - min(volumen)
A <- R / floor(k) 

limit_zoom <- A 
datos_zoom <- datos %>%
  filter(!is.na(Unintentional.Release..Barrels.)) %>%
  filter(Unintentional.Release..Barrels. <= limit_zoom)
p_zoom <- ggplot(datos_zoom, aes(x = Unintentional.Release..Barrels.)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white", alpha = 0.8) +
  scale_x_continuous(labels = scales::comma_format(suffix = " bbl")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  
  labs(
    title = "Gráfica 1: Distribución de Liberación involuntaria", 
    x = "Volumen (bbl)",
    y = "Cantidad"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

print(p_zoom)

4 Análisis del Rango Principal (Cantidad Local)

Se seleccionó el intervalo de 0 a 2.0 barriles para el análisis detallado debido a que la distribución global presenta una asimetría positiva extrema, concentrando más del 90% de los eventos en el primer intervalo del histograma general. Este enfoque permite ‘desempaquetar’ la información oculta en la barra de mayor frecuencia, revelando la verdadera estructura probabilística de los incidentes más habituales y permitiendo una visualización clara de la Moda y la Mediana que de otro modo quedarían enmascaradas por la escala de los valores atípicos.

library(ggplot2)
library(dplyr)
library(scales)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels.
volumen <- na.omit(variable_interes)
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Unintentional.Release..Barrels.)) %>%
  filter(Unintentional.Release..Barrels. <= limit_zoom)

p_zoom_5barras <- ggplot(datos_zoom, aes(x = Unintentional.Release..Barrels.)) +
  geom_histogram(bins = 5, fill = "steelblue", color = "white", alpha = 0.8) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1, suffix = " bbl")) +
  
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  
  labs(
    title = "Gráfica No2: Distribución detallada en el rango de mayor ocurrencia (0 - 2 bbl)", 
    x = "Volumen (bbl)",
    y = "Cantidad"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras)

5 Cantidad relativa local del rango principal

library(ggplot2)
library(dplyr)
library(scales)
datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels.
volumen <- na.omit(variable_interes)
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Unintentional.Release..Barrels.)) %>%
  filter(Unintentional.Release..Barrels. <= limit_zoom)

p_zoom_5barras_pct <- ggplot(datos_zoom, aes(x = Unintentional.Release..Barrels.)) +
  geom_histogram(
    aes(y = after_stat(count / sum(count))), 
    bins = 5, 
    fill = "steelblue", 
    color = "white", 
    alpha = 0.8
  ) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1, suffix = " bbl")) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1), 
    expand = expansion(mult = c(0, 0.05))
  ) +
  
  labs(
    title = "Gráfica No3 : Distribución porcentual local del rango principal", 
    x = "Volumen (bbl)",
    y = "Porcentaje (%)"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras_pct)

6 Cantidad relativa global del rango principal

library(ggplot2)
library(dplyr)
library(scales)

datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels.
limit_zoom <- 2

datos_zoom <- datos %>%
  filter(!is.na(Unintentional.Release..Barrels.)) %>%
  filter(Unintentional.Release..Barrels. <= limit_zoom)

p_zoom_5barras_100 <- ggplot(datos_zoom, aes(x = Unintentional.Release..Barrels.)) +
  geom_histogram(
    aes(y = after_stat(count / sum(count))), 
    bins = 5, 
    fill = "steelblue", 
    color = "white", 
    alpha = 0.8
  ) +
  scale_x_continuous(breaks = seq(0, 2, by = 0.4), labels = number_format(accuracy = 0.1, suffix = " bbl")) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    limits = c(0, 1), 
    expand = expansion(mult = c(0, 0.05))
  ) +
  
  labs(
    title = "Gráfica No4 : Distribución porcentual global del rango principal", 
    x = "Volumen (bbl)",
    y = "Porcentaje (%)"
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.x = element_text(hjust = 0.5)
  )

print(p_zoom_5barras_100)

7 Ojivas convinadas de la cantidad absoluta

library(ggplot2)
library(dplyr)
library(scales)
datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels.
volumen <- na.omit(variable_interes)
datos_local <- volumen[volumen <= 2]
k_local <- 6
min_val <- 0
max_val <- 2
R_local <- max_val - min_val
A_local <- R_local / k_local 
Li_num <- seq(min_val, max_val - A_local, length.out = k_local)
Ls_num <- Li_num + A_local
ni_local <- numeric(k_local)
for(i in 1:k_local){
  if(i == k_local){
    ni_local[i] <- sum(datos_local >= Li_num[i] & datos_local <= max_val)
  } else {
    ni_local[i] <- sum(datos_local >= Li_num[i] & datos_local < Ls_num[i])
  }
}

Niasc <- cumsum(ni_local)
Nidsc <- rev(cumsum(rev(ni_local)))
datos_asc <- data.frame(
  x = c(min_val, Ls_num),  
  y = c(0, Niasc),         
  Tipo = "Ascendente"
)

datos_dsc <- data.frame(
  x = c(Li_num, max_val),  
  y = c(Nidsc, 0),         
  Tipo = "Descendente"
)

datos_ojivas_plot <- rbind(datos_asc, datos_dsc)
p_ojiva_cruzada_solida <- ggplot(datos_ojivas_plot, aes(x = x, y = y, color = Tipo, linetype = Tipo)) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2) +
  
  scale_x_continuous(
    labels = scales::number_format(accuracy = 0.01, suffix = " bbl"),
    breaks = scales::pretty_breaks(n = 6)
  ) +
  scale_color_manual(values = c("Ascendente" = "black", "Descendente" = "blue")) +
  scale_linetype_manual(values = c("Ascendente" = "solid", "Descendente" = "solid")) +
  
  labs(
    title = "Gráfica 5: Distribución Acumulada Comparativa del rango principal",
    x = "Volumen (bbl)",
    y = "Cantidad Acumulada",
    color = NULL,
    linetype = NULL
  ) +
  
  theme_bw() + 
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    legend.position = c(0.85, 0.5),
    legend.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black")
  )

print(p_ojiva_cruzada_solida)

#Esta gráfica presenta las Ojivas de Frecuencia Acumulada para la variable de Volumen (medido en barriles), mostrando cómo se distribuyen los datos a lo largo del rango de 0 a 2 bbl. La línea negra (Ascendente) indica la cantidad de observaciones que están por debajo de un volumen determinado (acumulación "menor que"), mientras que la línea azul (Descendente) muestra cuántas observaciones restan por encima de ese mismo valor (acumulación "mayor que"). El punto donde ambas curvas se cruzan es visualmente importante porque representa la mediana aproximada del conjunto de datos, dividiendo la muestra en dos partes iguales.

8 Diagrama de cajas

library(ggplot2)
library(dplyr)
datos <- read.csv("database-_1_.csv")
variable_interes <- datos$Unintentional.Release..Barrels.
volumen <- na.omit(variable_interes)
variable_box <- volumen[volumen <= 2] # Filtro Local
par(mar = c(5, 2, 4, 2)) 
b <- boxplot(variable_box, 
        horizontal = TRUE, 
        col = "skyblue", 
        border = "gray30",         
        medcol = "red",            
        boxwex = 0.6,              
        outline = FALSE,           
        main = "Gráfica 6: Distribución de Volumen del rango principal",
        xlab = "Volumen (bbl)",
        xaxt = "n",  
        yaxt = "n",  
        frame = FALSE 
)

limite_visible <- b$stats[5] 
puntos_eje <- pretty(c(0, limite_visible))             
axis(1, at = puntos_eje, labels = format(puntos_eje, nsmall = 2), col = "gray30", col.axis = "gray30")
grid(nx = NULL, ny = NA, col = "lightgray", lty = "dotted", lwd = 1)

#Esta gráfica nos permite ver rápidamente dónde se concentran la mayoría de los datos y qué tan dispersos están. La caja azul representa el núcleo de la información (el 50% central de los casos), indicando que la mayoría de los volúmenes habituales oscilan entre 0.25 y 1.00 bbl.La línea roja dentro de la caja marca la mediana (aprox. 0.50 bbl), que es el punto de equilibrio exacto: la mitad de los registros son menores a esa cifra y la otra mitad mayores. Además, al ver que la línea punteada de la derecha es más larga que la de la izquierda, podemos concluir que existe una mayor dispersión en los volúmenes altos; es decir, los datos se 'estiran' más hacia los 2 barriles que hacia el 0.

9 Tabla estadísticos local

library(e1071)

## 
## Adjuntando el paquete: 'e1071'

## The following object is masked from 'package:ggplot2':
## 
##     element

library(knitr)
datos <- read.csv("database-_1_.csv")
variable_raw <- na.omit(datos$Unintentional.Release..Barrels.)
variable_analisis <- variable_raw[variable_raw <= 2]
n <- length(variable_analisis)
k_local <- 6
min_val <- 0
max_val <- 2
A_local <- (max_val - min_val) / k_local
Li <- seq(min_val, max_val - A_local, length.out = k_local)
Ls <- Li + A_local
MC <- (Li + Ls) / 2

ni <- numeric(k_local)
for(i in 1:k_local){
  if(i == k_local) ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis <= max_val)
  else ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis < Ls[i])
}

media_agrupada <- sum(MC * ni) / sum(ni)
desviacion_estandar <- sd(variable_analisis)
error_estandar <- desviacion_estandar / sqrt(n)
margen_error <- 1.96 * error_estandar
ic_inferior <- media_agrupada - margen_error
ic_superior <- media_agrupada + margen_error
texto_media_intervalo <- paste0("[", format(round(ic_inferior, 3), nsmall=3), " - ", format(round(ic_superior, 3), nsmall=3), "]")
ri <- min(variable_analisis)
rs <- max(variable_analisis)
mediana <- median(variable_analisis)
t <- table(variable_analisis)
Mo <- as.numeric(names(t)[which.max(t)])
cv <- (desviacion_estandar / media_agrupada) * 100
As <- skewness(variable_analisis)
K <- kurtosis(variable_analisis)
Tabla_local <- data.frame(
  "Volumen (Local 0-2)",
  paste(format(ri, nsmall=2), "bbl"),
  paste(format(rs, nsmall=2), "bbl"),
  texto_media_intervalo,  
  paste(format(round(mediana, 3), big.mark=","), "bbl"),
  paste(format(round(Mo, 3), big.mark=","), "bbl"),
  paste(format(round(desviacion_estandar, 3), big.mark=","), "bbl"),
  paste(round(cv, 2), "%"),
  round(As, 2),
  round(K, 2)
)

colnames(Tabla_local) <- c("Variable","Min","Max","Media (IC 95%)","Mediana","Moda","Desv. S","CV","As","K")
kable(Tabla_local, format = "markdown", caption = "Tabla 1: Indicadores Locales (con Intervalo de Confianza para la Media).")

Tabla 1: Indicadores Locales (con Intervalo de Confianza para la Media).

Variable	Min	Max	Media (IC 95%)	Mediana	Moda	Desv. S	CV	As	K
Volumen (Local 0-2)	0.00 bbl	2.00 bbl	[0.659 - 0.719]	0.5 bbl	1 bbl	0.568 bbl	82.47 %	1.04	0.15

10 Tabla estadísticos global

library(e1071)
library(knitr)
datos <- read.csv("database-_1_.csv")
variable_analisis <- na.omit(datos$Unintentional.Release..Barrels.)
n <- length(variable_analisis)
k_global <- floor(1 + 3.322 * log10(n))
R <- max(variable_analisis) - min(variable_analisis)
A_global <- R / k_global
Li <- seq(min(variable_analisis), max(variable_analisis) - A_global, length.out = k_global)
if(max(Li) + A_global < max(variable_analisis)) { 
   Li <- c(Li, tail(Li, 1) + A_global) 
   k_global <- k_global + 1
}
Ls <- Li + A_global
MC <- (Li + Ls) / 2

ni <- numeric(length(MC))
for(i in 1:length(MC)){
  if(i == length(MC)) ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis <= (max(variable_analisis) + 0.001))
  else ni[i] <- sum(variable_analisis >= Li[i] & variable_analisis < Ls[i])
}
media_agrupada <- sum(MC * ni) / sum(ni)
desviacion_estandar <- sd(variable_analisis)
error_estandar <- desviacion_estandar / sqrt(n)
margen_error <- 1.96 * error_estandar
ic_inferior <- media_agrupada - margen_error
ic_superior <- media_agrupada + margen_error
texto_media_intervalo <- paste0("[", format(round(ic_inferior, 2), big.mark=","), " - ", format(round(ic_superior, 2), big.mark=","), "]")
ri <- min(variable_analisis)
rs <- max(variable_analisis)
mediana <- median(variable_analisis)
t <- table(variable_analisis)
Mo <- as.numeric(names(t)[which.max(t)])
cv <- (desviacion_estandar / media_agrupada) * 100
As <- skewness(variable_analisis)
K <- kurtosis(variable_analisis)
Tabla_global <- data.frame(
  "Volumen (Global)",
  paste(format(ri, nsmall=2), "bbl"),
  paste(format(rs, big.mark=","), "bbl"),
  texto_media_intervalo, 
  paste(format(round(mediana, 2), big.mark=","), "bbl"),
  paste(format(round(Mo, 2), big.mark=","), "bbl"),
  paste(format(round(desviacion_estandar, 2), big.mark=","), "bbl"),
  paste(round(cv, 2), "%"),
  round(As, 2),
  round(K, 2)
)

colnames(Tabla_global) <- c("Variable","Min","Max","Media (IC 95%)","Mediana","Moda","Desv. S","CV","As","K")
kable(Tabla_global, format = "markdown", caption = "Tabla 2: Indicadores Globales (con Intervalo de Confianza para la Media).")

Tabla 2: Indicadores Globales (con Intervalo de Confianza para la Media).

Variable	Min	Max	Media (IC 95%)	Mediana	Moda	Desv. S	CV	As	K
Volumen (Global)	0.00 bbl	30,565 bbl	[1,341.27 - 1,442.75]	2 bbl	1 bbl	1,368.67 bbl	98.32 %	13.61	227.81

11 Valores Atípicos global

# --- ANÁLISIS GLOBAL ---
variable_global <- na.omit(datos$Unintentional.Release..Barrels.)
stats_outliers_global <- boxplot.stats(variable_global)$out
num_outliers_global <- length(stats_outliers_global)
minimooutliers_global <- if(num_outliers_global > 0) min(stats_outliers_global) else NA
maximooutliers_global <- if(num_outliers_global > 0) max(stats_outliers_global) else NA

cat("\n--- Análisis de Outliers (GLOBAL) ---\n")

## 
## --- Análisis de Outliers (GLOBAL) ---

cat("Número de valores atípicos:", num_outliers_global, "\n")

## Número de valores atípicos: 497

cat("Mínimo Outlier:", if(!is.na(minimooutliers_global)) paste(format(minimooutliers_global, big.mark=","), "bbl") else "Ninguno", "\n")

## Mínimo Outlier: 50 bbl

cat("Máximo Outlier:", if(!is.na(maximooutliers_global)) paste(format(maximooutliers_global, big.mark=","), "bbl") else "Ninguno", "\n")

## Máximo Outlier: 30,565 bbl

12 Valores atípicos del rango principal

# --- ANÁLISIS LOCAL (ZOOM 0-2) ---
variable_local <- variable_global[variable_global <= 2]
stats_outliers_local <- boxplot.stats(variable_local)$out
num_outliers_local <- length(stats_outliers_local)
minimooutliers_local <- if(num_outliers_local > 0) min(stats_outliers_local) else NA
maximooutliers_local <- if(num_outliers_local > 0) max(stats_outliers_local) else NA

cat("\n--- Análisis de Outliers (LOCAL 0-2 bbl) ---\n")

## 
## --- Análisis de Outliers (LOCAL 0-2 bbl) ---

cat("Número de valores atípicos:", num_outliers_local, "\n")

## Número de valores atípicos: 0

cat("Mínimo Outlier:", if(!is.na(minimooutliers_local)) paste(format(minimooutliers_local, big.mark=","), "bbl") else "Ninguno", "\n")

## Mínimo Outlier: Ninguno

cat("Máximo Outlier:", if(!is.na(maximooutliers_local)) paste(format(maximooutliers_local, big.mark=","), "bbl") else "Ninguno", "\n")

## Máximo Outlier: Ninguno

13 Conclusión

#El análisis de la variableLiberación involuntaria,  muestra que el Volumen de Derrame revela una distribución con un sesgo positivo extremo, caracterizada por una dualidad operativa clara: mientras que la inmensa mayoría de los incidentes son eventos de baja magnitud (concentrados masivamente entre 0 y 2 barriles, lo que sitúa la Moda y la Mediana en valores mínimos), el promedio (Media) se encuentra severamente inflado por la presencia de valores atípicos catastróficos. Esto concluye que el "derrame típico" es casi insignificante en volumen, pero la gestión de riesgos no debe basarse en el promedio, sino en mitigar esa pequeña probabilidad de eventos extremos que representan el mayor impacto ambiental y económico.