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

1. Librerías

library(readxl)
library(dplyr)
library(knitr)
library(kableExtra)
library(DT)
library(plotly)

2. Datos

ruta <- file.choose()
datos <- read_excel(ruta)

datos <- datos %>%
  mutate(
    AVG_PRODUCTION         = as.numeric(AVG_PRODUCTION),
    CUMULATIVE_PRODUCTION  = as.numeric(CUMULATIVE_PRODUCTION),
    YEARS_ACTIVE           = as.numeric(YEARS_ACTIVE)
  ) %>%
  na.omit()

nrow(datos)
## [1] 22399

3. Selección de las variables (causa y efecto)

Efecto (Y): AVG_PRODUCTION — Causas (X1, X2): CUMULATIVE_PRODUCTION, YEARS_ACTIVE

Justificación:

  • Precedencia temporal: tanto la producción acumulada como los años activos transcurren u ocurren de forma independiente al promedio de producción; un pozo acumula producción y cumple años de actividad sin que ello dependa de cuál sea su promedio. El promedio de producción, en cambio, es una cantidad derivada — se calcula a partir de lo que el pozo ha producido y del tiempo que lleva activo.
  • Mecanismo físico: el promedio de producción es, por definición, una función de cuánto se ha extraído en total y de cuánto tiempo ha estado el pozo operando (AVG_PRODUCTION ≈ CUMULATIVE_PRODUCTION / tiempo). Es decir, es una consecuencia matemática/física de las otras dos variables, no al revés.
  • Direccionalidad: invertir los roles (usar AVG_PRODUCTION como variable independiente) implicaría que la producción acumulada o los años activos dependen del promedio, lo cual no tiene sentido físico: el promedio no “genera” tiempo de actividad ni producción acumulada.

Por estas razones, CUMULATIVE_PRODUCTION (X1) y YEARS_ACTIVE (X2) se toman como variables independientes, y AVG_PRODUCTION (Y) como variable dependiente.

str(datos[, c("CUMULATIVE_PRODUCTION", "YEARS_ACTIVE", "AVG_PRODUCTION")])
## tibble [22,399 × 3] (S3: tbl_df/tbl/data.frame)
##  $ CUMULATIVE_PRODUCTION: num [1:22399] 47225 275063 82624 7544 681006 ...
##  $ YEARS_ACTIVE         : num [1:22399] 55 55 47 20 28 55 20 48 48 55 ...
##  $ AVG_PRODUCTION       : num [1:22399] 8.59e+09 5.00e+09 1.76e+09 3.77e+02 2.43e+09 ...
##  - attr(*, "na.action")= 'omit' Named int [1:11200] 22400 22401 22402 22403 22404 22405 22406 22407 22408 22409 ...
##   ..- attr(*, "names")= chr [1:11200] "22400" "22401" "22402" "22403" ...
summary(datos[, c("CUMULATIVE_PRODUCTION", "YEARS_ACTIVE", "AVG_PRODUCTION")])
##  CUMULATIVE_PRODUCTION  YEARS_ACTIVE   AVG_PRODUCTION     
##  Min.   :     1        Min.   : 1.00   Min.   :1.000e+00  
##  1st Qu.: 13753        1st Qu.:11.00   1st Qu.:4.562e+06  
##  Median : 54264        Median :20.00   Median :1.678e+09  
##  Mean   :150368        Mean   :25.73   Mean   :2.512e+09  
##  3rd Qu.:184136        3rd Qu.:39.00   3rd Qu.:3.872e+09  
##  Max.   :985283        Max.   :89.00   Max.   :9.999e+09

4. Tabla de estadísticos descriptivos

estadisticos <- data.frame(
  Variable = c("CUMULATIVE_PRODUCTION (X1)", "YEARS_ACTIVE (X2)", "AVG_PRODUCTION (Y)"),
  Media    = c(mean(datos$CUMULATIVE_PRODUCTION), mean(datos$YEARS_ACTIVE), mean(datos$AVG_PRODUCTION)),
  Minimo   = c(min(datos$CUMULATIVE_PRODUCTION), min(datos$YEARS_ACTIVE), min(datos$AVG_PRODUCTION)),
  Maximo   = c(max(datos$CUMULATIVE_PRODUCTION), max(datos$YEARS_ACTIVE), max(datos$AVG_PRODUCTION)),
  Rango    = c(max(datos$CUMULATIVE_PRODUCTION) - min(datos$CUMULATIVE_PRODUCTION),
               max(datos$YEARS_ACTIVE) - min(datos$YEARS_ACTIVE),
               max(datos$AVG_PRODUCTION) - min(datos$AVG_PRODUCTION))
)

kable(estadisticos, digits = 2, format.args = list(big.mark = ","),
      col.names = c("Variable", "Media", "Mínimo", "Máximo", "Rango"),
      align = "lrrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
Variable Media Mínimo Máximo Rango
CUMULATIVE_PRODUCTION (X1) 1.503677e+05 1.00 985,283 985,282
YEARS_ACTIVE (X2) 2.573000e+01 1.00 89 88
AVG_PRODUCTION (Y) 2.512220e+09 0.75 9,999,133,333 9,999,133,332

4.1 Tabla condensada por YEARS_ACTIVE (todos los valores, celdas contraíbles)

condensado <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    n_valores = n(),
    valores_cumulative = paste(CUMULATIVE_PRODUCTION, collapse = ", "),
    valores_avg        = paste(AVG_PRODUCTION, collapse = ", "),
    .groups   = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

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

datatable(condensado_html, escape = FALSE, rownames = FALSE,
          colnames = c("Años Activos", "Producción Acumulada (valores)", "Producción Promedio (valores)"),
          class = "stripe hover compact",
          options = list(pageLength = 10))

4.2 Tabla de tríadas (X1, X2, Y) usando la media aritmética por YEARS_ACTIVE

tabla_media <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    CUMULATIVE_PRODUCTION = mean(CUMULATIVE_PRODUCTION),
    AVG_PRODUCTION        = mean(AVG_PRODUCTION),
    .groups = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

datatable(tabla_media, rownames = FALSE,
          colnames = c("Años Activos", "Media Aritmética de Producción Acumulada", "Media Aritmética de Producción Promedio"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = c("CUMULATIVE_PRODUCTION", "AVG_PRODUCTION"), digits = 2)
write.csv(condensado[, c("YEARS_ACTIVE", "valores_cumulative", "valores_avg")], "tabla_condensada_por_anios.csv", row.names = FALSE)
write.csv(tabla_media, "tabla_media_por_anios.csv", row.names = FALSE)

5. Gráficas (nube de puntos)

plot(datos$CUMULATIVE_PRODUCTION, datos$AVG_PRODUCTION,
     main = "Producción promedio en función de la producción acumulada",
     xlab = "Producción Acumulada (X1)",
     ylab = "Producción Promedio (Y)",
     col  = "#1b2a4a",
     pch  = 19,
     cex  = 1.2)
grid(col = "gray88")

plot(tabla_media$YEARS_ACTIVE, tabla_media$AVG_PRODUCTION,
     main = "Producción promedio (media aritmética) en función de los años activos",
     xlab = "Años Activos (X2)",
     ylab = "Media Aritmética de Producción Promedio (Y)",
     col  = "#5b6b8c",
     pch  = 19,
     cex  = 1.2)
grid(col = "gray88")

5.1 Gráficas en 3D

5.1.1 Nube de puntos 3D — datos originales (sin agrupar)

plot_ly(
  data = datos,
  x = ~CUMULATIVE_PRODUCTION,
  y = ~YEARS_ACTIVE,
  z = ~AVG_PRODUCTION,
  type = "scatter3d",
  mode = "markers",
  marker = list(size = 3, color = "#5b6b8c", opacity = 0.7)
) %>%
  layout(
    title = "Nube de puntos 3D — Datos originales",
    scene = list(
      xaxis = list(title = "Producción Acumulada (X1)"),
      yaxis = list(title = "Años Activos (X2)"),
      zaxis = list(title = "Producción Promedio (Y)")
    )
  )

Esta gráfica usa todos los registros originales, uno por cada pozo, sin ningún tipo de agrupación ni resumen estadístico.

5.1.2 Nube de puntos 3D — con media aritmética por YEARS_ACTIVE

plot_ly(
  data = tabla_media,
  x = ~CUMULATIVE_PRODUCTION,
  y = ~YEARS_ACTIVE,
  z = ~AVG_PRODUCTION,
  type = "scatter3d",
  mode = "markers",
  marker = list(size = 5, color = "#1b2a4a", opacity = 0.9)
) %>%
  layout(
    title = "Nube de puntos 3D — Media Aritmética por Años Activos",
    scene = list(
      xaxis = list(title = "Media Aritmética de Producción Acumulada (X1)"),
      yaxis = list(title = "Años Activos (X2)"),
      zaxis = list(title = "Media Aritmética de Producción Promedio (Y)")
    )
  )

Esta gráfica usa la tabla resumida del punto 4.2, donde cada YEARS_ACTIVE está representado por un solo punto (la media aritmética de CUMULATIVE_PRODUCTION y de AVG_PRODUCTION para ese año). Al reducir el ruido de los datos repetidos, la forma de la relación potencial se aprecia con más claridad.

6. Conjetura

Ambas nubes de puntos muestran un crecimiento que parte cerca de cero y se acelera, sin bajar ni aplanarse. Se conjetura un modelo potencial multivariante:

\[y = a \cdot x_1^{b_1} \cdot x_2^{b_2}\]

7. Cálculo de parámetros (lm en R, regresión múltiple sobre logaritmos)

\[\ln(y) = \ln(a) + b_1 \cdot \ln(x_1) + b_2 \cdot \ln(x_2)\]

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

modelo_log <- lm(log(AVG_PRODUCTION) ~ log(CUMULATIVE_PRODUCTION) + log(YEARS_ACTIVE),
                  data = datos_potencial)
summary(modelo_log)
## 
## Call:
## lm(formula = log(AVG_PRODUCTION) ~ log(CUMULATIVE_PRODUCTION) + 
##     log(YEARS_ACTIVE), data = datos_potencial)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6438  -2.0712   0.8386   3.2026  10.5629 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 8.17420    0.17021   48.02   <2e-16 ***
## log(CUMULATIVE_PRODUCTION)  0.21439    0.02038   10.52   <2e-16 ***
## log(YEARS_ACTIVE)           2.75884    0.03922   70.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.625 on 22396 degrees of freedom
## Multiple R-squared:  0.3183, Adjusted R-squared:  0.3182 
## F-statistic:  5228 on 2 and 22396 DF,  p-value: < 2.2e-16
b1 <- unname(coef(modelo_log)[2])
b2 <- unname(coef(modelo_log)[3])
a  <- unname(exp(coef(modelo_log)[1]))

kable(data.frame(Parámetro = c("a", "b1 (exponente X1)", "b2 (exponente X2)"),
                  Valor = round(c(a, b1, b2), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
Parámetro Valor
a 3548.2265
b1 (exponente X1) 0.2144
b2 (exponente X2) 2.7588

8. Comparación del modelo de regresión sobre la realidad

prediccion <- a * datos_potencial$CUMULATIVE_PRODUCTION^b1 * datos_potencial$YEARS_ACTIVE^b2

plot(datos_potencial$AVG_PRODUCTION, prediccion,
     col  = "#5b6b8c",
     pch  = 19,
     cex  = 1.2,
     main = "Producción promedio real vs. producción promedio estimada por el modelo",
     xlab = "Producción Promedio Real (Y)",
     ylab = "Producción Promedio Estimada (ŷ)")
abline(a = 0, b = 1, col = "#1b2a4a", lwd = 2)
grid(col = "gray88")

legend("topleft",
       legend = c("Pares (real, estimado)", "Ajuste perfecto (y = ŷ)"),
       col    = c("#5b6b8c", "#1b2a4a"),
       pch    = c(19, NA),
       lty    = c(NA, 1),
       lwd    = c(NA, 2),
       bty    = "n")

curva_x1 <- datos_potencial %>%
  arrange(CUMULATIVE_PRODUCTION)

plot(curva_x1$CUMULATIVE_PRODUCTION, curva_x1$AVG_PRODUCTION,
     col  = "#5b6b8c", pch = 19, cex = 1.1,
     main = "Producción promedio en función de la producción acumulada (X2 fijo en su media aritmética)",
     xlab = "Producción Acumulada (X1)",
     ylab = "Producción Promedio (Y)")
grid(col = "gray88")
x2_media <- mean(datos_potencial$YEARS_ACTIVE)
curve(a * x^b1 * x2_media^b2, add = TRUE, col = "#1b2a4a", lwd = 2)

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

9. Test de Pearson y bondad de ajuste

El coeficiente r se calcula por separado para cada causa frente a la producción promedio, sobre los datos originales. El R² mostrado corresponde al modelo de regresión múltiple usado en el paso 7 (modelo_log, sobre log(x1), log(x2) y log(y)), de donde salen los parámetros a, b1 y b2.

r_x1 <- cor(datos_potencial$CUMULATIVE_PRODUCTION, datos_potencial$AVG_PRODUCTION)
r_x2 <- cor(datos_potencial$YEARS_ACTIVE, datos_potencial$AVG_PRODUCTION)
r2   <- summary(modelo_log)$r.squared

resultado_x1 <- ifelse(abs(r_x1) > 0.7, "Aceptado", "Rechazado")
resultado_x2 <- ifelse(abs(r_x2) > 0.7, "Aceptado", "Rechazado")

tabla_test <- data.frame(
  Test      = c("Correlación de Pearson: X1 (CUMULATIVE_PRODUCTION) vs Y",
                "Correlación de Pearson: X2 (YEARS_ACTIVE) vs Y"),
  `r`       = round(c(r_x1, r_x2), 4),
  `|r|`     = round(c(abs(r_x1), abs(r_x2)), 4),
  `R² (modelo múltiple)` = round(rep(r2, 2), 4),
  Criterio  = "|r| > 0.7",
  Resultado = c(resultado_x1, resultado_x2),
  check.names = FALSE
)

kable(tabla_test, align = "lrrrrc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  column_spec(6, bold = TRUE,
              color = ifelse(tabla_test$Resultado == "Aceptado", "#1b2a4a", "#6b6b6b"))
Test r &#124;r&#124; R² (modelo múltiple) Criterio Resultado
Correlación de Pearson: X1 (CUMULATIVE_PRODUCTION) vs Y 0.0608 0.0608 0.3183 &#124;r&#124; > 0.7 Rechazado
Correlación de Pearson: X2 (YEARS_ACTIVE) vs Y 0.2305 0.2305 0.3183 &#124;r&#124; > 0.7 Rechazado

10. Estimación

nuevos_x1 <- c(10000, 50000, 100000, 200000, 300000)
nuevos_x2 <- c(10, 20, 30, 40, 50)

estimaciones <- expand.grid(CUMULATIVE_PRODUCTION = nuevos_x1, YEARS_ACTIVE = nuevos_x2) %>%
  mutate(Produccion_Promedio_Estimada = a * CUMULATIVE_PRODUCTION^b1 * YEARS_ACTIVE^b2) %>%
  arrange(CUMULATIVE_PRODUCTION, YEARS_ACTIVE)

kable(estimaciones, digits = 2, format.args = list(big.mark = ","),
      col.names = c("Producción Acumulada (X1)", "Años Activos (X2)", "Producción Promedio Estimada"),
      align = "lrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
Producción Acumulada (X1) Años Activos (X2) Producción Promedio Estimada
1e+04 10 14,669,430
1e+04 20 99,290,566
1e+04 30 303,889,904
1e+04 40 672,051,779
1e+04 50 1,243,833,722
5e+04 10 20,714,053
5e+04 20 140,203,815
5e+04 30 429,109,486
5e+04 40 948,974,578
5e+04 50 1,756,362,557
1e+05 10 24,032,711
1e+05 20 162,666,270
1e+05 30 497,858,346
1e+05 40 1,101,012,512
1e+05 50 2,037,754,431
2e+05 10 27,883,060
2e+05 20 188,727,499
2e+05 30 577,621,657
2e+05 40 1,277,408,879
2e+05 50 2,364,228,903
3e+05 10 30,415,332
3e+05 20 205,867,269
3e+05 30 630,079,846
3e+05 40 1,393,420,036
3e+05 50 2,578,942,402

11. Conclusiones

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

\[y = 3548.226 \cdot x_1^{0.214} \cdot x_2^{2.759}\]

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

Cuando la producción acumulada es de 150,367.7 y los años activos son 26, se espera una producción promedio de 366,125,094.

La producción promedio está influenciada en un 31.8% conjuntamente por la producción acumulada y los años activos (según el R² del modelo múltiple sobre logaritmos), y en un 68.2% por otros factores.

De forma individual, el test de Pearson dio r = 0.061 (Rechazado) para la producción acumulada, y r = 0.231 (Rechazado) para los años activos.