library(readr) # lectura del csv
library(dplyr) # manipulación de datos
library(ggplot2) # gráficos
library(gt) # tablas con diseño premium
library(stringr) # manejo de texto
library(DT) # tabla interactiva paginada
cat("Librerías cargadas: readr, dplyr, ggplot2, gt, stringr, DT")Librerías cargadas: readr, dplyr, ggplot2, gt, stringr, DT
## DATASET ##
# Archivo con las variables YEARS_ACTIVE y CUMULATIVE_PRODUCTION.
# IMPORTANTE: ajusta esta ruta a la ubicación real del archivo en tu equipo.
ruta_archivo <- "oil_and_gas_leases_data___2_.csv"
datos <- read_csv(ruta_archivo, 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)
## Estructura de los datos
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 ...
Se definieron los Años Activos
(YEARS_ACTIVE) como variable independiente / causa
(x), ya que representan el tiempo durante el cual un pozo ha
estado en operación, constituyendo un factor clave en la acumulación de
producción.
La Producción Acumulada
(CUMULATIVE_PRODUCTION) actúa como variable
dependiente / efecto (y), debido a que refleja la
cantidad total de hidrocarburos extraídos a lo largo del tiempo de vida
del pozo.
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 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.
Esta relación busca modelar el comportamiento de la producción en función del tiempo. Desde el punto de vista operativo, a medida que aumentan los años activos, la producción acumulada tiende a incrementarse, aunque no necesariamente de forma lineal, debido a factores como declinación de producción, mantenimiento y condiciones del yacimiento.
El dataset depurado contiene 47757 pozos, pero solo 89 valores
distintos de Años Activos (x): al haber miles de pozos
que comparten el mismo año de actividad, trabajar con los datos crudos
implica manejar valores de x fuertemente repetidos. Por ello, tanto para
la variable independiente x (Años Activos) como para la
dependiente y (Producción Acumulada) se calcula su
media por cada año (floor(YEARS_ACTIVE)),
obteniendo así un único par (x̄, ȳ) por año — sin valores repetidos — que
es el que se emplea en el resto del análisis.
tabla_xy_completa <- datos %>%
mutate(Año = floor(YEARS_ACTIVE)) %>%
group_by(Año) %>%
summarise(
n_pozos = n(),
x_medio = mean(YEARS_ACTIVE),
y_medio = mean(CUMULATIVE_PRODUCTION),
.groups = "drop"
) %>%
arrange(Año)
cat("Total de pares (x̄, ȳ) obtenidos, uno por año:", nrow(tabla_xy_completa), "\n")Total de pares (x̄, ȳ) obtenidos, uno por año: 89
# Tabla condensada: para cada Año Activo se listan, en una celda
# contraíble (expandir/colapsar), todos los valores individuales de
# Producción Acumulada de los pozos que comparten ese año.
tabla_condensada <- datos %>%
mutate(`Años Activos` = floor(YEARS_ACTIVE)) %>%
group_by(`Años Activos`) %>%
summarise(
`Producción Acumulada (valores)` = sprintf(
'<details><summary>%d valor(es)</summary><div style="max-height:180px; overflow-y:auto; padding:6px 4px; font-size:12px; line-height:1.6;">%s</div></details>',
n(),
paste(format(round(CUMULATIVE_PRODUCTION, 2), big.mark = ","), collapse = ", ")
),
.groups = "drop"
) %>%
arrange(`Años Activos`)
datatable(
tabla_condensada,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.1: Tabla condensada (todos los valores, celdas contraíbles)"
),
escape = FALSE,
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = 0))
)
)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 = round(c(min_x, min_y), 2),
Maximo = round(c(max_x, max_y), 2),
Rango = round(c(max_x - min_x, max_y - min_y), 2),
Media = round(c(mean(datos$YEARS_ACTIVE), mean(datos$CUMULATIVE_PRODUCTION)), 2)
)
datatable(
resumen_general,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.2: Resumen General de las Variables"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
)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")
tabla_distribucion <- 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 = c(nrow(parte1), nrow(parte2), nrow(parte3))
)
datatable(
tabla_distribucion,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.3: Distribución de Pozos por Parte"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
)Con miles de pozos por parte, la nube de puntos es puro ruido. Por
eso, en lugar de agrupar por clases de rango, se agrupan los pozos
por cada año exacto de actividad
(floor(YEARS_ACTIVE)) y se calcula el promedio de
ambas variables —el par (x̄, ȳ)— para cada año. Estas parejas,
una por cada año representado dentro de la parte, son las que se usan
para ajustar los modelos de cada parte.
# Las tres partes se obtienen ahora dividiendo directamente la tabla
# resumida (tabla_xy_completa), es decir, sobre los pares (x̄, ȳ) ya
# promediados por año, y no sobre los pozos crudos.
tabla_xy_completa <- tabla_xy_completa %>%
mutate(parte = case_when(
x_medio <= limite1 ~ "Parte 1",
x_medio <= limite2 ~ "Parte 2",
TRUE ~ "Parte 3"
))
tabla1 <- tabla_xy_completa %>% filter(parte == "Parte 1") %>% select(-parte)
tabla2 <- tabla_xy_completa %>% filter(parte == "Parte 2") %>% select(-parte)
tabla3 <- tabla_xy_completa %>% filter(parte == "Parte 3") %>% select(-parte)datatable(
tabla1,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.4.1: Pares (x̄, ȳ) — Parte 1 (", nrow(tabla1), " años promediados)"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) %>%
formatRound(columns = c("x_medio", "y_medio"), digits = 2)datatable(
tabla2,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.4.2: Pares (x̄, ȳ) — Parte 2 (", nrow(tabla2), " años promediados)"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) %>%
formatRound(columns = c("x_medio", "y_medio"), digits = 2)datatable(
tabla3,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 4.4.3: Pares (x̄, ȳ) — Parte 3 (", nrow(tabla3), " años promediados)"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) %>%
formatRound(columns = c("x_medio", "y_medio"), digits = 2)Primero la nube de puntos completa, construida a partir de la tabla resumida (pares x̄, ȳ promediados por año), coloreada por parte.
colores <- c("Parte 1" = "#2E86AB", "Parte 2" = "#E67E22", "Parte 3" = "#C0392B")
par(mar = c(5, 5, 4, 9))
plot(tabla_xy_completa$x_medio, tabla_xy_completa$y_medio,
col = colores[tabla_xy_completa$parte], pch = 16, cex = 1,
xlab = "Años Activos (X)", ylab = "Producción Acumulada (Y)",
main = "Gráfica N°1: Nube de puntos (tabla resumida), dividida en 3 partes",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
abline(v = c(limite1, limite2), lty = 2, col = "gray40")
legend("topright", inset = c(-0.32, 0), legend = names(colores), col = colores,
pch = 16, bty = "n", xpd = TRUE)
box()Ya con los pares (x̄, ȳ) que se usarán en la regresión.
par(mar = c(5, 5, 4, 2))
plot(tabla1$x_medio, tabla1$y_medio, pch = 19, col = "#2E86AB",
xlab = "X medio (Años Activos)", ylab = "Y medio (Producción Acumulada)",
main = "Gráfica N°2: Nube de Puntos — Parte 1",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()par(mar = c(5, 5, 4, 2))
plot(tabla2$x_medio, tabla2$y_medio, pch = 19, col = "#E67E22",
xlab = "X medio (Años Activos)", ylab = "Y medio (Producción Acumulada)",
main = "Gráfica N°3: Nube de Puntos — Parte 2",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()par(mar = c(5, 5, 4, 2))
plot(tabla3$x_medio, tabla3$y_medio, pch = 19, col = "#C0392B",
xlab = "X medio (Años Activos)", ylab = "Y medio (Producción Acumulada)",
main = "Gráfica N°4: Nube de Puntos — Parte 3",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()Parte 1 (pozos jóvenes, 0 - 30 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, 30 - 59 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 - 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}\]
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 = 7490.256 * x^ 0.972
Call:
lm(formula = log(y_medio) ~ log(x_medio), data = tabla1)
Residuals:
Min 1Q Median 3Q Max
-0.26168 -0.06562 -0.01314 0.02253 0.43217
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.92136 0.10032 88.93 <0.0000000000000002 ***
log(x_medio) 0.97196 0.03867 25.14 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1734 on 27 degrees of freedom
Multiple R-squared: 0.959, Adjusted R-squared: 0.9575
F-statistic: 631.9 on 1 and 27 DF, p-value: < 0.00000000000000022
El modelo potencial (\(y = a \cdot x^{b}\)) queda definido por dos parámetros, a y b:
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 = 1557110 + -60147.85 *x + 654.2263 *x^2
Call:
lm(formula = y_medio ~ x_medio + I(x_medio^2), data = tabla2)
Residuals:
Min 1Q Median 3Q Max
-60489 -18498 -10897 17185 66573
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1557110.1 166220.4 9.368 0.000000000566 ***
x_medio -60147.8 7665.0 -7.847 0.000000019474 ***
I(x_medio^2) 654.2 85.8 7.625 0.000000033495 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31440 on 27 degrees of freedom
Multiple R-squared: 0.7113, Adjusted R-squared: 0.69
F-statistic: 33.27 on 2 and 27 DF, p-value: 0.00000005189
El modelo polinómico de grado 2 (\(y = a + bx + cx^{2}\)) queda definido por tres parámetros, a, b y c:
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 = 56824.06 * e^( 0.0259 * x )
Call:
lm(formula = log(y_medio) ~ x_medio, data = tabla3)
Residuals:
Min 1Q Median 3Q Max
-0.22259 -0.05647 -0.01596 0.06611 0.23948
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.947715 0.151249 72.38 < 0.0000000000000002 ***
x_medio 0.025944 0.002017 12.87 0.000000000000283 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0956 on 28 degrees of freedom
Multiple R-squared: 0.8553, Adjusted R-squared: 0.8501
F-statistic: 165.5 on 1 and 28 DF, p-value: 0.0000000000002833
El modelo exponencial (\(y = a \cdot e^{bx}\)) queda definido por dos parámetros, a y b:
par(mar = c(5, 5, 4, 9))
plot(tabla1$x_medio, tabla1$y_medio, col = "#2E86AB", pch = 16,
xlab = "X (Años Activos)", ylab = "Y (Producción Acumulada)",
main = "Gráfica N°5: Modelo Potencial vs. Realidad — Parte 1",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
curve(a1 * x^b1, add = TRUE, col = "#1F2A33", lwd = 3)
legend("topright", inset = c(-0.32, 0), legend = c("Datos (x̄, ȳ)", "Modelo Potencial"),
col = c("#2E86AB", "#1F2A33"), pch = c(16, NA), lwd = c(NA, 3), bty = "n", xpd = TRUE)
box()par(mar = c(5, 5, 4, 9))
plot(tabla2$x_medio, tabla2$y_medio, col = "#E67E22", pch = 16,
xlab = "X (Años Activos)", ylab = "Y (Producción Acumulada)",
main = "Gráfica N°6: Modelo Polinómico vs. Realidad — Parte 2",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
curve(a2 + b2*x + c2*x^2, add = TRUE, col = "#1F2A33", lwd = 3)
legend("topright", inset = c(-0.32, 0), legend = c("Datos (x̄, ȳ)", "Modelo Polinómico"),
col = c("#E67E22", "#1F2A33"), pch = c(16, NA), lwd = c(NA, 3), bty = "n", xpd = TRUE)
box()par(mar = c(5, 5, 4, 9))
plot(tabla3$x_medio, tabla3$y_medio, col = "#C0392B", pch = 16,
xlab = "X (Años Activos)", ylab = "Y (Producción Acumulada)",
main = "Gráfica N°7: Modelo Exponencial vs. Realidad — Parte 3",
cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
curve(a3 * exp(b3 * x), add = TRUE, col = "#1F2A33", lwd = 3)
legend("topright", inset = c(-0.32, 0), legend = c("Datos (x̄, ȳ)", "Modelo Exponencial"),
col = c("#C0392B", "#1F2A33"), pch = c(16, NA), lwd = c(NA, 3), bty = "n", xpd = TRUE)
box()El coeficiente de correlación de Pearson, calculado en R mediante
cor(x, y), debe superar 0.7 o -0.7 para
aceptar el modelo, dado que toma valores entre 1 y -1. Para los modelos
monótonos (potencial y exponencial) se prueba dicho coeficiente 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; en este caso el coeficiente
de determinación (\(R^2\))
aparece directamente en el summary() del modelo lineal.
# 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)
)
)
datatable(
tabla_pearson,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 9: Test de Bondad de Ajuste por Parte (umbral |r| o R > 0.7)"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) %>%
formatRound(columns = c("r_o_R", "R2", "p_valor_modelo"), digits = 4)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.
# 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)
tabla_estimaciones <- 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)
)
datatable(
tabla_estimaciones,
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left; font-size: 16px; font-weight: 700; color:#1F2A33;",
"Tabla 10: Estimaciones Puntuales de Producción Acumulada"
),
rownames = FALSE,
class = "display compact stripe hover",
options = list(
pageLength = 10,
dom = "ltip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
)