Importamos el archivo “database.xls” que contiene el histórico de
incidentes en oleoductos. Utilizamos clean_names() para
asegurar que el manejo de las columnas sea eficiente.
library(readxl)
library(janitor)
library(tidyverse)
library(knitr)
library(DT)
datos <- read_excel("database.xls") %>% clean_names()
Se identifican y extraen las columnas correspondientes al volumen de pérdida no intencional y los costos totales asociados.
v_barrels <- datos$unintentional_release_barrels
v_costs <- datos$all_costs
Definimos nuestras variables de estudio: Barriles (X) como causa y Costo Total (Y) como efecto.
x_indep <- v_barrels
y_depen <- v_costs
Consolidamos los datos filtrando valores NA y ceros. Presentamos una tabla interactiva para facilitar la exploración de los datos.
tabla_pares <- data.frame(X = x_indep, Y = y_depen) %>%
filter(!is.na(X), !is.na(Y), X > 0, Y > 0)
datatable(tabla_pares,
colnames = c("Barriles (X)", "Costos USD (Y)"),
options = list(pageLength = 10))
La nube de puntos en escala logarítmica permite visualizar la relación potencial de forma lineal.
ggplot(tabla_pares, aes(x = X, y = Y)) +
geom_point(alpha = 0.4, color = "dodgerblue4") +
scale_x_log10() + scale_y_log10() +
labs(title = "Nube de Puntos (Log-Log)", x = "Log10: Barriles", y = "Log10: Costos") +
theme_minimal()
Tras observar la tendencia lineal en logaritmos, conjeturamos un modelo de la forma: \[Y = \alpha \cdot X^{\beta}\] Este modelo es ideal para representar costos que crecen de forma no lineal ante desastres ambientales. # 7. CALCULAR PARÁMETROS (LOGARITMO X Y LOGARITMO Y) Procedemos a la linealización mediante logaritmos neperianos para obtener los parámetros mediante el método de mínimos cuadrados.
log_x <- log(tabla_pares$X)
log_y <- log(tabla_pares$Y)
modelo_lineal <- lm(log_y ~ log_x)
beta <- coef(modelo_lineal)[2]
alpha <- exp(coef(modelo_lineal)[1])
Para el modelo ajustado, se han obtenido los siguientes valores fundamentales:
Superponemos la curva del modelo calculado sobre los datos reales. Se utilizan escalas logarítmicas para visualizar correctamente todos los órdenes de magnitud.
x_curva <- seq(min(tabla_pares$X), max(tabla_pares$X), length.out = 500)
df_curva <- data.frame(X = x_curva, Y = alpha * (x_curva^beta))
ggplot(tabla_pares, aes(x = X, y = Y)) +
geom_point(alpha = 0.6, color = "steelblue", size = 2) +
geom_line(data = df_curva, color = "firebrick", size = 1.2) +
scale_x_log10() + scale_y_log10() +
labs(title = "Ajuste del Modelo Potencial",
subtitle = paste("Ecuación: Y =", format(round(alpha, 2), big.mark = ","), "* X^", round(beta, 2))) +
theme_minimal()
Calculamos el coeficiente de correlación de Pearson sobre los datos transformados para determinar la fuerza de la relación lineal entre las variables logarítmicas.
r_pearson <- cor(log_x, log_y)
cat("Correlación de Pearson (r):", round(r_pearson, 4))
## Correlación de Pearson (r): 0.5111
En este paso, generamos el resumen estadístico detallado del modelo lineal. La función summary() nos permite validar la calidad del ajuste mediante el Error Estándar Residual, los p-valores (para confirmar la significancia de los coeficientes) y el Estadístico F. Además, extraemos el coeficiente de determinación \(R^2\) para cuantificar qué porcentaje de la variabilidad de los costos es explicada por el volumen de barriles.
# Generamos el resumen estadístico con las variables del paso 7
summary(modelo_lineal)
##
## Call:
## lm(formula = log_y ~ log_x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8162 -1.1939 0.0235 1.2398 6.9810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.50934 0.04398 216.23 <2e-16 ***
## log_x 0.45996 0.01478 31.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.061 on 2740 degrees of freedom
## Multiple R-squared: 0.2612, Adjusted R-squared: 0.2609
## F-statistic: 968.8 on 1 and 2740 DF, p-value: < 2.2e-16
R²: 0.2612 ( 26.12 %)
Finalmente, utilizamos el modelo para realizar una estimación práctica. ¿Cuál sería el costo esperado para un incidente de 500 barriles?
x_estimar <- 500
y_estimado <- alpha * (x_estimar^beta)
cat("Para un derrame de", x_estimar, "barriles, el costo estimado es de $",
format(round(y_estimado, 2), big.mark = ","))
## Para un derrame de 500 barriles, el costo estimado es de $ 235,116.9
Entre la variable independiente Volumen Liberado (X) y la variable dependiente Costos Totales (Y) existe una relación matemática de tipo regresión potencial, la cual indica que el impacto económico crece de forma no lineal conforme aumenta la magnitud del evento.
Esta relación se expresa mediante la fórmula del modelo: \[Y = e^{9.5093} \cdot X^{0.4600}\] (donde el exponente \(\beta = 0.4600\) determina la curvatura de la respuesta económica), sujeta a las restricciones de incluir únicamente valores de liberación y costos mayores a cero, y tras haber aplicado un Test de coeficiente de Pearson de 0.5111.
Finalmente, el modelo permite realizar una estimación técnica en la que, ante un escenario de 500 barriles liberados, el costo económico total proyectado es de $ 235,116.90.