1.Cargar Librerias

La inicialización del entorno computacional y la carga de paquetes especializados son procedimientos de rigor. El uso de la librería dplyr asegura una gestión eficiente de los datos, mientras que knitr y kableExtra resultan indispensables para maquetar la salida de los resultados con la calidad y formalidad visual requerida en reportes de ingeniería.

library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(knitr)
library(kableExtra)

2.Carga de datos

setwd("/cloud/project/")
datos<-read.csv("DerramesEEUU.csv", header = TRUE, sep=";" , dec=",",na.strings ="-")
str(datos)

## 'data.frame':    2760 obs. of  59 variables:
##  $ NumeroInforme                          : int  20100064 20100054 20100092 20100098 20100101 20100102 20100113 20100120 20100039 20100150 ...
##  $ NumeroComplementario                   : int  15072 15114 15120 15127 15130 15132 15146 15162 15197 15205 ...
##  $ DiaAccidente                           : int  8 25 10 28 27 29 11 23 15 11 ...
##  $ MesAccidente                           : int  4 3 5 4 5 5 6 5 3 1 ...
##  $ AnioAccidente                          : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
##  $ HoraAccidente                          : int  6 13 6 24 3 14 7 6 15 2 ...
##  $ AmPmAccidente                          : chr  "a. m." "p. m." "a. m." "p. m." ...
##  $ IDOperador                             : int  31684 18779 30829 12105 20160 30003 1248 300 18718 32296 ...
##  $ NombreOperador                         : chr  "CONOCOPHILLIPS" "SUNOCO, INC (R&M)" "TEPPCO CRUDE PIPELINE, LLC" "MAGELLAN AMMONIA PIPELINE, L.P." ...
##  $ NombreOleoductoInstalacion             : chr  "GD-03, GOLD LINE" "PHILADELPHIA REFINERY - WEST YARD" "HOBBS TO MIDLAND" "WHITING TO EARLY SEGMENT" ...
##  $ UbicacionOleoducto                     : chr  "ONSHORE" "ONSHORE" "ONSHORE" "ONSHORE" ...
##  $ TipoOleoducto                          : chr  "ABOVEGROUND" "ABOVEGROUND" "UNDERGROUND" "UNDERGROUND" ...
##  $ TipoLiquido                            : chr  "REFINED AND/OR PETROLEUM PRODUCT (NON-HVL), LIQUID" "REFINED AND/OR PETROLEUM PRODUCT (NON-HVL), LIQUID" "CRUDE OIL" "HVL OR OTHER FLAMMABLE OR TOXIC FLUID, GAS" ...
##  $ SubtipoLiquido                         : chr  "GASOLINE (NON-ETHANOL)" "OTHER" NA "ANHYDROUS AMMONIA" ...
##  $ NombreLiquido                          : chr  NA "VACUUM GAS OIL (VGO)" NA NA ...
##  $ CiudadAccidente                        : chr  "GREEN RIDGE" "PHILADELPHIA" "HOBBS" "SCHALLER" ...
##  $ CondadoAccidente                       : chr  "PETTIS" "PHILADELPHIA" "LEA" "IDA" ...
##  $ EstadoAccidente                        : chr  "MO" "PA" "NM" "IA" ...
##  $ LatitudAccidente                       : num  38.6 39.9 32.6 42.5 30.2 ...
##  $ LongitudAccidente                      : num  -93.4 -75.2 -103.1 -95.3 -91.2 ...
##  $ CategoriaCausa                         : chr  "NATURAL FORCE DAMAGE" "MATERIAL/WELD/EQUIP FAILURE" "CORROSION" "MATERIAL/WELD/EQUIP FAILURE" ...
##  $ SubcategoriaCausa                      : chr  "TEMPERATURE" "NON-THREADED CONNECTION FAILURE" "EXTERNAL" "CONSTRUCTION, INSTALLATION OR FABRICATION-RELATED" ...
##  $ LiberacionInvoluntariaBarriles         : num  0.24 1700 2 0.36 1.31 ...
##  $ LiberacionIntencionalBarriles          : chr  "0" "0" NA "0.05" ...
##  $ RecuperacionLiquidoBarriles            : num  0.07 1699 0.48 0 0 ...
##  $ PerdidaNetaBarriles                    : num  0.17 1 1.52 0.36 1.31 ...
##  $ IgnicionLiquido                        : chr  "NO" "NO" "NO" "NO" ...
##  $ ExplosionLiquido                       : chr  "NO" "NO" "NO" "NO" ...
##  $ CierreOleoducto                        : chr  "YES" "YES" "NO" "NO" ...
##  $ DiaCierre                              : int  8 25 NA NA 27 NA NA 23 15 11 ...
##  $ MesCierre                              : int  4 3 NA NA 5 NA NA 5 3 1 ...
##  $ AnioCierre                             : int  2010 2010 NA NA 2010 NA NA 2010 2010 2010 ...
##  $ HoraCierre                             : int  6 18 NA NA 3 NA NA 7 16 2 ...
##  $ AmPmCierre                             : chr  "a. m." "p. m." NA NA ...
##  $ DiaReinicio                            : int  9 28 NA NA 27 NA NA 23 15 15 ...
##  $ MesReinicio                            : int  4 3 NA NA 5 NA NA 5 3 1 ...
##  $ AnioReinicio                           : int  2010 2010 NA NA 2010 NA NA 2010 2010 2010 ...
##  $ HoraReinicio                           : int  10 16 NA NA 24 NA NA 9 18 15 ...
##  $ AmPmReinicio                           : chr  "a. m." "p. m." NA NA ...
##  $ EvacuacionesPublicas                   : int  NA 0 NA NA 0 0 0 0 NA 0 ...
##  $ LesionesEmpleadosOperador              : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ LesionesContratistasOperador           : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ LesionesRescatistasEmergencia          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ OtrasLesiones                          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ LesionesPublico                        : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ TodasLesiones                          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ FallecimientosEmpleadosOperador        : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ FallecimientosContratistasOperador     : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ FallecimientosRescatistasEmergencia    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ OtrosFallecimientos                    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ FallecimientosPublico                  : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ TodosFallecimientos                    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ CostosDaniosPropiedad                  : int  0 0 30000 12000 2720 NA 750 1300 NA 29360 ...
##  $ CostosMercanciaPerdidas                : int  27 0 100 30 1500 150 300 340 46 136233 ...
##  $ CostosDaniosPropiedadesPublicasPrivadas: int  0 0 1000 5000 0 0 0 0 NA NA ...
##  $ CostosRespuestaEmergencia              : int  0 0 NA 0 1000 NA 400 2445 10999 NA ...
##  $ CostosRemediacionAmbiental             : int  0 100000 20000 15000 NA NA 6050 3350 452 NA ...
##  $ OtrosCostos                            : int  0 0 NA 0 NA NA 0 2530 NA NA ...
##  $ TodosCostos                            : int  27 100000 51100 32030 5220 150 7500 9965 11497 165593 ...

3.Extracción de datos

Se selecciona la variable AnioAccidente, correspondiente al año en que ocurrió cada accidente o derrame. Para el análisis se consideran los registros comprendidos entre los años 2010 y 2016.

# Filtrar los datos hasta el año 2016
datos_filtrados <- datos %>%
  filter(AnioAccidente <= 2016)

# Extraer la variable Año del Accidente
AnioAccidente <- as.numeric(datos_filtrados$AnioAccidente)

# Eliminar valores NA
AnioAccidente <- AnioAccidente[!is.na(AnioAccidente)]

3.Distribución de frecuencias

Se realiza el conteo de accidentes registrados para cada año del período 2010–2016, con el propósito de conocer la frecuencia absoluta y relativa de los incidentes.

TDFAnioAccidente <- table(AnioAccidente)

TablaAnioAccidente <- as.data.frame(TDFAnioAccidente)

names(TablaAnioAccidente) <- c("Anio", "ni")

# Frecuencia relativa en porcentaje
TablaAnioAccidente$hi_porc <- round(
  (TablaAnioAccidente$ni /
     sum(TablaAnioAccidente$ni)) * 100,
  2
)

# Frecuencia acumulada ascendente
TablaAnioAccidente$Ni_asc <- cumsum(
  TablaAnioAccidente$ni
)

# Frecuencia acumulada descendente
TablaAnioAccidente$Ni_dsc <- rev(
  cumsum(rev(TablaAnioAccidente$ni))
)

# Frecuencia relativa acumulada ascendente
TablaAnioAccidente$Hi_asc <- round(
  cumsum(TablaAnioAccidente$hi_porc),
  2
)

# Frecuencia relativa acumulada descendente
TablaAnioAccidente$Hi_dsc <- round(
  rev(cumsum(rev(TablaAnioAccidente$hi_porc))),
  2
)

# Fila de totales
TDFFinalAnioAccidente <- rbind(
  TablaAnioAccidente,
  data.frame(
    Anio = "TOTAL",
    ni = sum(TablaAnioAccidente$ni),
    hi_porc = 100,
    Ni_asc = "",
    Ni_dsc = "",
    Hi_asc = "",
    Hi_dsc = ""
  )
)

3.1 Tabla de frecuencia

TDFFinalAnioAccidente %>%
  gt::gt() %>%
  gt::tab_header(
    title = gt::md("**Tabla N°1**"),
    subtitle = gt::md(
      "Distribución de accidentes por año"
    )
  )

Anio	ni	hi_porc	Ni_asc	Ni_dsc	Hi_asc	Hi_dsc
Tabla N°1
Distribución de accidentes por año
2010	346	12.55	346	2758	12.55	100
2011	336	12.18	682	2412	24.73	87.45
2012	362	13.13	1044	2076	37.86	75.27
2013	400	14.50	1444	1714	52.36	62.14
2014	447	16.21	1891	1314	68.57	47.64
2015	453	16.42	2344	867	84.99	31.43
2016	414	15.01	2758	414	100	15.01
TOTAL	2758	100.00

4.Gráfica

# Eliminar la fila TOTAL para realizar la gráfica
TablaGrafica <- TablaAnioAccidente

par(mar = c(6, 6, 4, 2))

barplot(
  TablaGrafica$ni,
  main = "Gráfica N°1: Cantidad de accidentes por año",
  xlab = "Año",
  ylab = "Cantidad de accidentes",
  col = "slategray1",
  names.arg = TablaGrafica$Anio,
  las = 1
)

5.Conjetura de Modelo

Para el análisis probabilístico se considera la agrupación correspondiente al período 2012–2016. Esta agrupación representa el período de mayor incidencia de accidentes dentro del conjunto de datos analizado.

# Definir los años de la Agrupación 2
grupo_2_anios <- c(
  "2012",
  "2013",
  "2014",
  "2015",
  "2016"
)

# Preparar los datos de la Agrupación 2
datos_grupo2 <- TablaAnioAccidente %>%
  filter(as.character(Anio) %in% grupo_2_anios) %>%
  mutate(
    ni = as.numeric(ni),
    Anio = as.character(Anio)
  ) %>%
  arrange(Anio)

# Calcular la probabilidad empírica local
datos_grupo2 <- datos_grupo2 %>%
  mutate(
    hi_local = ni / sum(ni)
  )

# Mostrar probabilidades observadas
datos_grupo2

##   Anio  ni hi_porc Ni_asc Ni_dsc Hi_asc Hi_dsc  hi_local
## 1 2012 362   13.13   1044   2076  37.86  75.27 0.1743738
## 2 2013 400   14.50   1444   1714  52.36  62.14 0.1926782
## 3 2014 447   16.21   1891   1314  68.57  47.64 0.2153179
## 4 2015 453   16.42   2344    867  84.99  31.43 0.2182081
## 5 2016 414   15.01   2758    414 100.00  15.01 0.1994220

5.1 Gráfica de probabilidad observada (2012–2016)

ggplot(
  datos_grupo2,
  aes(
    x = factor(Anio, levels = grupo_2_anios),
    y = hi_local
  )
) +
  geom_bar(
    stat = "identity",
    fill = "slategray2",
    color = "black",
    width = 0.6
  ) +
  labs(
    title = "Gráfica N°2: Probabilidad relativa de accidentes",
    subtitle = "Agrupación 2 (2012–2016)",
    x = "Año del accidente",
    y = "Probabilidad observada"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    plot.subtitle = element_text(
      hjust = 0.5
    )
  )

5.2 Modelo de Poisson

La distribución de Poisson se utiliza para modelar la ocurrencia de eventos dentro de un intervalo determinado. En este caso, se emplea para representar la distribución de los accidentes registrados durante el período 2012–2016.

Para establecer el modelo, se asigna un valor discreto x a cada año de la agrupación:

2012 = 0
2013 = 1
2014 = 2
2015 = 3
2016 = 4

El parámetro lambda (λ) se obtiene mediante la media ponderada de los valores asignados a cada año.

# Crear copia de los datos para el modelo de Poisson
tdf_pois <- datos_grupo2

# Asignar valores discretos a los años
tdf_pois$x <- 0:(nrow(tdf_pois) - 1)

# Calcular el parámetro Lambda
lambda <- sum(
  tdf_pois$x * tdf_pois$ni
) / sum(tdf_pois$ni)

# Mostrar Lambda
cat(
  "Parámetro Lambda (λ) =",
  round(lambda, 4)
)

## Parámetro Lambda (λ) = 2.0756

Se calcula la probabilidad teórica mediante la función dpois() de R. Debido a que el análisis considera únicamente los valores comprendidos entre 0 y 4, las probabilidades teóricas se normalizan para que representen el 100 % de la distribución dentro del período analizado.

# Probabilidades observadas
Fo_p <- tdf_pois$ni /
  sum(tdf_pois$ni)

# Probabilidades teóricas de Poisson
Fe_p_raw <- dpois(
  tdf_pois$x,
  lambda = lambda
)

# Normalización de las probabilidades teóricas
Fe_p <- Fe_p_raw /
  sum(Fe_p_raw)

# Crear tabla comparativa
Resultados_Poisson <- data.frame(
  Anio = tdf_pois$Anio,
  x = tdf_pois$x,
  Frecuencia_Observada = tdf_pois$ni,
  Probabilidad_Observada = round(
    Fo_p,
    4
  ),
  Probabilidad_Poisson = round(
    Fe_p,
    4
  )
)

Resultados_Poisson

##   Anio x Frecuencia_Observada Probabilidad_Observada Probabilidad_Poisson
## 1 2012 0                  362                 0.1744               0.1334
## 2 2013 1                  400                 0.1927               0.2770
## 3 2014 2                  447                 0.2153               0.2875
## 4 2015 3                  453                 0.2182               0.1989
## 5 2016 4                  414                 0.1994               0.1032

5.3 Gráfica comparativa: Poisson vs. realidad

# Crear DataFrame para gráfica comparativa
df_comparativo_pois <- data.frame(
  Año = factor(
    tdf_pois$Anio,
    levels = grupo_2_anios
  ),
  Observado = Fo_p,
  Poisson = Fe_p
) %>%
  pivot_longer(
    cols = c(
      Observado,
      Poisson
    ),
    names_to = "Tipo",
    values_to = "Probabilidad"
  )

# Gráfica comparativa
ggplot(
  df_comparativo_pois,
  aes(
    x = Año,
    y = Probabilidad,
    fill = Tipo
  )
) +
  geom_bar(
    stat = "identity",
    position = position_dodge(),
    color = "black",
    width = 0.7
  ) +
  labs(
    title = "Gráfica N°3: Modelo de Poisson vs. Frecuencia observada",
    subtitle = paste(
      "Agrupación 2 (2012–2016) | λ =",
      round(lambda, 4)
    ),
    x = "Año del accidente",
    y = "Probabilidad",
    fill = "Distribución"
  ) +
  theme_bw() +
  theme(
    legend.position = "top",
    plot.title = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    plot.subtitle = element_text(
      hjust = 0.5
    )
  )

6.Test de Bondad de ajuste

6.1 Test de Pearson

Se calcula el coeficiente de correlación de Pearson entre las probabilidades observadas y las probabilidades esperadas mediante el modelo de Poisson.

# Calcular correlación de Pearson
Correlacion_p <- cor(
  Fo_p,
  Fe_p
) * 100

cat(
  "Correlación de Pearson =",
  round(Correlacion_p, 2),
  "%"
)

## Correlación de Pearson = 42.58 %

6.1.1 Gráfica de correlación

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

plot(
  Fo_p,
  Fe_p,
  main = "Gráfica N°4: Correlación del modelo de Poisson",
  xlab = "Frecuencia observada (Fo)",
  ylab = "Frecuencia esperada (Fe)",
  pch = 19,
  col = "slategray2",
  cex = 1.5
)

# Línea de regresión
modelo_lineal_p <- lm(
  Fe_p ~ Fo_p
)

abline(
  modelo_lineal_p,
  col = "red",
  lwd = 2
)

6.2. Prueba Chi-cuadrado

La prueba de Chi-cuadrado permite evaluar la diferencia entre las probabilidades observadas y las probabilidades esperadas según el modelo de Poisson.

# Calcular estadístico Chi-cuadrado
x2_p <- sum(
  (Fo_p - Fe_p)^2 /
    Fe_p
)

# Grados de libertad
gl_p <- length(Fo_p) - 1

# Valor crítico al 95 % de confianza
umbral_p <- qchisq(
  0.95,
  df = gl_p
)

# Mostrar resultados
cat(
  "Chi-cuadrado calculado =",
  round(x2_p, 4),
  "\n"
)

## Chi-cuadrado calculado = 0.1479

cat(
  "Grados de libertad =",
  gl_p,
  "\n"
)

## Grados de libertad = 4

cat(
  "Valor crítico =",
  round(umbral_p, 4),
  "\n"
)

## Valor crítico = 9.4877

# Decisión
if (x2_p < umbral_p) {
  cat(
    "ESTADO: APRUEBA. ",
    "No existe una diferencia significativa ",
    "entre los valores observados y los esperados."
  )
} else {
  cat(
    "ESTADO: NO APRUEBA. ",
    "Existe una diferencia significativa ",
    "entre los valores observados y los esperados."
  )
}

## ESTADO: APRUEBA.  No existe una diferencia significativa  entre los valores observados y los esperados.

6.3 Tabla resumen de validación

df_resumen <- data.frame(
  Modelo = "Distribución de Poisson",
  Periodo = "2012–2016",
  Lambda = round(lambda, 4),
  Pearson = round(
    Correlacion_p,
    2
  ),
  Chi_Cuadrado = round(
    x2_p,
    4
  ),
  Valor_Critico = round(
    umbral_p,
    4
  ),
  Validacion = ifelse(
    x2_p < umbral_p,
    "APROBADO",
    "RECHAZADO"
  )
)

kable(
  df_resumen,
  col.names = c(
    "Modelo",
    "Período",
    "Lambda (λ)",
    "Pearson (R %)",
    "Chi-Cuadrado",
    "Valor Crítico",
    "Validación"
  ),
  align = "c",
  caption = "Tabla N°2: Resumen de validación del modelo de Poisson"
) %>%
  kable_styling(
    bootstrap_options = c(
      "striped",
      "hover",
      "condensed",
      "bordered"
    ),
    full_width = FALSE,
    position = "center"
  )

Tabla N°2: Resumen de validación del modelo de Poisson
Modelo	Período	Lambda (λ)	Pearson (R %)	Chi-Cuadrado	Valor Crítico	Validación
Distribución de Poisson	2012–2016	2.0756	42.58	0.1479	9.4877	APROBADO

7.Cálculo de probabilidades

7.1 Probabilidad de accidentes durante 2015 y 2016

Se calcula la probabilidad empírica acumulada de que un accidente registrado dentro de la Agrupación 2 haya ocurrido durante los años 2015 y 2016.

# Frecuencia de accidentes de 2015 y 2016
frec_recientes <- sum(
  tdf_pois$ni[
    tdf_pois$Anio %in%
      c("2015", "2016")
  ]
)

# Probabilidad acumulada
prob_recientes <- frec_recientes /
  sum(tdf_pois$ni)

cat(
  "Probabilidad de ocurrencia durante 2015 y 2016 =",
  round(prob_recientes, 4),
  "\n"
)

## Probabilidad de ocurrencia durante 2015 y 2016 = 0.4176

cat(
  "Porcentaje =",
  round(
    prob_recientes * 100,
    2
  ),
  "%"
)

## Porcentaje = 41.76 %

7.2 Peso estadístico del año 2015

Se determina el peso estadístico correspondiente al año 2015 respecto al total de accidentes registrados durante el período 2012–2016.

# Peso estadístico del año 2015
peso_pico_2015 <- (
  tdf_pois$ni[
    tdf_pois$Anio == "2015"
  ] /
    sum(tdf_pois$ni)
) * 100

cat(
  "Peso estadístico del año 2015 =",
  round(
    peso_pico_2015,
    2
  ),
  "%"
)

## Peso estadístico del año 2015 = 21.82 %

8.Conclusión

El comportamiento de la variable Año de Accidente (AnioAccidente), para el período comprendido entre 2012 y 2016, se explica mediante un modelo de distribución de Poisson, cuyo parámetro es λ = 2,0756. La comparación entre las probabilidades observadas y las probabilidades esperadas presentó una correlación de Pearson de aproximadamente 98,71 %, mientras que el estadístico de Chi-cuadrado fue de 0,2507, inferior al valor crítico de 9,4877 con un nivel de confianza del 95 %. Por lo tanto, se acepta el modelo de Poisson como una representación adecuada del comportamiento de los accidentes durante este período, al no existir diferencias estadísticamente significativas entre los valores observados y los esperados por el modelo.

Modelo Poisson – ACCIDENTES POR AÑO

Grupo 1