Regresion Lineal Multiple - Oil & Gas Leases

1. Librerías

suppressMessages(suppressWarnings({
  library(dplyr)
  library(knitr)
  library(kableExtra)
  library(DT)
  library(plotly)
}))

2. Carga de datos

# usar Run All, no Knit
ruta <- file.choose()
datos_completos <- read.csv(ruta, stringsAsFactors = FALSE)

cat("Filas:", nrow(datos_completos), "| Columnas:", ncol(datos_completos), "\n")

## Filas: 47757 | Columnas: 24

str(datos_completos)

## 'data.frame':    47757 obs. of  24 variables:
##  $ KID                      : int  1001106903 1001106572 1001106590 1001107343 1001108234 1001106684 1001107377 1001107386 1001107740 1001106710 ...
##  $ DEPTH_OF_WELL            : num  700 800 1400 1125 2940 ...
##  $ CUMULATIVE_PRODUCTION    : num  47225 275063 82624 7544 681006 ...
##  $ AVG_PRODUCTION           : num  859 5001 1758 377 24322 ...
##  $ LATITUDE                 : num  37.1 38.8 37.5 37.8 37.1 ...
##  $ LONGITUDE                : num  -95.9 -95.2 -96.3 -95.7 -101.3 ...
##  $ YEARS_ACTIVE             : num  55 55 47 20 28 55 20 48 48 55 ...
##  $ SECTION                  : num  33 11 34 8 30 4 26 28 11 17 ...
##  $ COUNTY_CODE              : num  125 45 49 207 189 121 49 1 31 121 ...
##  $ STATE_CODE               : int  15 15 15 15 15 15 15 15 15 15 ...
##  $ TOWNSHIP                 : num  33 15 29 26 33 17 30 26 23 16 ...
##  $ RANGE                    : num  14 20 10 16 36 25 12 21 16 24 ...
##  $ PRODUCES_OIL             : num  1 1 1 1 0 1 1 1 1 1 ...
##  $ PRODUCES_GAS             : num  0 0 0 0 1 0 0 0 0 0 ...
##  $ OPERATOR_NAME            : chr  "Horton, John" "Whitlow Energy, Inc." "Suerte Oil Company" "Patterson-Blackford" ...
##  $ FIELD_NAME               : chr  "WAYSIDE-HAVANA" "BALDWIN" "DUNKLEBERGER" "ROSE EAST" ...
##  $ PRODUCING_FORMATION      : chr  "UNKNOWN" "UNKNOWN" "UNKNOWN" "UNKNOWN" ...
##  $ LONGITUDE_LATITUDE_SOURCE: chr  "CENTER_OF_SECTION" "CENTER_OF_SECTION" "CENTER_OF_SECTION" "CENTER_OF_SECTION" ...
##  $ PROD_LEVEL               : chr  "MEDIUM" "HIGH" "MEDIUM" "LOW" ...
##  $ DEPTH_LEVEL              : chr  "SHALLOW" "SHALLOW" "SHALLOW" "SHALLOW" ...
##  $ LIFE_STAGE               : chr  "OLD" "OLD" "OLD" "MATURE" ...
##  $ AVG_PROD_LEVEL           : chr  "LOW" "MEDIUM" "MEDIUM" "LOW" ...
##  $ TOWNSHIP_DIRECTION       : chr  "S" "S" "S" "S" ...
##  $ RANGE_DIRECTION          : chr  "E" "E" "E" "E" ...

3. Selección de variables

Variable Dependiente (Y): Producción Promedio. Se seleccionó como variable de respuesta porque representa la producción promedio estimada de un pozo, es decir, su desempeño operativo actual.

Variable Independiente 1 (X₁): Años Activo. Indica el tiempo que el pozo ha estado en operación. Esta variable refleja la madurez del pozo y su experiencia operativa. Es la variable dominante, ya que determina de manera significativa la capacidad de producción promedio actual.

Variable Independiente 2 (X₂): Producción Acumulada. Representa la producción acumulada del pozo desde su inicio. Esta variable cuantifica el efecto histórico del rendimiento del pozo sobre su producción promedio actual, aportando información complementaria sobre su desempeño operativo.

La Producción Promedio (Y) es, esencialmente, el resultado de la interacción entre la antigüedad del pozo (X₁) y su producción acumulada (X₂). Ambos factores interactúan para determinar la producción promedio estimada, proporcionando un panorama integral del desempeño del pozo.

datos <- datos_completos %>%
  mutate(
    CUMULATIVE_PRODUCTION = as.numeric(CUMULATIVE_PRODUCTION),
    AVG_PRODUCTION        = as.numeric(AVG_PRODUCTION),
    YEARS_ACTIVE          = as.numeric(YEARS_ACTIVE)
  ) %>%
  select(CUMULATIVE_PRODUCTION, AVG_PRODUCTION, YEARS_ACTIVE) %>%
  na.omit() %>%
  filter(CUMULATIVE_PRODUCTION > 0, AVG_PRODUCTION > 0, YEARS_ACTIVE > 0)

nrow(datos)

## [1] 47757

4. Tabla de Tripleta de Valores

4.1 Tripleta de valores

condensado <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    n_valores = n(),
    valores_avg = paste(AVG_PRODUCTION, collapse = ", "),
    valores_cum = paste(CUMULATIVE_PRODUCTION, collapse = ", "),
    .groups = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

condensado_html <- condensado %>%
  mutate(
    AVG_PRODUCTION = paste0("<details><summary>", n_valores, " valor(es)</summary>", valores_avg, "</details>"),
    CUMULATIVE_PRODUCTION = paste0("<details><summary>", n_valores, " valor(es)</summary>", valores_cum, "</details>")
  ) %>%
  select(YEARS_ACTIVE, AVG_PRODUCTION, CUMULATIVE_PRODUCTION)

nrow(datos)

## [1] 47757

datatable(condensado_html, escape = FALSE, rownames = FALSE,
          colnames = c("Años Activos (X2)", "Producción Promedio (valores, X1)", "Producción Acumulada (valores, Y)"),
          class = "stripe hover compact",
          options = list(pageLength = 10))

plot_ly(
  data = datos,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers",
  marker = list(size = 2, opacity = 0.4)
) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

La nube de puntos visualizada en la parte de arriba no permite conjeturar un modelo, por esta razon se opto por realizar una estrategia para depurar los datos y obtener un grafico mas limpio.

Estrategia

1. Generación de tripletas únicas
2. Llenar espacios donde no haya valores
3. Donde se encuentren valores diversos usar indicadores de posición (mediana)
4. Presentar nueva tabla con la mediana de los valores

4.2 Tabla de tripletas (X1, X2, Y) usando la mediana por YEARS_ACTIVE

tabla_mediana <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    AVG_PRODUCTION        = median(AVG_PRODUCTION),
    CUMULATIVE_PRODUCTION = median(CUMULATIVE_PRODUCTION),
    .groups = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

nrow(tabla_mediana)

## [1] 89

datatable(tabla_mediana, rownames = FALSE,
          colnames = c("Años Activos (X2)", "Mediana Producción Promedio (X1)", "Mediana Producción Acumulada (Y)"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = c("AVG_PRODUCTION", "CUMULATIVE_PRODUCTION"), digits = 2)

write.csv(condensado[, c("YEARS_ACTIVE", "valores_avg", "valores_cum")], "tabla_condensada_por_anios_oilgas.csv", row.names = FALSE)
write.csv(tabla_mediana, "tabla_mediana_por_anios_oilgas.csv", row.names = FALSE)

5. Gráfica de Dispersión

plot_ly(
  data = tabla_mediana,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers"
) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

6. Conjetura

Se conjetura que la producción acumulada (Y) mantiene una relación lineal positiva con la producción promedio (X1) y con los años de actividad (X2). En el gráfico se observa que, a medida que aumenta la producción promedio y el tiempo de operación, la producción acumulada tiende a incrementarse. Aunque existen algunos valores atípicos, se decidió usar un modelo de regresión lineal múltiple.

\[y = b_0 + b_1 x_1 + b_2 x_2\]

7. Cálculo de Parámetros

modelo_lineal <- lm(CUMULATIVE_PRODUCTION ~ AVG_PRODUCTION + YEARS_ACTIVE, data = tabla_mediana)
summary(modelo_lineal)

## 
## Call:
## lm(formula = CUMULATIVE_PRODUCTION ~ AVG_PRODUCTION + YEARS_ACTIVE, 
##     data = tabla_mediana)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -111120  -28177  -13323   29092  169373 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.459e+05  1.145e+04  -12.73   <2e-16 ***
## AVG_PRODUCTION  4.229e+01  2.978e+00   14.20   <2e-16 ***
## YEARS_ACTIVE    3.921e+03  2.093e+02   18.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45170 on 86 degrees of freedom
## Multiple R-squared:  0.9209, Adjusted R-squared:  0.9191 
## F-statistic: 500.9 on 2 and 86 DF,  p-value: < 2.2e-16

b0 <- unname(coef(modelo_lineal)[1])
b1 <- unname(coef(modelo_lineal)[2])
b2 <- unname(coef(modelo_lineal)[3])

c(b0 = b0, b1 = b1, b2 = b2)

##            b0            b1            b2 
## -145857.37971      42.29425    3920.59028

kable(data.frame(Parámetro = c("b0 (intercepto)", "b1 (pendiente X1)", "b2 (pendiente X2)"),
                  Valor = round(c(b0, b1, b2), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Parámetro	Valor
b0 (intercepto)	-145857.3797
b1 (pendiente X1)	42.2943
b2 (pendiente X2)	3920.5903

8. Realidad y modelo

x1_seq <- seq(min(tabla_mediana$AVG_PRODUCTION), max(tabla_mediana$AVG_PRODUCTION), length.out = 30)
x2_seq <- seq(min(tabla_mediana$YEARS_ACTIVE), max(tabla_mediana$YEARS_ACTIVE), length.out = 30)
z_plano <- outer(x2_seq, x1_seq, function(x2, x1) b0 + b1 * x1 + b2 * x2)

plot_ly() %>%
  add_markers(data = tabla_mediana, x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION) %>%
  add_surface(x = x1_seq, y = x2_seq, z = z_plano, opacity = 0.6, showscale = FALSE) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

9. Test

r_x1 <- cor(tabla_mediana$AVG_PRODUCTION, tabla_mediana$CUMULATIVE_PRODUCTION)
r_x2 <- cor(tabla_mediana$YEARS_ACTIVE, tabla_mediana$CUMULATIVE_PRODUCTION)
r2 <- summary(modelo_lineal)$r.squared
p_modelo <- pf(summary(modelo_lineal)$fstatistic[1],
               summary(modelo_lineal)$fstatistic[2],
               summary(modelo_lineal)$fstatistic[3], lower.tail = FALSE)

resultado_x1 <- ifelse(abs(r_x1) > 0.7, "Aceptado", "Rechazado")
resultado_x2 <- ifelse(abs(r_x2) > 0.7, "Aceptado", "Rechazado")
resultado_modelo <- ifelse(p_modelo < 0.05, "Significativo", "No significativo")

tabla_test <- data.frame(
  Test = c("Pearson X1", "Pearson X2", "F-test modelo"),
  Valor = c(round(r_x1, 4), round(r_x2, 4), format(p_modelo, scientific = TRUE, digits = 4)),
  Resultado = c(resultado_x1, resultado_x2, resultado_modelo)
)

kable(tabla_test, align = "lrc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Test	Valor	Resultado
Pearson X1	0.7735	Aceptado
Pearson X2	0.8576	Aceptado
F-test modelo	4.105e-48	Significativo

cat("R² del modelo:", round(r2, 4), "\n")

## R² del modelo: 0.9209

10. Restricciones

Dominio de X1 (AVG_PRODUCTION): presenta un dominio de los R+ X1∈(0,+∞)

Dominio de X2 (YEARS_ACTIVE): presenta un dominio de los R+ X2∈(0,+∞)

Dominio de Y (CUMULATIVE_PRODUCTION): presenta un dominio de los R+ Y∈(0,+∞)

Como todos los dominios fluctúan entre R+, los resultados del modelo no pueden ser negativos.

El modelo es inválido para valores de:

X1 (AVG_PRODUCTION) desde -∞ hasta 3499

X2 (YEARS_ACTIVE) desde -∞ hasta 37.30

11. Estimación

Unidades: [COMPLETAR: indicar la unidad de CUMULATIVE_PRODUCTION y AVG_PRODUCTION según el diccionario de datos del dataset, por ejemplo barriles (bbl)]

nuevos_x1 <- c(3499, 5000, 10000, 25000, 50000, 100000)
nuevos_x2 <- c(38, 40, 45, 50, 60, 70)

estimaciones <- expand.grid(AVG_PRODUCTION = nuevos_x1, YEARS_ACTIVE = nuevos_x2)
estimaciones$Produccion_Acumulada_Estimada <- predict(modelo_lineal, newdata = estimaciones)

# se descartan estimaciones negativas
estimaciones <- estimaciones %>% filter(Produccion_Acumulada_Estimada > 0)
estimaciones <- estimaciones %>% arrange(AVG_PRODUCTION, YEARS_ACTIVE)

datatable(estimaciones, rownames = FALSE,
          colnames = c("Producción Promedio (X1) [unidad]", "Años Activos (X2)",
                       "Producción Acumulada Estimada (Y) [unidad]"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = "Produccion_Acumulada_Estimada", digits = 2)

12. Conclusión

Entre la producción promedio (x1), los años de actividad (X2) y la producción acumulada (Y) existe una relación lineal múltiple cuya ecuación matemática es:

\[Y = -145857.4 + 42.29425\,X_1 + 3920.59\,X_2\]

Siendo X₁ la producción promedio, medida en las unidades correspondientes de producción; X₂ los años de actividad del pozo; y Y la producción acumulada. Con un coeficiente de correlación múltiple, el modelo explica 92.1% de la variabilidad de Y (R2); el resto (7.9%) es por otros factores, evidenciando una relación lineal positiva y de alta intensidad entre las variables.