1. Librerías

library(readxl)
library(dplyr)
library(gt)
library(plotly)
library(DT)

2. Carga de datos

datos <- read_excel("dataset_mundial_petro.xlsx")
cat("Número de registros (filas crudas):", nrow(datos), "\n")
## Número de registros (filas crudas): 49212
cat("Número de variables disponibles:", ncol(datos), "\n")
## Número de variables disponibles: 32

El archivo trae múltiples filas por yacimiento (una fila distinta por cada registro de reservas/producción asociado a un mismo campo), por eso más adelante se agrupa por Unit ID antes de construir la tabla del modelo.

3. Selección de variables

Y = PRODUCTION_START_YEAR (efecto: año en que el yacimiento empieza a producir)

Antes se usaban solo dos variables independientes (Discovery year y FID Year). Se decidió incorporar dos variables adicionales porque resultó más conveniente para tener un modelo más completo, con lógica causal y no solo eligiendo las de mejor correlación:

  • X1 = DISCOVERY_YEAR (Discovery year): año en que se descubrió el yacimiento.
  • X2 = FID_YEAR (FID Year): año de la Decisión Final de Inversión.
  • X3 = STATUS_YEAR (Status year): año del último estado operativo registrado del yacimiento (activo, cerrado, en construcción, etc.).
  • X4 = QUANTITY (Quantity (converted)): volumen de reservas/producción reportado para el yacimiento.

Un yacimiento se descubre (X1), luego recibe aprobación de inversión (X2), y a lo largo de su vida operativa acumula estados y volúmenes reportados (X3, X4) antes o junto con el arranque de producción (Y). Por eso X1, X2, X3 y X4 se tratan como variables explicativas y Y como el efecto.

FID Year y Production start year llegan como texto porque incluyen valores como “2025 (expected)”; se extrae el año numérico de esos casos con una función auxiliar.

extraer_anio <- function(v) as.numeric(sub("^(\\d{4}).*", "\\1", as.character(v)))

datos <- datos %>%
  mutate(
    X1_raw = as.numeric(`Discovery year`),
    X2_raw = extraer_anio(`FID Year`),
    X3_raw = as.numeric(`Status year`),
    X4_raw = as.numeric(`Quantity (converted)`),
    Y_raw  = extraer_anio(`Production start year`)
  )

4. Dominios de las variables

Antes de ajustar cualquier modelo, se define el dominio de cada variable, ya que el modelo solo se interpreta dentro de estos rangos, nunca fuera de ellos.

resumen_dominios <- data.frame(
  Variable  = c("X1 (DISCOVERY_YEAR)", "X2 (FID_YEAR)", "X3 (STATUS_YEAR)",
                "X4 (QUANTITY)", "Y (PRODUCTION_START_YEAR)"),
  Dominio_teorico = c(
    "Enteros positivos (Z+): un año calendario no puede ser negativo, fraccionario ni anterior al origen de la industria petrolera",
    "Enteros positivos (Z+): puede incluir años próximos porque hay proyectos con FID ya aprobado a futuro",
    "Enteros positivos (Z+): fecha de la última actualización de estado del yacimiento",
    "Reales no negativos (R >= 0): un volumen de reservas o producción físico no puede ser negativo",
    "Enteros positivos (Z+): año calendario de arranque de producción; no existen registros de años negativos"
  ),
  Rango_observado = c(
    paste0("[", min(datos$X1_raw, na.rm = TRUE), ", ", max(datos$X1_raw, na.rm = TRUE), "]"),
    paste0("[", min(datos$X2_raw, na.rm = TRUE), ", ", max(datos$X2_raw, na.rm = TRUE), "]"),
    paste0("[", min(datos$X3_raw, na.rm = TRUE), ", ", max(datos$X3_raw, na.rm = TRUE), "]"),
    paste0("[0, ", round(max(datos$X4_raw, na.rm = TRUE), 0), "]"),
    paste0("[", min(datos$Y_raw, na.rm = TRUE), ", ", max(datos$Y_raw, na.rm = TRUE), "]")
  )
) %>%
  gt() %>%
  tab_header(title = md("**Dominio teórico y rango observado de cada variable**"))

resumen_dominios
Dominio teórico y rango observado de cada variable
Variable Dominio_teorico Rango_observado
X1 (DISCOVERY_YEAR) Enteros positivos (Z+): un año calendario no puede ser negativo, fraccionario ni anterior al origen de la industria petrolera [1869, 2023]
X2 (FID_YEAR) Enteros positivos (Z+): puede incluir años próximos porque hay proyectos con FID ya aprobado a futuro [1975, 2029]
X3 (STATUS_YEAR) Enteros positivos (Z+): fecha de la última actualización de estado del yacimiento [1989, 2024]
X4 (QUANTITY) Reales no negativos (R >= 0): un volumen de reservas o producción físico no puede ser negativo [0, 1]
Y (PRODUCTION_START_YEAR) Enteros positivos (Z+): año calendario de arranque de producción; no existen registros de años negativos [1896, 2040]

Sobre Y en particular: como Y representa años, su dominio se restringe a los enteros positivos dentro del rango observado. El modelo solo se usa dentro de ese rango, y ahí Y siempre se mantiene positivo.

5. Tabla de valores (agrupada y balanceada)

En la versión anterior, algunas variables tenían muchos más datos disponibles que otras (por ejemplo, FID Year solo tenía unos pocos cientos de valores frente a miles del resto), lo que generaba tablas desbalanceadas. Esto ocurre porque el archivo trae varias filas por yacimiento (una por cada registro de reservas). Para corregirlo:

  1. Se agrupa por Unit ID (identificador único de yacimiento) y se toma un solo valor representativo de cada variable por yacimiento.
  2. Solo después se filtran los yacimientos que tienen las cinco variables completas a la vez (X1, X2, X3, X4, Y), de modo que todas las filas de la tabla final tengan exactamente la misma cantidad de datos disponibles.
agregado <- datos %>%
  group_by(`Unit ID`) %>%
  summarise(
    DISCOVERY_YEAR        = first(X1_raw),
    FID_YEAR               = first(X2_raw),
    STATUS_YEAR             = first(X3_raw),
    QUANTITY                = suppressWarnings(max(X4_raw, na.rm = TRUE)),
    PRODUCTION_START_YEAR   = first(Y_raw),
    .groups = "drop"
  ) %>%
  mutate(QUANTITY = ifelse(is.infinite(QUANTITY), NA, QUANTITY))

tabla_original <- agregado %>% select(-`Unit ID`)
cat("Yacimientos únicos disponibles:", nrow(tabla_original), "\n")
## Yacimientos únicos disponibles: 8334

5.1 Depuración: solo yacimientos con las 4 variables y Y completas

tabla <- tabla_original %>%
  filter(!is.na(DISCOVERY_YEAR), !is.na(FID_YEAR), !is.na(STATUS_YEAR),
         !is.na(QUANTITY), !is.na(PRODUCTION_START_YEAR)) %>%
  filter(DISCOVERY_YEAR > 0, FID_YEAR > 0, STATUS_YEAR > 0,
         QUANTITY >= 0, PRODUCTION_START_YEAR > 0)

cat("Yacimientos con las 5 variables completas y válidas:", nrow(tabla), "\n")
## Yacimientos con las 5 variables completas y válidas: 248
# Verificación de balance: todas las columnas deben tener el mismo número de datos
sapply(tabla, function(col) sum(!is.na(col)))
##        DISCOVERY_YEAR              FID_YEAR           STATUS_YEAR 
##                   248                   248                   248 
##              QUANTITY PRODUCTION_START_YEAR 
##                   248                   248
datatable(tabla,
          options = list(pageLength = 10),
          colnames = c("Año de Descubrimiento (X1)", "Año de FID (X2)",
                       "Año del Último Estado (X3)", "Cantidad de Reservas/Producción (X4)",
                       "Año de Inicio de Producción (Y)"))

6. Gráfica de dispersión

Con 4 variables independientes ya no se puede ver todo en un solo gráfico 3D. Primero se muestra cada X por separado contra Y (más fácil de leer que una matriz de dispersión completa), y luego un 3D con las dos variables de mayor peso (X1, X2) frente a Y.

QUANTITY tiene valores muy dispares (desde 0 hasta varios millones), así que solo para verla en el gráfico se usa una escala logarítmica; el modelo y los datos originales no se modifican.

tabla <- tabla %>% mutate(QUANTITY_LOG = log10(QUANTITY + 1))

g1 <- plot_ly(tabla, x = ~DISCOVERY_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X1 Discovery"), yaxis = list(title = "Y"))

g2 <- plot_ly(tabla, x = ~FID_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X2 FID"), yaxis = list(title = "Y"))

g3 <- plot_ly(tabla, x = ~STATUS_YEAR, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X3 Status"), yaxis = list(title = "Y"))

g4 <- plot_ly(tabla, x = ~QUANTITY_LOG, y = ~PRODUCTION_START_YEAR,
              type = "scatter", mode = "markers",
              marker = list(size = 5, color = "#2E86C1", opacity = 0.6)) %>%
  layout(xaxis = list(title = "X4 Quantity (log10)"), yaxis = list(title = "Y"))

subplot(g1, g2, g3, g4, nrows = 2, shareX = FALSE, shareY = TRUE,
        titleX = TRUE, titleY = TRUE, margin = 0.06) %>%
  layout(showlegend = FALSE)
plot_ly(tabla, x = ~DISCOVERY_YEAR, y = ~FID_YEAR, z = ~PRODUCTION_START_YEAR,
        type = "scatter3d", mode = "markers",
        marker = list(size = 4, color = "#2E86C1")) %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

7. Conjetura

Los puntos crecen de forma conjunta con X2 principalmente, y X1, X3 y X4 aportan variación adicional sin romper la tendencia general. Se prueba con regresión lineal múltiple de 4 variables:

\[y = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_3 + b_4 x_4\]

8. Cálculo de parámetros

modelo <- lm(PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + STATUS_YEAR + QUANTITY,
             data = tabla)
summary(modelo)
## 
## Call:
## lm(formula = PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + 
##     STATUS_YEAR + QUANTITY, data = tabla)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.4225  -1.2296  -0.0468   1.2834   7.6489 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -804.13349  308.54181  -2.606  0.00972 ** 
## DISCOVERY_YEAR   -0.01017    0.01351  -0.753  0.45208    
## FID_YEAR          0.99450    0.01839  54.079  < 2e-16 ***
## STATUS_YEAR       0.41479    0.15682   2.645  0.00870 ** 
## QUANTITY         -0.31762    1.39034  -0.228  0.81949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.899 on 243 degrees of freedom
## Multiple R-squared:  0.9584, Adjusted R-squared:  0.9577 
## F-statistic:  1400 on 4 and 243 DF,  p-value: < 2.2e-16
b0 <- coef(modelo)[1]
b1 <- coef(modelo)[2]
b2 <- coef(modelo)[3]
b3 <- coef(modelo)[4]
b4 <- coef(modelo)[5]

data.frame(
  Parámetro = c("b0 (intercepto)", "b1 (pendiente X1)", "b2 (pendiente X2)",
                "b3 (pendiente X3)", "b4 (pendiente X4)"),
  Valor     = c(round(b0, 4), round(b1, 4), round(b2, 4), round(b3, 4), round(b4, 6))
) %>%
  gt() %>%
  tab_header(title = md("**Parámetros del Modelo**"))
Parámetros del Modelo
Parámetro Valor
b0 (intercepto) -804.133500
b1 (pendiente X1) -0.010200
b2 (pendiente X2) 0.994500
b3 (pendiente X3) 0.414800
b4 (pendiente X4) -0.317618

Sobre el signo de b1: X1 (DISCOVERY_YEAR) sale con pendiente negativa; el peso explicativo principal lo lleva X2. Esto no es un problema porque, como se definió en la sección 4, el modelo solo se usa dentro del rango de X1 (1938 a 2022), y en ese rango Y se mantiene siempre positivo.

9. Realidad y modelo

Con 4 variables no existe un único “plano” 3D que represente el modelo completo. Se muestran dos evidencias complementarias:

(a) Corte parcial: el plano ajustado para X1 y X2, dejando X3 y X4 fijas en su valor promedio (para poder graficar en 3D).

x3_medio <- mean(tabla$STATUS_YEAR)
x4_medio <- mean(tabla$QUANTITY)

x1_seq <- seq(min(tabla$DISCOVERY_YEAR), max(tabla$DISCOVERY_YEAR), length.out = 25)
x2_seq <- seq(min(tabla$FID_YEAR), max(tabla$FID_YEAR), length.out = 25)
malla  <- outer(x1_seq, x2_seq, function(a, b) b0 + b1 * a + b2 * b + b3 * x3_medio + b4 * x4_medio)

plot_ly() %>%
  add_markers(data = tabla, x = ~DISCOVERY_YEAR, y = ~FID_YEAR, z = ~PRODUCTION_START_YEAR,
              marker = list(size = 4, color = "#2E86C1"), name = "Datos reales") %>%
  add_surface(x = x1_seq, y = x2_seq, z = t(malla), opacity = 0.6,
              showscale = FALSE, name = "Corte del modelo (X3, X4 en su media)") %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

(b) Predicho vs. observado: la forma estándar de evaluar el ajuste cuando hay más de 2 variables independientes.

tabla$Y_predicho <- predict(modelo)

plot_ly(tabla, x = ~PRODUCTION_START_YEAR, y = ~Y_predicho,
        type = "scatter", mode = "markers",
        marker = list(size = 6, color = "#2E86C1")) %>%
  add_trace(x = ~PRODUCTION_START_YEAR, y = ~PRODUCTION_START_YEAR,
             mode = "lines", line = list(color = "red", dash = "dash"),
             name = "Ajuste perfecto") %>%
  layout(title = "Y observado vs. Y predicho",
         xaxis = list(title = "Y observado"),
         yaxis = list(title = "Y predicho"))

Mientras más cerca estén los puntos de la línea roja, mejor es el ajuste del modelo.

10. Test

pearson_x1 <- cor(tabla$DISCOVERY_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x2 <- cor(tabla$FID_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x3 <- cor(tabla$STATUS_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x4 <- cor(tabla$QUANTITY, tabla$PRODUCTION_START_YEAR)
f_pvalue   <- pf(summary(modelo)$fstatistic[1],
                 summary(modelo)$fstatistic[2],
                 summary(modelo)$fstatistic[3], lower.tail = FALSE)
r2         <- summary(modelo)$r.squared

data.frame(
  Test      = c("Pearson X1 (DISCOVERY_YEAR)", "Pearson X2 (FID_YEAR)",
                "Pearson X3 (STATUS_YEAR)", "Pearson X4 (QUANTITY)", "F-test modelo"),
  Valor     = c(round(pearson_x1, 4), round(pearson_x2, 4), round(pearson_x3, 4),
                round(pearson_x4, 4), format(f_pvalue, scientific = TRUE, digits = 4)),
  Resultado = c(ifelse(abs(pearson_x1) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(abs(pearson_x2) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(abs(pearson_x3) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(abs(pearson_x4) > 0.7, "Supera 0.7", "No supera 0.7"),
                ifelse(f_pvalue < 0.05, "Significativo", "No significativo"))
) %>%
  gt() %>%
  tab_header(title = md("**Resumen del Test**"))
Resumen del Test
Test Valor Resultado
Pearson X1 (DISCOVERY_YEAR) 0.6021 No supera 0.7
Pearson X2 (FID_YEAR) 0.9783 Supera 0.7
Pearson X3 (STATUS_YEAR) 0.4043 No supera 0.7
Pearson X4 (QUANTITY) 0.2207 No supera 0.7
F-test modelo 1.946e-166 Significativo
 
cat("R² del modelo:", round(r2, 4), "\n")
## R² del modelo: 0.9584

Nota honesta: de las 4 variables, solo X2 (FID_YEAR) supera 0.7 de correlación con Y por sí sola; X1, X3 y X4 no llegan a ese nivel individualmente. Aun así, el modelo conjunto explica 95.8% de la variabilidad de Y (R² = 0.9584), y el F-test confirma que el modelo en conjunto es significativo.

11. Restricciones

  • El modelo es válido dentro del dominio muestral definido en la sección 4:
    • X1 (DISCOVERY_YEAR): 1938 – 2022
    • X2 (FID_YEAR): 1975 – 2029
    • X3 (STATUS_YEAR): 2018 – 2024
    • X4 (QUANTITY): 0 – 1
    • Y (PRODUCTION_START_YEAR): 1975 – 2034 (siempre entero positivo, como corresponde a un año)
  • No se debe extrapolar el modelo fuera de estos rangos: al ser Y un año calendario, su dominio real son los enteros positivos, y el modelo solo tiene sentido dentro del rango observado.
  • La relación podría dejar de comportarse de forma lineal fuera de este rango (por ejemplo, por cambios regulatorios o tecnológicos entre décadas).
  • Se usaron únicamente yacimientos con las 5 variables completas y reales (sin relleno artificial); no sirve para predecir yacimientos con datos faltantes en X1, X2, X3 o X4.

12. Estimación

escenarios <- data.frame(
  Escenario = c("Típico (medianas)", "Descubrimiento antiguo", "Reservas grandes",
                "FID reciente", "Estado más antiguo"),
  DISCOVERY_YEAR = c(median(tabla$DISCOVERY_YEAR), min(tabla$DISCOVERY_YEAR),
                      median(tabla$DISCOVERY_YEAR), median(tabla$DISCOVERY_YEAR),
                      median(tabla$DISCOVERY_YEAR)),
  FID_YEAR        = c(median(tabla$FID_YEAR), median(tabla$FID_YEAR),
                      median(tabla$FID_YEAR), max(tabla$FID_YEAR),
                      median(tabla$FID_YEAR)),
  STATUS_YEAR     = c(median(tabla$STATUS_YEAR), median(tabla$STATUS_YEAR),
                      median(tabla$STATUS_YEAR), median(tabla$STATUS_YEAR),
                      min(tabla$STATUS_YEAR)),
  QUANTITY        = c(median(tabla$QUANTITY), median(tabla$QUANTITY),
                      max(tabla$QUANTITY), median(tabla$QUANTITY),
                      median(tabla$QUANTITY))
)

escenarios$PRODUCTION_START_YEAR_ESTIMADA <- predict(modelo, newdata = escenarios)

escenarios %>%
  mutate(PRODUCTION_START_YEAR_ESTIMADA = round(PRODUCTION_START_YEAR_ESTIMADA, 2)) %>%
  gt() %>%
  tab_header(title = md("**Estimaciones del Modelo**"))
Estimaciones del Modelo
Escenario DISCOVERY_YEAR FID_YEAR STATUS_YEAR QUANTITY PRODUCTION_START_YEAR_ESTIMADA
Típico (medianas) 2004 2018.5 2022 1 2021.28
Descubrimiento antiguo 1938 2018.5 2022 1 2021.95
Reservas grandes 2004 2018.5 2022 1 2021.28
FID reciente 2004 2029.0 2022 1 2031.72
Estado más antiguo 2004 2018.5 2018 1 2019.62

13. Conclusión

Con DISCOVERY_YEAR (X1), FID_YEAR (X2), STATUS_YEAR (X3) y QUANTITY (X4) se obtiene el modelo de regresión lineal múltiple:

\[y = -804.13 + -0.0102x_1 + 0.9945x_2 + 0.4148x_3 + -0.317618x_4\]

Interpretación de cada pendiente (manteniendo las demás variables constantes):

  • Por cada año adicional en X1 (descubrimiento), Y cambia en -0.0102.
  • Por cada año adicional en X2 (FID), Y cambia en 0.9945.
  • Por cada año adicional en X3 (último estado), Y cambia en 0.4148.
  • Por cada unidad adicional de X4 (reservas/producción), Y cambia en -0.317618.

El modelo explica 95.8% de la variabilidad de Y (R² = 0.9584). El F-test resulta significativo (p = 1.95e-166), lo que confirma que, en conjunto, las 4 variables sí aportan información útil para explicar Y — aunque, tomadas de forma individual, X1, X3 y X4 no superen 0.7 de correlación con Y.

El dominio de Y se definió como los enteros positivos dentro del rango observado (1975–2034), y el modelo se interpreta únicamente dentro de ese dominio.