1. Librerías

library(dplyr)
library(knitr)
library(kableExtra)
library(DT)

2. Carga de datos

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 Acumulada. Se seleccionó como variable dependiente porque representa la cantidad total de producción obtenida por un pozo durante toda su vida útil. Es el resultado que se desea analizar y explicar a partir de otra variable.

Variable Independiente (X): Años Activo. Se eligió como variable independiente porque indica el tiempo que un pozo ha permanecido en operación. En general, mientras más años ha estado activo un pozo, mayor es la producción que puede acumular, por lo que se espera que esta variable influya directamente sobre la producción acumulada.

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

nrow(datos)
## [1] 47757
str(datos)
## 'data.frame':    47757 obs. of  2 variables:
##  $ YEARS_ACTIVE         : num  55 55 47 20 28 55 20 48 48 55 ...
##  $ CUMULATIVE_PRODUCTION: num  47225 275063 82624 7544 681006 ...
summary(datos)
##   YEARS_ACTIVE  CUMULATIVE_PRODUCTION
##  Min.   : 1.0   Min.   :     1       
##  1st Qu.:11.0   1st Qu.: 13030       
##  Median :19.0   Median : 51074       
##  Mean   :24.3   Mean   :144072       
##  3rd Qu.:36.0   3rd Qu.:172967       
##  Max.   :89.0   Max.   :985283

4. Tabla de Tripleta de Valores

4.1 Tabla de valores

nrow(datos)
## [1] 47757
condensado <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    n_valores = n(),
    valores   = paste(CUMULATIVE_PRODUCTION, collapse = ", "),
    .groups   = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

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

datatable(condensado_html, escape = FALSE, rownames = FALSE,
          colnames = c("Años Activos", "Producción Acumulada (valores)"),
          class = "stripe hover compact",
          options = list(pageLength = 10))
plot(datos$YEARS_ACTIVE, datos$CUMULATIVE_PRODUCTION,
     main = "Nube de puntos (todos los valores, sin agrupar)",
     xlab = "Años Activos (X)",
     ylab = "Producción Acumulada (Y)",
     col  = "#5b6b8c",
     pch  = 19,
     cex  = 0.6)
grid(col = "gray88")

La nube de puntos visualizada en la parte de arriba no permite conjeturar un modelo con claridad, 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 pares únicos
  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 pares (X, Y) usando la mediana

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

datatable(tabla_mediana, rownames = FALSE,
          colnames = c("Años Activos", "Mediana de Producción Acumulada"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = "CUMULATIVE_PRODUCTION", digits = 2)
write.csv(condensado[, c("YEARS_ACTIVE", "valores")], "tabla_condensada_por_anios.csv", row.names = FALSE)
write.csv(tabla_mediana, "tabla_mediana_por_anios.csv", row.names = FALSE)

5. Gráfica de Dispersión

plot(tabla_mediana$YEARS_ACTIVE, tabla_mediana$CUMULATIVE_PRODUCTION,
     main = "Producción acumulada en función de los años activos",
     xlab = "Años Activos (X)",
     ylab = "Mediana de Producción Acumulada (Y)",
     col  = "#1b2a4a",
     pch  = 19,
     cex  = 1.2)
grid(col = "gray88")

6. Conjetura

La curva crece cerca de cero y se acelera con los años, sin bajar ni aplanarse. Se conjetura un modelo potencial:

\[y = a \cdot x^{b}\]

7. Cálculo de Parámetros

\[y = a \cdot x^{b}\]

El ajuste se realizó mediante regresión no lineal utilizando los datos originales, sin aplicar transformaciones logarítmicas. Esto permitió que el modelo se adaptara mejor a los valores altos de la variable Y, logrando un ajuste más representativo del comportamiento observado en los datos.

datos_potencial <- tabla_mediana %>%
  filter(YEARS_ACTIVE > 0, CUMULATIVE_PRODUCTION > 0)

# Valores iniciales (a, b) a partir del ajuste log-log, usados solo
# como punto de partida para el algoritmo no lineal
modelo_log <- lm(log(CUMULATIVE_PRODUCTION) ~ log(YEARS_ACTIVE), data = datos_potencial)
b_inicial <- unname(coef(modelo_log)[2])
a_inicial <- unname(exp(coef(modelo_log)[1]))

modelo_nls <- nls(
  CUMULATIVE_PRODUCTION ~ a * YEARS_ACTIVE^b,
  data  = datos_potencial,
  start = list(a = a_inicial, b = b_inicial)
)
summary(modelo_nls)
## 
## Formula: CUMULATIVE_PRODUCTION ~ a * YEARS_ACTIVE^b
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a   16.058     11.723    1.37    0.174    
## b    2.321      0.169   13.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 64910 on 87 degrees of freedom
## 
## Number of iterations to convergence: 30 
## Achieved convergence tolerance: 4.216e-06
a <- unname(coef(modelo_nls)["a"])
b <- unname(coef(modelo_nls)["b"])

# R² del nls, calculado sobre la escala original (1 - SSres/SStot)
pred_nls   <- predict(modelo_nls)
ss_res     <- sum((datos_potencial$CUMULATIVE_PRODUCTION - pred_nls)^2)
ss_tot     <- sum((datos_potencial$CUMULATIVE_PRODUCTION - mean(datos_potencial$CUMULATIVE_PRODUCTION))^2)
r2_nls     <- 1 - ss_res / ss_tot

kable(data.frame(Parámetro = c("a", "b", "R² (escala original)"),
                  Valor = round(c(a, b, r2_nls), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
Parámetro Valor
a 16.0581
b 2.3210
R² (escala original) 0.8349

8. Realidad y modelo

plot(datos_potencial$YEARS_ACTIVE, datos_potencial$CUMULATIVE_PRODUCTION,
     col  = "#5b6b8c",
     pch  = 19,
     cex  = 1.2,
     main = "Modelo de producción acumulada (mediana) en función de los años activos",
     xlab = "Años Activos (X)",
     ylab = "Mediana de Producción Acumulada (Y)")
grid(col = "gray88")

curve(a * x^b, add = TRUE, col = "#1b2a4a", lwd = 2)

legend("topleft",
       legend = c("Datos reales", paste0("y = ", round(a, 2), " x^", round(b, 2))),
       col    = c("#5b6b8c", "#1b2a4a"),
       pch    = c(19, NA),
       lty    = c(NA, 1),
       lwd    = c(NA, 2),
       bty    = "n")

9. Test

El coeficiente de correlación de Pearson (r) se calculó utilizando los datos originales de años activos y producción acumulada. Por su parte, el coeficiente de determinación (R²) se obtuvo a partir del modelo de regresión no lineal (nls), ajustado con los datos originales, a partir del cual se estimaron los parámetros a y b.

r  <- cor(datos_potencial$YEARS_ACTIVE, datos_potencial$CUMULATIVE_PRODUCTION)
r2 <- r2_nls

resultado_test <- ifelse(abs(r) > 0.7, "Aceptado", "Rechazado")

tabla_test <- data.frame(
  Test      = "Correlación de Pearson (r = cor(x, y))",
  `r`       = round(r, 4),
  `abs(r)`  = round(abs(r), 4),
  `R²`      = round(r2, 4),
  Criterio  = "abs(r) > 0.7",
  Resultado = resultado_test,
  check.names = FALSE
)

kable(tabla_test, align = "lrrrrc", format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  column_spec(6, bold = TRUE,
              color = ifelse(resultado_test == "Aceptado", "#1b2a4a", "#6b6b6b"))
Test r abs(r) R² Criterio Resultado
Correlación de Pearson (r = cor(x, y)) 0.8576 0.8576 0.8349 abs(r) > 0.7 Aceptado

10. Restricciones

Dominios de las variables:

tabla_dominios <- data.frame(
  Variable = c("Y (CUMULATIVE_PRODUCTION)", "X (YEARS_ACTIVE)"),
  Dominio  = c("R+ : Y ∈ (0, +∞)", "Z+ : X ∈ 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, +∞)
X (YEARS_ACTIVE) Z+ : X ∈ Z+}

¿Existen valores de X que, utilizando la ecuación del modelo, den valores de Y fuera de su dominio?

En este caso no se presentan restricciones para el modelo, ya que el dominio de la variable independiente, Años Activos, está formado por los números enteros positivos. Al evaluar el valor mínimo del dominio (X=1), se obtiene un valor positivo de la variable dependiente, Producción Acumulada, el cual pertenece a su dominio. Además, al tratarse de una ecuación potencial con exponente positivo, la función es creciente, por lo que a medida que aumentan los años activos también aumenta la producción acumulada. En consecuencia, los valores de Y siempre serán reales positivos, cumpliendo con el dominio definido para la variable dependiente.

x_limite <- 0
cat("X =", x_limite, "\n")
## X = 0
x_min_ajuste <- min(datos_potencial$YEARS_ACTIVE)
cat("Menor valor de X usado en el ajuste:", x_min_ajuste, "\n")
## Menor valor de X usado en el ajuste: 1

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)", "X (YEARS_ACTIVE)"),
  Valido   = c("Y > 0", paste0("X ≥ ", x_min_ajuste)),
  Invalido = c("Y ≤ 0", paste0("X < ", x_min_ajuste, "  (0 ; ", x_min_ajuste, ")"))
)

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 = "5cm") %>%
  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
X (YEARS_ACTIVE) X ≥ 1 X < 1 (0 ; 1)

11. Estimación

¿Cuánta producción acumulada se espera para un pozo con 10, 30 y 60 años activos?

Se responde reemplazando cada valor de X en la ecuación del modelo \(y = a \cdot x^{b}\):

  • Para X = 10 años activos: \(y = 16.058 \cdot 10^{2.321} = 3,362.53\) de producción acumulada.

  • Para X = 30 años activos: \(y = 16.058 \cdot 30^{2.321} = 43,057.84\) de producción acumulada.

  • Para X = 60 años activos: \(y = 16.058 \cdot 60^{2.321} = 215,146.8\) de producción acumulada.

nuevos_anios <- c(5, 10, 15, 20, 25, 30)

estimaciones <- data.frame(
  YEARS_ACTIVE = nuevos_anios,
  Produccion_Estimada = a * nuevos_anios^b
)

kable(estimaciones, digits = 2, format.args = list(big.mark = ","),
      col.names = c("Años Activos", "Producción Acumulada Estimada"),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
Años Activos Producción Acumulada Estimada
5 672.95
10 3,362.53
15 8,617.27
20 16,801.54
25 28,201.66
30 43,057.84

12. Conclusión

Entre los años activos y la producción acumulada existe una relación de tipo potencial, cuya ecuación matemática es

\[y = 16.058 \cdot x^{2.321}\]

siendo y la producción acumulada y x los años activos, y donde el modelo es válido para x > 0 (años activos positivos), rango en el que fue ajustado.

Cuando los años activos son de 45, se espera una producción acumulada de 110,345.7.

La producción acumulada está influenciada en un 83.5% por los años activos, y en un 16.5% por otros factores.

El test de Pearson sobre los datos originales (r = 0.858) fue Aceptado.