Relación entre Años Activos y Producción Acumulada (Análisis por Partes)

Usar Run All (Code → Run Region → Run All), no Knit, porque el chunk de carga usa file.choose().

1. Librerías

suppressMessages(suppressWarnings({
  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. Tablas

4.1 Tabla de valores

nrow(datos)

## [1] 47757

condensado <- datos %>%
  mutate(Año = floor(YEARS_ACTIVE)) %>%
  group_by(Año) %>%
  summarise(
    n_valores = n(),
    valores   = paste(round(CUMULATIVE_PRODUCTION, 2), collapse = ", "),
    .groups   = "drop"
  ) %>%
  arrange(Año)

condensado_html <- condensado %>%
  mutate(
    CUMULATIVE_PRODUCTION = paste0(
      "<details><summary>", n_valores, " valor(es)</summary>", valores, "</details>"
    )
  ) %>%
  select(Año, 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.4)
grid(col = "gray88")

La nube de puntos visualizada en la parte de arriba no permite conjeturar un modelo con claridad: al haber miles de pozos que comparten el mismo año de actividad, los puntos se amontonan y no dejan ver ninguna tendencia. Por esta razón se optó por realizar una estrategia para depurar los datos y obtener un gráfico más limpio.

Estrategia

1. Agrupar los pozos por cada año exacto de actividad (floor(YEARS_ACTIVE))
2. Calcular la media de X y de Y dentro de cada año, obteniendo un único par (x̄, ȳ) por año
3. Dividir el eje de X en tres partes iguales (según su valor máximo), ya que la nube completa no sigue un único patrón en todo su rango
4. Conjeturar y ajustar un modelo por separado en cada parte

4.2 Tabla de pares (X, Y) usando la media

tabla_media <- datos %>%
  mutate(Año = floor(YEARS_ACTIVE)) %>%
  group_by(Año) %>%
  summarise(
    x_medio = mean(YEARS_ACTIVE),
    y_medio = mean(CUMULATIVE_PRODUCTION),
    .groups = "drop"
  ) %>%
  arrange(Año)

datatable(tabla_media, rownames = FALSE,
          colnames = c("Año", "X̄ (Años Activos)", "Ȳ (Producción Acumulada)"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = c("x_medio", "y_medio"), digits = 2)

4.3 División en tres partes según el valor máximo de X

max_x <- max(tabla_media$x_medio)
min_x <- min(tabla_media$x_medio)

limite1 <- max_x / 3
limite2 <- 2 * max_x / 3

tabla_media <- tabla_media %>%
  mutate(parte = case_when(
    x_medio <= limite1 ~ "Parte 1",
    x_medio <= limite2 ~ "Parte 2",
    TRUE ~ "Parte 3"
  ))

tabla1 <- tabla_media %>% filter(parte == "Parte 1") %>% select(-parte)
tabla2 <- tabla_media %>% filter(parte == "Parte 2") %>% select(-parte)
tabla3 <- tabla_media %>% filter(parte == "Parte 3") %>% select(-parte)

data.frame(
  Parte   = c("Parte 1", "Parte 2", "Parte 3"),
  Rango_X = c(paste0("[", round(min_x,0), " - ", round(limite1,0), "]"),
              paste0("(", round(limite1,0), " - ", round(limite2,0), "]"),
              paste0("(", round(limite2,0), " - ", round(max_x,0), "]")),
  n_años  = c(nrow(tabla1), nrow(tabla2), nrow(tabla3))
) %>%
  kable(col.names = c("Parte", "Rango de X (años)", "Años promediados (n)"),
        align = "lcc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")

Parte	Rango de X (años)	Años promediados (n)
Parte 1	[1 - 30]	29
Parte 2	(30 - 59]	30
Parte 3	(59 - 89]	30

write.csv(condensado[, c("Año", "valores")], "tabla_condensada_por_anios.csv", row.names = FALSE)
write.csv(tabla_media, "tabla_media_por_anios.csv", row.names = FALSE)

5. Gráfica de Dispersión

colores <- c("Parte 1" = "#5b6b8c", "Parte 2" = "#8a7a5b", "Parte 3" = "#1b2a4a")

par(mar = c(5, 5, 4, 9))
plot(tabla_media$x_medio, tabla_media$y_medio,
     col = colores[tabla_media$parte], pch = 19, cex = 1.1,
     main = "Producción acumulada en función de los años activos, dividida en 3 partes",
     xlab = "X̄ (Años Activos)",
     ylab = "Ȳ (Producción Acumulada)",
     cex.main = 0.9, frame.plot = FALSE)
grid(col = "gray88")
abline(v = c(limite1, limite2), lty = 2, col = "gray50")
legend("topright", inset = c(-0.30, 0), legend = names(colores), col = colores,
       pch = 19, bty = "n", xpd = TRUE)
box()

6. Conjetura

Al observar la nube de puntos completa, esta resulta compleja: no sigue un único patrón a lo largo de todo el rango de X. Por esta razón se conjetura un modelo distinto para cada una de las tres partes:

• Parte 1 (pozos jóvenes, 1 - 30 años): la nube muestra un crecimiento que se va acelerando a medida que aumentan los años activos, sin señales de aplanamiento — un patrón de crecimiento compuesto. Se conjetura un modelo exponencial: \[y = a \cdot e^{bx}\]

• Parte 2 (pozos maduros, 30 - 59 años): la nube no es monótona: sube, se estabiliza y vuelve a subir (posible efecto de reacondicionamientos a mitad de vida del pozo). Ningún modelo monótono puede seguir esa curvatura, así que se conjetura un modelo polinómico de grado 2: \[y = a + bx + cx^{2}\]

• Parte 3 (pozos muy antiguos, 59 - 89 años): el crecimiento vuelve a acelerarse de forma marcada al final. Se conjetura un modelo exponencial: \[y = a \cdot e^{bx}\]

7. Cálculo de Parámetros

Parte 1 — Modelo exponencial

modelo1 <- lm(log(y_medio) ~ x_medio, data = tabla1)
a1 <- unname(exp(coef(modelo1)[1]))
b1 <- unname(coef(modelo1)[2])
cat("y =", round(a1, 3), "* e^(", round(b1, 4), "* x )\n")

## y = 19947.68 * e^( 0.0939 * x )

kable(data.frame(Parámetro = c("a", "b"), Valor = round(c(a1, b1), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Parámetro	Valor
a	19947.6829
b	0.0939

Parte 2 — Modelo polinómico (grado 2)

modelo2 <- lm(y_medio ~ x_medio + I(x_medio^2), data = tabla2)
a2 <- unname(coef(modelo2)[1])
b2 <- unname(coef(modelo2)[2])
c2 <- unname(coef(modelo2)[3])

signo_b2 <- ifelse(b2 >= 0, "+", "-")
signo_c2 <- ifelse(c2 >= 0, "+", "-")
b2_txt <- paste(signo_b2, round(abs(b2), 2))
c2_txt <- paste(signo_c2, round(abs(c2), 4))

cat("y =", round(a2,2), b2_txt, "*x", c2_txt, "*x^2\n")

## y = 1557110 - 60147.85 *x + 654.2263 *x^2

kable(data.frame(Parámetro = c("a", "b", "c"), Valor = round(c(a2, b2, c2), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Parámetro	Valor
a	1557110.0645
b	-60147.8532
c	654.2263

Parte 3 — Modelo exponencial

modelo3 <- lm(log(y_medio) ~ x_medio, data = tabla3)
a3 <- unname(exp(coef(modelo3)[1]))
b3 <- unname(coef(modelo3)[2])
cat("y =", round(a3, 3), "* e^(", round(b3, 4), "* x )\n")

## y = 56824.06 * e^( 0.0259 * x )

kable(data.frame(Parámetro = c("a", "b"), Valor = round(c(a3, b3), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Parámetro	Valor
a	56824.0626
b	0.0259

8. Realidad y modelo

plot(tabla1$x_medio, tabla1$y_medio, col = "#5b6b8c", pch = 19, cex = 1.2,
     main = "Parte 1 — Modelo exponencial vs. realidad",
     xlab = "X̄ (Años Activos)", ylab = "Ȳ (Producción Acumulada)",
     cex.main = 0.9, frame.plot = FALSE)
grid(col = "gray88")
curve(a1 * exp(b1 * x), add = TRUE, col = "#1b2a4a", lwd = 2)
legend("topleft", legend = c("Datos (x̄, ȳ)", paste0("y=", round(a1,2), "e^(", round(b1,4), "x)")),
       col = c("#5b6b8c", "#1b2a4a"), pch = c(19, NA), lty = c(NA, 1), lwd = c(NA, 2), bty = "n")
box()

plot(tabla2$x_medio, tabla2$y_medio, col = "#8a7a5b", pch = 19, cex = 1.2,
     main = "Parte 2 — Modelo polinómico vs. realidad",
     xlab = "X̄ (Años Activos)", ylab = "Ȳ (Producción Acumulada)",
     cex.main = 0.9, frame.plot = FALSE)
grid(col = "gray88")
curve(a2 + b2*x + c2*x^2, add = TRUE, col = "#1b2a4a", lwd = 2)
legend("topleft", legend = c("Datos (x̄, ȳ)", "Modelo Polinómico"),
       col = c("#8a7a5b", "#1b2a4a"), pch = c(19, NA), lty = c(NA, 1), lwd = c(NA, 2), bty = "n")
box()

plot(tabla3$x_medio, tabla3$y_medio, col = "#1b2a4a", pch = 19, cex = 1.2,
     main = "Parte 3 — Modelo exponencial vs. realidad",
     xlab = "X̄ (Años Activos)", ylab = "Ȳ (Producción Acumulada)",
     cex.main = 0.9, frame.plot = FALSE)
grid(col = "gray88")
curve(a3 * exp(b3 * x), add = TRUE, col = "#5b6b8c", lwd = 2)
legend("topleft", legend = c("Datos (x̄, ȳ)", paste0("y=", round(a3,2), "e^(", round(b3,4), "x)")),
       col = c("#1b2a4a", "#5b6b8c"), pch = c(19, NA), lty = c(NA, 1), lwd = c(NA, 2), bty = "n")
box()

9. Test

El coeficiente de correlación de Pearson (r) se calcula sobre la escala en que cada modelo es lineal (X vs. log(Y) para los exponenciales). Para el modelo polinómico, al no ser lineal en una sola escala, se usa el coeficiente de correlación múltiple R = √(R²), que cumple el mismo rol de bondad de ajuste entre 0 y 1.

r1 <- cor(tabla1$x_medio, log(tabla1$y_medio))
R2_2 <- summary(modelo2)$r.squared
r2 <- sqrt(R2_2)
r3 <- cor(tabla3$x_medio, log(tabla3$y_medio))

resultado1 <- ifelse(abs(r1) > 0.7, "Aceptado", "Rechazado")
resultado2 <- ifelse(r2 > 0.7, "Aceptado", "Rechazado")
resultado3 <- ifelse(abs(r3) > 0.7, "Aceptado", "Rechazado")

tabla_test <- data.frame(
  Parte      = c("Parte 1 (Exponencial)", "Parte 2 (Polinómico)", "Parte 3 (Exponencial)"),
  r_o_R      = round(c(r1, r2, r3), 4),
  r_absoluto = round(c(abs(r1), r2, abs(r3)), 4),
  Criterio   = "r absoluto > 0.7",
  Resultado  = c(resultado1, resultado2, resultado3),
  check.names = FALSE
)

kable(tabla_test,
      col.names = c("Parte", "r / R", "r absoluto", "Criterio", "Resultado"),
      align = "lcccc",
      format = "html", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  column_spec(5, bold = TRUE,
              color = ifelse(tabla_test$Resultado == "Aceptado", "#1b2a4a", "#6b6b6b"))

Parte	r / R	r absoluto	Criterio	Resultado
Parte 1 (Exponencial)	0.9506	0.9506	r absoluto > 0.7	Aceptado
Parte 2 (Polinómico)	0.8434	0.8434	r absoluto > 0.7	Aceptado
Parte 3 (Exponencial)	0.9248	0.9248	r absoluto > 0.7	Aceptado

10. Restricciones

Dominios de las variables:

tabla_dominios <- data.frame(
  Parte   = c("Parte 1", "Parte 2", "Parte 3"),
  Modelo  = c("Exponencial", "Polinómico (grado 2)", "Exponencial"),
  Dominio = c("R+ : X,Y ∈ (0, +∞)", "R+ : X,Y ∈ (0, +∞)", "R+ : X,Y ∈ (0, +∞)")
)

kable(tabla_dominios,
      col.names = c("Parte", "Modelo", "Dominio"),
      align = "llc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")

Parte	Modelo	Dominio
Parte 1	Exponencial	R+ : X,Y ∈ (0, +∞)
Parte 2	Polinómico (grado 2)	R+ : X,Y ∈ (0, +∞)
Parte 3	Exponencial	R+ : X,Y ∈ (0, +∞)

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

En los tres modelos no se identificaron restricciones dentro del rango de X observado (1 a 89 años), ya que la variable independiente (Años Activos) solo toma valores positivos y, para ese rango, las tres ecuaciones producen siempre valores de Y dentro de su dominio (Y > 0).

11. Estimación

Ingresa un valor de Años Activos (X) y elige a qué parte pertenece (según el rango de X); la calculadora cambia automáticamente de ecuación según la parte elegida. No se permiten valores negativos de X.

Calculadora Interactiva de Estimación

Selecciona la parte / modelo:
Parte 1 — Exponencial (jóvenes, 1 - 29 años) Parte 2 — Polinómico grado 2 (maduros, 30 - 59 años) Parte 3 — Exponencial (muy antiguos, 60 - 89 años)

X (Años Activos):
Calcular Y

12. Conclusión

Entre los años activos de un pozo (X) y su producción acumulada (Y) existe relación, pero no es la misma en todo el rango — por eso la nube se dividió en tres partes.

Parte 1 (pozos jóvenes, 1 - 30 años): relación exponencial

\[y = 19947.68 \cdot e^{0.0939x}\]

con r = 0.951 (Aceptado).

Parte 2 (pozos maduros, 30 - 59 años): relación polinómica de grado 2

\[y = 1557110 - 60147.85x + 654.2263x^2\]

con R = 0.843 (Aceptado). La curvatura refleja posibles reactivaciones/reacondicionamientos a mitad de vida del pozo.

Parte 3 (pozos muy antiguos, 59 - 89 años): relación exponencial

\[y = 56824.06 \cdot e^{0.0259x}\]

con r = 0.925 (Aceptado).

En conjunto, esto confirma que el tiempo de actividad del pozo es un buen predictor de su producción acumulada, aunque la forma de esa relación cambia con la madurez del pozo, lo que justifica el análisis por partes en lugar de un único modelo global.