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 razón se optó por realizar una estrategia para depurar los datos y obtener un gráfico más 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 Tripletas (X1, X2, Y) usando la mediana

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

Dominios de las variables

tabla_dominios <- data.frame(
  Variable = c("Y (CUMULATIVE_PRODUCTION)", "X1 (AVG_PRODUCTION)", "X2 (YEARS_ACTIVE)"),
  Dominio = c("R+ : Y ∈ (0, +∞)", "R+ : X1 ∈ (0, +∞)", "Z+ : X2 ∈ {0 , Z+}")
)

kable(tabla_dominios,
      col.names = c("Variable", "Dominio"),
      align = "lc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  column_spec(1, bold = TRUE, width = "6cm") %>%
  column_spec(2, width = "6cm") %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")
Variable Dominio
Y (CUMULATIVE_PRODUCTION) R+ : Y ∈ (0, +∞)
X1 (AVG_PRODUCTION) R+ : X1 ∈ (0, +∞)
X2 (YEARS_ACTIVE) Z+ : X2 ∈ {0 , Z+}

¿Existen valores de X1 y X2 que, utilizando la ecuación del modelo, nos den valores fuera del dominio de Y?

Para responder esta pregunta se especifica primero qué representa cada variable:

  • Y (CUMULATIVE_PRODUCTION): variable dependiente del modelo; representa la producción acumulada del pozo. Su dominio es R+, es decir, valores reales mayores a 0.
  • X1 (AVG_PRODUCTION): primera variable independiente; representa la producción promedio del pozo. Su dominio es R+, es decir, valores reales mayores a 0.
  • X2 (YEARS_ACTIVE): segunda variable independiente; representa los años de actividad del pozo. Su dominio es Z+, es decir, valores enteros mayores o iguales a 0.

Para encontrar los valores de X1 y X2 que producen resultados de Y fuera de su dominio, se despeja la ecuación del modelo \(y = b_0 + b_1 x_1 + b_2 x_2\) igualando Y a 0 (el límite del dominio de Y): primero se despeja X1, dejando X2 en cero, y luego se despeja X2, dejando X1 en cero.

Despeje para X1 (con X2 = 0):

\[0 = b_0 + b_1 X_1 \;\Rightarrow\; X_1 = -\frac{b_0}{b_1}\]

x1_limite <- -b0 / b1
cat("X1 =", round(x1_limite, 4), "\n")
## X1 = 3448.633

Despeje para X2 (con X1 = 0):

\[0 = b_0 + b_2 X_2 \;\Rightarrow\; X_2 = -\frac{b_0}{b_2}\]

x2_limite <- -b0 / b2
cat("X2 =", round(x2_limite, 4), "\n")
## X2 = 37.2029

Estos valores marcan el punto exacto en el que Y = 0: por debajo de ellos el modelo arroja valores de Y fuera de su dominio (Y ≤ 0). Por lo tanto, sí existen valores de X1 y X2 —los menores a los calculados arriba— que, usando la ecuación del modelo, dan como resultado un Y fuera de su dominio. A partir de estos límites se construyen las restricciones prácticas del modelo.

Con base en el análisis realizado, se determinaron los rangos en los que el modelo presenta un desempeño adecuado, así como los intervalos en los que deja de ser efectivo.

Restricciones del modelo

tabla_restricciones <- data.frame(
  Variable = c("Y (CUMULATIVE_PRODUCTION)", "X1 (AVG_PRODUCTION)", "X2 (YEARS_ACTIVE)"),
  Valido = c("Y > 0", "X1 ≥ 3499", "X2 ≥ 38"),
  Invalido = c("Y ≤ 0", "X1 < 3499 (0 ; 3499)", "X2 < 38 (0 ; 38)")
)

kable(tabla_restricciones,
      col.names = c("Variable", "Rango válido", "Rango inválido"),
      align = "lcc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  column_spec(1, bold = TRUE, width = "6cm") %>%
  column_spec(2, width = "4cm") %>%
  column_spec(3, width = "5cm") %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")
Variable Rango válido Rango inválido
Y (CUMULATIVE_PRODUCTION) Y > 0 Y ≤ 0
X1 (AVG_PRODUCTION) X1 ≥ 3499 X1 < 3499 (0 ; 3499)
X2 (YEARS_ACTIVE) X2 ≥ 38 X2 < 38 (0 ; 38)

11. Estimación

Calculadora interactiva del modelo

Ingresa un valor de X1 (Producción Promedio) y X2 (Años Activos) para estimar Y (Producción Acumulada) usando la ecuación del modelo. El resultado se marca como inválido si los valores no respetan el dominio de las variables, si están fuera del rango en el que el modelo es confiable, o si el resultado de Y es negativo.

Calculadora del modelo

Y = -145857.3797 + 42.2943 · X1 + 3920.5903 · X2

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 (R²); el resto (7.9%) es por otros factores, evidenciando una relación lineal positiva y de alta intensidad entre las variables.