library(openxlsx)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gt)

1.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 ...
# 1.1 Extracción de datos
AnioAccidente <- as.numeric(datos$AnioAccidente)

# ELIMINAR SOLO 2017
AnioAccidente <- AnioAccidente[AnioAccidente != 2017]

2.Distribución de Frecuencias

TDFAnioAccidente <- table(AnioAccidente)
TablaAnioAccidente <- as.data.frame(TDFAnioAccidente)
names(TablaAnioAccidente) <- c("Anio","ni")

TablaAnioAccidente$hi_porc <- round((TablaAnioAccidente$ni / sum(TablaAnioAccidente$ni))*100,2)
TablaAnioAccidente$Ni_asc <- cumsum(TablaAnioAccidente$ni)
TablaAnioAccidente$Ni_dsc <- rev(cumsum(rev(TablaAnioAccidente$ni)))
TablaAnioAccidente$Hi_asc <- round(cumsum(TablaAnioAccidente$hi_porc),3)
TablaAnioAccidente$Hi_dsc <- round(rev(cumsum(rev(TablaAnioAccidente$hi_porc))),3)

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

Tabla

TDFFinalAnioAccidente %>%
  gt() %>%
  tab_header(
    title=md("**Tabla N°1**"),
    subtitle=md("Distribución de accidentes por año")
  )
Tabla N°1
Distribución de accidentes por año
Anio ni hi_porc Ni_asc Ni_dsc Hi_asc Hi_dsc
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

Gráfica

par(mar=c(6,6,4,2))
barplot(TablaAnioAccidente$ni,
        main="Gráfica No.1: Cantidad de accidentes por año",
        xlab="Año",
        ylab="Cantidad",
        col="slategray1",
        names.arg=TablaAnioAccidente$Anio,
        las=1)

3. PRIMER PERIODO (2010–2013) → MODELO UNIFORME

Accidentes_2010_2013 <- subset(AnioAccidente,
                               AnioAccidente>=2010 & AnioAccidente<=2013)

TDF_2010_2013 <- table(Accidentes_2010_2013)
Tabla_2010_2013 <- as.data.frame(TDF_2010_2013)
names(Tabla_2010_2013) <- c("Anio","Frecuencia")

Tabla_2010_2013$hi1 <- Tabla_2010_2013$Frecuencia/sum(Tabla_2010_2013$Frecuencia)

GRÁFICA PROBABILIDAD

barplot(Tabla_2010_2013$hi1,
        main="Gráfica N°2: Probabilidad (2010–2013)",
        xlab="Año",
        ylab="Probabilidad",
        names.arg=Tabla_2010_2013$Anio,
        col="slategray2")

3.1 MODELO UNIFORME

Fo <- Tabla_2010_2013$Frecuencia
k <- length(Fo)
total_accidentes <- sum(Fo)

Fe <- rep(total_accidentes/k,k)

barplot(rbind(Fo,Fe),
        beside=TRUE,
        col=c("slategray2","slategray4"),
        names.arg=Tabla_2010_2013$Anio,
        xlab="Año",
        ylab="Frecuencia",
        ylim=c(0,500))

title(main="Gráfica No.3: Modelo Uniforme vs Observado")

legend("topright",
       legend=c("Observado","Uniforme"),
       fill=c("slategray2","slategray4"))

3.2 TESTS MODELO UNIFORME

Pearson

Correlacion_U <- cor(Fo,Fe)*100
## Warning in cor(Fo, Fe): the standard deviation is zero
cat("\nCorrelación Uniforme =",Correlacion_U,"%\n")
## 
## Correlación Uniforme = NA %

Grafica correlacion

plot(Fo,Fe,
     main="Gráfica No.4: Correlación modelo Uniforme",
     xlab="Frecuencia Observada",
     ylab="Frecuencia Esperada",
     col="slategray2",
     pch=19)

abline(lm(Fe~Fo),col="red",lwd=2)

Chi-cuadrado

x2_u <- sum((Fo-Fe)^2/Fe)
gl_u <- k-1
umbral_u <- qchisq(0.95,gl_u)

cat("Chi2 =",x2_u,"\n")
## Chi2 = 6.570637
cat("Umbral =",umbral_u,"\n")
## Umbral = 7.814728

4. SEGUNDO PERIODO (2014–2016) → POISSON

Accidentes_2014_2017 <- subset(AnioAccidente,
                               AnioAccidente>=2014 & AnioAccidente<=2017)

TDF_2014_2017 <- table(Accidentes_2014_2017)
Tabla_2014_2017 <- as.data.frame(TDF_2014_2017)
names(Tabla_2014_2017) <- c("Anio","ni2")

Tabla_2014_2017$hi2 <- Tabla_2014_2017$ni2/sum(Tabla_2014_2017$ni2)

GRAFICA

barplot(Tabla_2014_2017$hi2,
        main="Gráfica N°5: Probabilidad (2014–2017)",
        xlab="Año",
        ylab="Probabilidad",
        names.arg=Tabla_2014_2017$Anio,
        col="slategray2")

4.1 MODELO POISSON

Tabla_2014_2017$Clase <- 1:nrow(Tabla_2014_2017)

lambda <- sum(Tabla_2014_2017$Clase*Tabla_2014_2017$ni2)/sum(Tabla_2014_2017$ni2)

Fe_p <- dpois(Tabla_2014_2017$Clase,lambda)
Fo_p <- Tabla_2014_2017$hi2

barplot(rbind(Fo_p,Fe_p),
        beside=TRUE,
        main="Gráfica No.6: Poisson vs Observado",
        xlab="Año",
        ylab="Probabilidad",
        names.arg=Tabla_2014_2017$Anio,
        col=c("slategray2","skyblue4"))

legend("topright",
       legend=c("Observado","Poisson"),
       fill=c("slategray2","skyblue4"))

4.2 TESTS POISSON

Pearson

Correlacion_p <- cor(Fo_p,Fe_p)*100
cat("\nCorrelación Poisson =",Correlacion_p,"%\n")
## 
## Correlación Poisson = 98.47304 %

Grafica correlacion

plot(Fo_p,Fe_p,
     main="Gráfica No.7: Correlación modelo Poisson",
     xlab="Frecuencia Observada",
     ylab="Frecuencia Esperada",
     col="slategray2",
     pch=19)

abline(lm(Fe_p~Fo_p),col="red",lwd=2)

Chi-cuadrado

x2_p <- sum((Fo_p-Fe_p)^2/Fe_p)
gl_p <- (nrow(Tabla_2014_2017)-1)-1
umbral_p <- qchisq(0.95,gl_p)

cat("Chi2 =",x2_p,"\n")
## Chi2 = 0.1414731
cat("Umbral =",umbral_p,"\n")
## Umbral = 3.841459