Regresión: Años activos vs. Producción acumulada de pozos petroleros
Análisis de datos - Oil and Gas Leases
2026-06-30
1 Librerías
2 Datos
datos <- read_csv("oil_and_gas_leases_data___2_.csv", show_col_types = FALSE)
# Rellenamiento de celdas faltantes (NA) con la media de cada variable numérica de interés
if (any(is.na(datos$YEARS_ACTIVE))) {
datos$YEARS_ACTIVE[is.na(datos$YEARS_ACTIVE)] <- mean(datos$YEARS_ACTIVE, na.rm = TRUE)
}
if (any(is.na(datos$CUMULATIVE_PRODUCTION))) {
datos$CUMULATIVE_PRODUCTION[is.na(datos$CUMULATIVE_PRODUCTION)] <- mean(datos$CUMULATIVE_PRODUCTION, na.rm = TRUE)
}
# Se descartan pozos con valores no positivos (no tienen sentido físico y
# además impiden aplicar logaritmos en los modelos exponencial/potencial)
datos <- datos %>%
filter(YEARS_ACTIVE > 0, CUMULATIVE_PRODUCTION > 0)
str(datos[, c("YEARS_ACTIVE", "CUMULATIVE_PRODUCTION")])## tibble [47,757 × 2] (S3: tbl_df/tbl/data.frame)
## $ YEARS_ACTIVE : num [1:47757] 55 55 47 20 28 55 20 48 48 55 ...
## $ CUMULATIVE_PRODUCTION: num [1:47757] 47225 275063 82624 7544 681006 ...3 Selección de las variables (causa y efecto)
Variable causa (X): YEARS_ACTIVE — años
que el pozo ha permanecido activo.
Variable efecto (Y):
CUMULATIVE_PRODUCTION — producción acumulada del
pozo (barriles/pies³ equivalentes).
Justificación: la producción acumulada es, por
definición, la suma de la producción a lo largo del tiempo de vida del
pozo. Por lo tanto el tiempo de actividad (YEARS_ACTIVE) es
la variable que antecede y explica la magnitud de lo acumulado
(efecto), y no al revés: el pozo no produce más años porque acumuló más,
sino que acumula más precisamente porque lleva más años en producción.
Esto satisface la relación causa → efecto exigida.
4 Tabla de pares de variables
4.1 Máximos, mínimos y separación de la variable
max_x <- max(datos$YEARS_ACTIVE)
min_x <- min(datos$YEARS_ACTIVE)
max_y <- max(datos$CUMULATIVE_PRODUCTION)
min_y <- min(datos$CUMULATIVE_PRODUCTION)
resumen_general <- data.frame(
Variable = c("YEARS_ACTIVE (X)", "CUMULATIVE_PRODUCTION (Y)"),
Minimo = c(min_x, min_y),
Maximo = c(max_x, max_y),
Rango = c(max_x - min_x, max_y - min_y),
Media = c(mean(datos$YEARS_ACTIVE), mean(datos$CUMULATIVE_PRODUCTION))
)
kable(resumen_general, digits = 2, caption = "Resumen general de las variables")| Variable | Minimo | Maximo | Rango | Media |
|---|---|---|---|---|
| YEARS_ACTIVE (X) | 1 | 89 | 88 | 24.3 |
| CUMULATIVE_PRODUCTION (Y) | 1 | 985283 | 985282 | 144072.2 |
4.2 División en tres partes según el valor máximo de X
Se toma el valor máximo de YEARS_ACTIVE (89 años) y se
divide el eje en tres tercios iguales, de modo que cada
parte agrupe pozos “jóvenes”, “maduros” y “muy antiguos”. Esto también
permite visualizar mejor la nube de puntos, que al graficarse completa
se satura por la gran cantidad de observaciones (47757 pozos).
limite1 <- max_x / 3
limite2 <- 2 * max_x / 3
datos <- datos %>%
mutate(parte = case_when(
YEARS_ACTIVE <= limite1 ~ "Parte 1",
YEARS_ACTIVE <= limite2 ~ "Parte 2",
TRUE ~ "Parte 3"
))
parte1 <- datos %>% filter(parte == "Parte 1")
parte2 <- datos %>% filter(parte == "Parte 2")
parte3 <- datos %>% filter(parte == "Parte 3")
kable(data.frame(
Parte = c("Parte 1", "Parte 2", "Parte 3"),
Rango_X = c(paste0("[", round(min_x,1), " - ", round(limite1,1), "]"),
paste0("(", round(limite1,1), " - ", round(limite2,1), "]"),
paste0("(", round(limite2,1), " - ", round(max_x,1), "]")),
n = c(nrow(parte1), nrow(parte2), nrow(parte3))
), caption = "Pozos por parte")| Parte | Rango_X | n |
|---|---|---|
| Parte 1 | [1 - 29.7] | 31814 |
| Parte 2 | (29.7 - 59.3] | 13200 |
| Parte 3 | (59.3 - 89] | 2743 |
4.3 Construcción de la tabla de pares (clases y medias)
Con miles de pozos por parte, la nube de puntos es puro ruido. Por
eso se agrupa YEARS_ACTIVE en clases (regla de Sturges:
k = 1 + 3.322*log10(n)), se calcula la
amplitud de cada clase (rango/k) y se
obtiene el par (x̄, ȳ) -promedio de cada variable- por
clase. Estas parejas (x̄, ȳ), todas del mismo tamaño de muestra por
clase, son las que se usan para ajustar los modelos.
construir_tabla <- function(sub_datos) {
n <- nrow(sub_datos)
k <- ceiling(1 + 3.322 * log10(n)) # número de clases (Sturges)
minx <- min(sub_datos$YEARS_ACTIVE)
maxx <- max(sub_datos$YEARS_ACTIVE)
amplitud <- (maxx - minx) / k
cortes <- seq(minx, maxx, length.out = k + 1)
sub_datos$clase <- cut(sub_datos$YEARS_ACTIVE, breaks = cortes, include.lowest = TRUE)
tabla <- sub_datos %>%
group_by(clase) %>%
summarise(
n_pozos = n(),
x_medio = mean(YEARS_ACTIVE),
y_medio = mean(CUMULATIVE_PRODUCTION),
.groups = "drop"
) %>%
filter(!is.na(x_medio)) %>%
arrange(x_medio)
attr(tabla, "k") <- k
attr(tabla, "amplitud") <- amplitud
tabla
}
tabla1 <- construir_tabla(parte1)
tabla2 <- construir_tabla(parte2)
tabla3 <- construir_tabla(parte3)kable(tabla1, digits = 2, caption = paste0("Parte 1 · k = ", attr(tabla1,"k"),
" clases · amplitud = ", round(attr(tabla1,"amplitud"), 2)))| clase | n_pozos | x_medio | y_medio |
|---|---|---|---|
| [1,2.75] | 3269 | 1.33 | 10987.78 |
| (2.75,4.5] | 2149 | 3.47 | 29459.15 |
| (4.5,6.25] | 1939 | 5.52 | 39367.34 |
| (6.25,8] | 2023 | 7.50 | 46118.11 |
| (8,9.75] | 960 | 9.00 | 50316.09 |
| (9.75,11.5] | 2999 | 10.57 | 74746.10 |
| (11.5,13.2] | 2936 | 12.48 | 75794.74 |
| (13.2,15] | 2839 | 14.50 | 81603.45 |
| (15,16.8] | 1240 | 16.00 | 93328.99 |
| (16.8,18.5] | 3222 | 17.47 | 115858.88 |
| (18.5,20.2] | 2670 | 19.40 | 136701.51 |
| (20.2,22] | 1615 | 21.45 | 144098.86 |
| (22,23.8] | 584 | 23.00 | 159312.19 |
| (23.8,25.5] | 1116 | 24.46 | 164375.88 |
| (25.5,27.2] | 969 | 26.53 | 195726.17 |
| (27.2,29] | 1284 | 28.50 | 295431.82 |
kable(tabla2, digits = 2, caption = paste0("Parte 2 · k = ", attr(tabla2,"k"),
" clases · amplitud = ", round(attr(tabla2,"amplitud"), 2)))| clase | n_pozos | x_medio | y_medio |
|---|---|---|---|
| [30,31.9] | 1421 | 30.48 | 371985.1 |
| (31.9,33.9] | 1099 | 32.53 | 265403.8 |
| (33.9,35.8] | 1159 | 34.54 | 243750.0 |
| (35.8,37.7] | 975 | 36.52 | 250284.0 |
| (37.7,39.7] | 922 | 38.51 | 201286.9 |
| (39.7,41.6] | 1274 | 40.51 | 166132.3 |
| (41.6,43.5] | 1222 | 42.53 | 160773.2 |
| (43.5,45.5] | 1296 | 44.45 | 166921.6 |
| (45.5,47.4] | 775 | 46.44 | 220442.7 |
| (47.4,49.3] | 592 | 48.44 | 209613.0 |
| (49.3,51.3] | 458 | 50.41 | 196486.0 |
| (51.3,53.2] | 370 | 52.50 | 199175.3 |
| (53.2,55.1] | 915 | 54.78 | 175533.1 |
| (55.1,57.1] | 371 | 56.56 | 291887.7 |
| (57.1,59] | 351 | 58.47 | 254173.2 |
kable(tabla3, digits = 2, caption = paste0("Parte 3 · k = ", attr(tabla3,"k"),
" clases · amplitud = ", round(attr(tabla3,"amplitud"), 2)))| clase | n_pozos | x_medio | y_medio |
|---|---|---|---|
| [60,62.2] | 1049 | 60.67 | 292429.7 |
| (62.2,64.5] | 154 | 63.45 | 299179.7 |
| (64.5,66.7] | 150 | 65.55 | 305899.6 |
| (66.7,68.9] | 197 | 67.63 | 380550.7 |
| (68.9,71.2] | 237 | 70.08 | 324237.0 |
| (71.2,73.4] | 179 | 72.41 | 379722.6 |
| (73.4,75.6] | 150 | 74.45 | 382004.9 |
| (75.6,77.8] | 126 | 76.45 | 383540.5 |
| (77.8,80.1] | 117 | 78.98 | 427609.9 |
| (80.1,82.3] | 90 | 81.54 | 485647.9 |
| (82.3,84.5] | 99 | 83.46 | 506358.8 |
| (84.5,86.8] | 89 | 85.46 | 568473.9 |
| (86.8,89] | 106 | 88.03 | 571222.5 |
5 Gráfica (nube de puntos)
Primero la nube de puntos completa (pozo por pozo), coloreada por parte, y luego cada parte por separado ya con los pares (x̄, ȳ) que se usarán en la regresión.
colores <- c("Parte 1" = "#2E86AB", "Parte 2" = "#E67E22", "Parte 3" = "#C0392B")
plot(datos$YEARS_ACTIVE, datos$CUMULATIVE_PRODUCTION,
col = colores[datos$parte], pch = 16, cex = 0.4,
xlab = "Años activo (X)", ylab = "Producción acumulada (Y)",
main = "Nube de puntos completa - dividida en 3 partes")
abline(v = c(limite1, limite2), lty = 2, col = "gray40")
legend("topright", legend = names(colores), col = colores, pch = 16, bty = "n")
par(mfrow = c(1, 3))
plot(tabla1$x_medio, tabla1$y_medio, pch = 19, col = "#2E86AB",
xlab = "X medio (años)", ylab = "Y medio (producción)", main = "Parte 1")
plot(tabla2$x_medio, tabla2$y_medio, pch = 19, col = "#E67E22",
xlab = "X medio (años)", ylab = "Y medio (producción)", main = "Parte 2")
plot(tabla3$x_medio, tabla3$y_medio, pch = 19, col = "#C0392B",
xlab = "X medio (años)", ylab = "Y medio (producción)", main = "Parte 3")
6 Conjetura (modelo matemático supuesto por parte)
Parte 1 (pozos jóvenes, 0 - 29.7 años): la nube crece rápido al inicio y luego se suaviza — forma típica de curva de declinación de producción (comportamiento hiperbólico/potencial de Arps). Se conjetura un modelo potencial: \(y = a \cdot x^{b}\).
Parte 2 (pozos maduros, 29.7 - 59.3 años): la nube no es monótona: sube, se estabiliza y vuelve a subir (efecto de reacondicionamientos/reactivaciones a mitad de vida del pozo). Ningún modelo monótono (log, exp, potencial) puede seguir esa curvatura, así que se conjetura obligatoriamente un modelo polinómico de grado 2: \(y = a + bx + cx^{2}\).
Parte 3 (pozos muy antiguos, 59.3 - 89 años): el crecimiento se acelera de forma marcada al final (pozos centenarios acumulan de forma compuesta). Se conjetura un modelo exponencial: \(y = a \cdot e^{bx}\).
7 Cálculo de parámetros
7.1 Parte 1 — Modelo potencial
modelo1 <- lm(log(y_medio) ~ log(x_medio), data = tabla1)
a1 <- exp(coef(modelo1)[1])
b1 <- coef(modelo1)[2]
cat("y =", round(a1, 3), "* x^", round(b1, 3))## y = 7404.869 * x^ 0.972##
## Call:
## lm(formula = log(y_medio) ~ log(x_medio), data = tabla1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21980 -0.12876 0.00091 0.04676 0.42970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.90989 0.13400 66.49 < 2e-16 ***
## log(x_medio) 0.97213 0.05143 18.90 2.31e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1662 on 14 degrees of freedom
## Multiple R-squared: 0.9623, Adjusted R-squared: 0.9596
## F-statistic: 357.2 on 1 and 14 DF, p-value: 2.311e-117.2 Parte 2 — Modelo polinómico (grado 2)
modelo2 <- lm(y_medio ~ x_medio + I(x_medio^2), data = tabla2)
a2 <- coef(modelo2)[1]; b2 <- coef(modelo2)[2]; c2 <- coef(modelo2)[3]
cat("y =", round(a2,2), "+", round(b2,2), "*x +", round(c2,4), "*x^2")## y = 1563103 + -60369.54 *x + 656.0594 *x^2##
## Call:
## lm(formula = y_medio ~ x_medio + I(x_medio^2), data = tabla2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49278 -21779 -9006 24140 45994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1563103.5 242185.0 6.454 3.14e-05 ***
## x_medio -60369.5 11158.7 -5.410 0.000157 ***
## I(x_medio^2) 656.1 124.8 5.257 0.000202 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32200 on 12 degrees of freedom
## Multiple R-squared: 0.7252, Adjusted R-squared: 0.6794
## F-statistic: 15.83 on 2 and 12 DF, p-value: 0.0004317.3 Parte 3 — Modelo exponencial
modelo3 <- lm(log(y_medio) ~ x_medio, data = tabla3)
a3 <- exp(coef(modelo3)[1])
b3 <- coef(modelo3)[2]
cat("y =", round(a3, 3), "* e^(", round(b3, 4), "* x )")## y = 58647.76 * e^( 0.0257 * x )##
## Call:
## lm(formula = log(y_medio) ~ x_medio, data = tabla3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.091705 -0.040116 0.006166 0.017453 0.131209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.979305 0.160172 68.55 7.91e-16 ***
## x_medio 0.025710 0.002137 12.03 1.13e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06487 on 11 degrees of freedom
## Multiple R-squared: 0.9294, Adjusted R-squared: 0.9229
## F-statistic: 144.7 on 1 and 11 DF, p-value: 1.133e-078 Comparación del modelo de regresión sobre la realidad
par(mfrow = c(1, 3))
plot(tabla1$x_medio, tabla1$y_medio, col = "#2E86AB", pch = 16,
xlab = "X (años)", ylab = "Y (producción)", main = "Parte 1: Potencial")
curve(a1 * x^b1, add = TRUE, col = "black", lwd = 2)
plot(tabla2$x_medio, tabla2$y_medio, col = "#E67E22", pch = 16,
xlab = "X (años)", ylab = "Y (producción)", main = "Parte 2: Polinómico")
curve(a2 + b2*x + c2*x^2, add = TRUE, col = "black", lwd = 2)
plot(tabla3$x_medio, tabla3$y_medio, col = "#C0392B", pch = 16,
xlab = "X (años)", ylab = "Y (producción)", main = "Parte 3: Exponencial")
curve(a3 * exp(b3 * x), add = TRUE, col = "black", lwd = 2)
9 Test de Pearson y bondad de ajuste
Para los modelos monótonos (potencial y exponencial) se prueba el coeficiente de Pearson sobre la escala en que el modelo es lineal (logarítmica). Para el modelo polinómico, al no ser lineal en una sola escala, se usa el coeficiente de correlación múltiple \(R = \sqrt{R^2}\), que cumple el mismo rol de “bondad de ajuste” entre 0 y 1.
# Parte 1: potencial -> se prueba log(y) vs log(x)
r1 <- cor.test(log(tabla1$x_medio), log(tabla1$y_medio))
# Parte 2: polinómico -> R multiple = sqrt(R^2) del modelo
R2_2 <- summary(modelo2)$r.squared
r2 <- sqrt(R2_2)
# Parte 3: exponencial -> se prueba x vs log(y)
r3 <- cor.test(tabla3$x_medio, log(tabla3$y_medio))
tabla_pearson <- data.frame(
Parte = c("Parte 1 (Potencial)", "Parte 2 (Polinómico)", "Parte 3 (Exponencial)"),
r_o_R = c(r1$estimate, r2, r3$estimate),
R2 = c(r1$estimate^2, R2_2, r3$estimate^2),
Supera_0.7 = c(abs(r1$estimate) > 0.7, r2 > 0.7, abs(r3$estimate) > 0.7),
p_valor_modelo = c(
pf(summary(modelo1)$fstatistic[1], summary(modelo1)$fstatistic[2], summary(modelo1)$fstatistic[3], lower.tail = FALSE),
pf(summary(modelo2)$fstatistic[1], summary(modelo2)$fstatistic[2], summary(modelo2)$fstatistic[3], lower.tail = FALSE),
pf(summary(modelo3)$fstatistic[1], summary(modelo3)$fstatistic[2], summary(modelo3)$fstatistic[3], lower.tail = FALSE)
)
)
kable(tabla_pearson, digits = 4, caption = "Test de bondad de ajuste por parte (umbral |r| o R > 0.7)")| Parte | r_o_R | R2 | Supera_0.7 | p_valor_modelo |
|---|---|---|---|---|
| Parte 1 (Potencial) | 0.9810 | 0.9623 | TRUE | 0e+00 |
| Parte 2 (Polinómico) | 0.8516 | 0.7252 | TRUE | 4e-04 |
| Parte 3 (Exponencial) | 0.9640 | 0.9294 | TRUE | 0e+00 |
Las tres partes superan el umbral de |r| > 0.7 (o
R > 0.7 en el caso del polinómico) y además el
estadístico F de cada modelo es significativo (p-valor < 0.05), por
lo que los tres modelos se aceptan.
10 Estimación
# Un ejemplo de estimación puntual dentro del rango de cada parte
x_est1 <- round(mean(range(tabla1$x_medio)), 1)
x_est2 <- round(mean(range(tabla2$x_medio)), 1)
x_est3 <- round(mean(range(tabla3$x_medio)), 1)
y_est1 <- a1 * x_est1^b1
y_est2 <- a2 + b2*x_est2 + c2*x_est2^2
y_est3 <- a3 * exp(b3 * x_est3)
kable(data.frame(
Parte = c("Parte 1", "Parte 2", "Parte 3"),
X_estimado = c(x_est1, x_est2, x_est3),
Y_estimado = round(c(y_est1, y_est2, y_est3), 1)
), caption = "Estimaciones puntuales de producción acumulada")| Parte | X_estimado | Y_estimado |
|---|---|---|
| Parte 1 | 14.9 | 102331.3 |
| Parte 2 | 44.5 | 175820.7 |
| Parte 3 | 74.3 | 396144.2 |
11 Conclusiones
- 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 fue necesario dividir la nube de puntos en tres partes según el valor máximo de X.
- Parte 1 (pozos jóvenes): relación potencial \(y = 7404.87 \cdot x^{0.972}\), con \(r = 0.981\). Coherente con curvas de declinación de producción típicas de pozos nuevos.
- Parte 2 (pozos maduros): relación polinómica de grado 2 \(y = 1.5631035\times 10^{6} + -6.03695\times 10^{4}x + 656.0594x^2\), con \(R = 0.852\). La curvatura refleja reactivaciones/reacondicionamientos a mitad de vida del pozo, algo que ningún modelo monótono podía capturar.
- Parte 3 (pozos muy antiguos): relación exponencial \(y = 5.864776\times 10^{4} \cdot e^{0.0257x}\), con \(r = 0.964\). La producción acumulada crece de forma compuesta en pozos con décadas de actividad.
- En los tres tramos el coeficiente de correlación (o su equivalente \(R\)) supera 0.7, y el test F de cada modelo resultó significativo, de modo que los tres modelos pasan el test de bondad de ajuste y se aceptan.
- 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.