1.Carga de datos

setwd("/cloud/project/")
datos<-read.csv("DerramesEEUU.csv", header = TRUE, sep=";" , dec=",")
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" "-" "ANHYDROUS AMMONIA" ...
##  $ NombreLiquido                          : chr  "-" "VACUUM GAS OIL (VGO)" "-" "-" ...
##  $ 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" "-" "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                              : chr  "8" "25" "-" "-" ...
##  $ MesCierre                              : chr  "4" "3" "-" "-" ...
##  $ AnioCierre                             : chr  "2010" "2010" "-" "-" ...
##  $ HoraCierre                             : chr  "6" "18" "-" "-" ...
##  $ AmPmCierre                             : chr  "a. m." "p. m." "-" "-" ...
##  $ DiaReinicio                            : chr  "9" "28" "-" "-" ...
##  $ MesReinicio                            : chr  "4" "3" "-" "-" ...
##  $ AnioReinicio                           : chr  "2010" "2010" "-" "-" ...
##  $ HoraReinicio                           : chr  "10" "16" "-" "-" ...
##  $ AmPmReinicio                           : chr  "a. m." "p. m." "-" "-" ...
##  $ EvacuacionesPublicas                   : chr  "-" "0" "-" "-" ...
##  $ LesionesEmpleadosOperador              : chr  "-" "-" "-" "-" ...
##  $ LesionesContratistasOperador           : chr  "-" "-" "-" "-" ...
##  $ LesionesRescatistasEmergencia          : chr  "-" "-" "-" "-" ...
##  $ OtrasLesiones                          : chr  "-" "-" "-" "-" ...
##  $ LesionesPublico                        : chr  "-" "-" "-" "-" ...
##  $ TodasLesiones                          : chr  "-" "-" "-" "-" ...
##  $ FallecimientosEmpleadosOperador        : chr  "-" "-" "-" "-" ...
##  $ FallecimientosContratistasOperador     : chr  "-" "-" "-" "-" ...
##  $ FallecimientosRescatistasEmergencia    : chr  "-" "-" "-" "-" ...
##  $ OtrosFallecimientos                    : chr  "-" "-" "-" "-" ...
##  $ FallecimientosPublico                  : chr  "-" "-" "-" "-" ...
##  $ TodosFallecimientos                    : chr  "-" "-" "-" "-" ...
##  $ CostosDaniosPropiedad                  : chr  "0" "0" "30000" "12000" ...
##  $ CostosMercanciaPerdidas                : chr  "27" "0" "100" "30" ...
##  $ CostosDaniosPropiedadesPublicasPrivadas: chr  "0" "0" "1000" "5000" ...
##  $ CostosRespuestaEmergencia              : chr  "0" "0" "-" "0" ...
##  $ CostosRemediacionAmbiental             : chr  "0" "100000" "20000" "15000" ...
##  $ OtrosCostos                            : chr  "0" "0" "-" "0" ...
##  $ TodosCostos                            : int  27 100000 51100 32030 5220 150 7500 9965 11497 165593 ...

2.Extracción de datos

AnioAccidente <- datos$AnioAccidente

# 2. Agrupacion de años manualmente
datos$grupo <- ifelse(datos$AnioAccidente %in% c(2010,2011), "2010-2011",
                      ifelse(datos$AnioAccidente %in% c(2012,2013), "2012-2013",
                             ifelse(datos$AnioAccidente %in% c(2014,2015), "2014-2015",
                                    ifelse(datos$AnioAccidente %in% c(2016,2017), "2016-2017", NA))))

2.1 Gráfica de la variable

2.1.1 Distribución de frecuencias

# Frecuencias simples
TDFAnioAccidente <- table(datos$grupo)
TablaAnioAccidente <- as.data.frame(TDFAnioAccidente)
names(TablaAnioAccidente) <- c("Anio","ni")
TablaAnioAccidente$hi <- round((TablaAnioAccidente$ni / sum(TablaAnioAccidente$ni)), 2)

2.1.2 Tabla de distribución de frecuencias

TDFFinalAnioAccidente<- rbind(TablaAnioAccidente, data.frame(
  Anio = "TOTAL",
  ni = sum(TablaAnioAccidente$ni),
  hi = 1
  ))

library(gt)
tabla_AnioAccidente <- TDFFinalAnioAccidente %>%
  gt() %>%
  cols_label(
    Anio = md("**Año**"),
    ni = md("**ni**"),
    hi = md("**hi**"),
  ) %>%
  tab_header(
    title = md("**Tabla N° 1**"),
    subtitle = md("**Distribución de accidentes en oleoductos por año en EE.UU. (2010-2017)**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 1")
  ) %>%
  tab_options(
    table.background.color = "white",
    row.striping.background_color = "white",
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    table.border.top.style = "solid",
    table.border.bottom.style = "solid",
    column_labels.font.weight = "bold",
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black"
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      rows = as.character(Anio) == "TOTAL"
    )
  )

tabla_AnioAccidente
Tabla N° 1
Distribución de accidentes en oleoductos por año en EE.UU. (2010-2017)
Año ni hi
2010-2011 682 0.25
2012-2013 762 0.28
2014-2015 900 0.33
2016-2017 416 0.15
TOTAL 2760 1.00
Autor: Grupo 1

2.1.3 Gráfica de frecuencia absoluta (Escala Local)

par(mar = c(6, 6, 4, 2)) 
barplot(
  TablaAnioAccidente$hi, 
  main = "Gráfica No.1: Distribución de probabilidad de accidentes
  por año en EE.UU.",
  xlab = "Año",
  ylab = "Probabilidad",
  col = "darkblue",
  names.arg = TablaAnioAccidente$Anio,
  las = 1,
  cex.main = 1.2,    
  cex.lab = 1.2,   
  cex.axis = 0.8,
  cex.names = 0.8
)

3. Conjetura del Modelo

Se considera que la variable AñoAccidente, podría seguir una distribución binomial. Bajo este modelo, se asume que cada año representa un “ensayo” independiente con una probabilidad constante p de que ocurra un accidente. Por lo tanto, el número total de accidentes observados en un periodo determinado es el resultado de una serie de ensayos independientes, donde el “éxito” (un accidente) tiene una probabilidad constante de ocurrir.

3.1 Definición de Hipótesis

  • Hipótesis nula (Ho): Los accidentes siguen una distribución binomial.
  • Hipótesis alternativa (H1): Los accidentes NO siguen una distribución binomial.

3.2 Parámetros del modelo Binomial

Número de categorías (años)

k <- length(TablaAnioAccidente$ni)

El Número de categorias es: 4

Asignamos valores de X: 0, 1, 2, …, a (k-1)

X <- 0:(k-1)

Media observada

media_observada <- sum(X * TablaAnioAccidente$ni) / sum(TablaAnioAccidente$ni)

Estimación de p

# Probabilidad de éxito por año
p <- media_observada / (k-1)

# Probabilidad de que no ocurra un accidente en ese “ensayo”.
q <- 1 - p

Probabilidades teóricas binomiales

P_binomial <- dbinom(X, size = k-1, prob = p)

3.3 Gráfico Observado vs Teórico

par(mar = c(6, 5, 6, 2), xpd = TRUE)
barplot(rbind(TablaAnioAccidente$hi, P_binomial),
        beside = TRUE,
        col = c("darkblue", "tomato"),
        names.arg = TablaAnioAccidente$Anio,
        main = "Modelo Binomial para AnioAccidente",
        xlab = "Año",
        ylab = "Probabilidad")

legend(x = 20, y = 0.35,
       legend = c("Observado", "Modelo"),
       fill = c("darkblue", "tomato"),
       bty = "o",
       y.intersp = 0.6,
       inset = 0.2,
       cex = 0.7)

3.3 Test

3.3.1 Test de Pearson (Correlación de frecuencias)

Fo <- TablaAnioAccidente$hi
Fe <- P_binomial

Correlacion <- cor(Fo,Fe)*100
Correlacion
## [1] 82.17718

Gráfico de correlación

plot(Fo,Fe,
     main = "Correlación Observado vs Teórico (Binomial)",
     xlab = "Teórico (Binomial)",
     ylab = "Observado",
     pch = 19, col = "blue")
abline(lm(Fe ~ Fo), col = "red", lwd = 2)
text(min(Fe), max(Fo),
     labels = paste("r =", round(Correlacion, 3)),
     pos = 4, col = "blue")

3.3.2 Test de Bondad de ajuste

Chi-cuadrado

x2 <- sum((Fo - Fe)^2 / Fe)

Estadístico Chi-cuadrado (X2): 0.1206164

Nivel de significancia y grados de libertad

nivel_significancia <- 0.05
grados_libertad <- k - 1

Valor crítico

Vc <- qchisq(1 - nivel_significancia, grados_libertad)

Valor crítico: 7.814728

3.3 Tabla resumen de test

4. Decisión del Modelo

if (x2 < Vc) {
  cat("Conclusión: No se rechaza H0, los accidentes podrían seguir una distribución binomial.")
} else {
  cat("Conclusión: Se rechaza H0, los accidentes NO siguen una distribución binomial.")
}
## Conclusión: No se rechaza H0, los accidentes podrían seguir una distribución binomial.

4.Cálculo de probabilidades

dbinom(3, size = length(X)-1, prob=p)
## [1] 0.09742803

5.Conclusion