0.Carga de librerías

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

1.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.

2.Selección de tres variables (Causa y Efecto)

Causa y Efecto:

Un yacimiento se descubre (X1), luego recibe aprobación de inversión (X2), y a lo largo de su vida operativa acumula un estado registrado (X3) antes o junto con el arranque de producción (Y). Por eso X1, X2 y X3 se tratan como variables explicativas (causa) y Y como el efecto:

Y = PRODUCTION_START_YEAR: año en que el yacimiento empieza a producir.
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.).

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`),
    Y_raw  = extraer_anio(`Production start year`)
  )

3.Tabla de Cuarteta de Valores

En la versión original, 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, primero se agrupa por Unit ID (identificador único de yacimiento) tomando un solo valor representativo por variable, y luego se filtran solo los yacimientos que tienen las cuatro variables completas a la vez (X1, X2, X3, 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),
    PRODUCTION_START_YEAR   = first(Y_raw),
    .groups = "drop"
  )

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

## Yacimientos únicos disponibles: 8334

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

cat("Yacimientos con las 4 variables completas y válidas:", nrow(tabla), "\n")

## Yacimientos con las 4 variables completas y válidas: 259

# 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 
##                   259                   259                   259 
## PRODUCTION_START_YEAR 
##                   259

datatable(tabla,
          options = list(pageLength = 10),
          colnames = c("Año de Descubrimiento (X1)", "Año de FID (X2)",
                       "Año del Último Estado (X3)",
                       "Año de Inicio de Producción (Y)"))

4.Gráfica de Dispersión

Primero se muestra cada X por separado contra Y (más fácil de leer una por una), y luego un 3D con las tres variables juntas.

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"))

subplot(g1, g2, g3, nrows = 1, shareX = FALSE, shareY = TRUE,
        titleX = TRUE, titleY = TRUE, margin = 0.05) %>%
  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 = ~STATUS_YEAR, colorscale = "Viridis", showscale = TRUE)) %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

El color de los puntos en el 3D representa a X3 (STATUS_YEAR), la tercera variable que no cabe en los ejes.

5.Conjetura

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

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

6.Parámetros

Modelo matemático:

\[Y = b_0 + b_1 X_1 + b_2 X_2 + b_3 X_3\]

Ajuste del modelo:

modelo <- lm(PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + STATUS_YEAR,
             data = tabla)
summary(modelo)

## 
## Call:
## lm(formula = PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR + 
##     STATUS_YEAR, data = tabla)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.474  -1.264   0.004   1.307   7.709 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -709.74418  286.59299  -2.476   0.0139 *  
## DISCOVERY_YEAR   -0.01276    0.01293  -0.986   0.3248    
## FID_YEAR          0.99566    0.01774  56.126   <2e-16 ***
## STATUS_YEAR       0.36935    0.14620   2.526   0.0121 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.855 on 255 degrees of freedom
## Multiple R-squared:  0.9584, Adjusted R-squared:  0.9579 
## F-statistic:  1960 on 3 and 255 DF,  p-value: < 2.2e-16

Parámetros:

b0 <- coef(modelo)[1]
b1 <- coef(modelo)[2]
b2 <- coef(modelo)[3]
b3 <- coef(modelo)[4]

data.frame(
  Parámetro = c("b0 (intercepto)", "b1 (pendiente X1)", "b2 (pendiente X2)",
                "b3 (pendiente X3)"),
  Valor     = c(round(b0, 4), round(b1, 4), round(b2, 4), round(b3, 4))
) %>%
  gt() %>%
  tab_header(title = md("**Parámetros del Modelo**"))

Parámetro	Valor
Parámetros del Modelo
b0 (intercepto)	-709.7442
b1 (pendiente X1)	-0.0128
b2 (pendiente X2)	0.9957
b3 (pendiente X3)	0.3693

Ecuación múltiple:

\[y = -709.74 + -0.0128x_1 + 0.9957x_2 + 0.3693x_3\]

Sobre el signo de b1: X1 (DISCOVERY_YEAR) sale con pendiente negativa; el peso explicativo principal lo lleva X2. Esto se revisa formalmente en la Sección 9 (Restricciones).

7.Gráfica de comparación de la realidad con el modelo

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)
x3_medio <- mean(tabla$STATUS_YEAR)
malla  <- outer(x1_seq, x2_seq, function(a, b) b0 + b1 * a + b2 * b + b3 * x3_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 = "Modelo (X3 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")
  ))

Como X3 también forma parte del modelo, la superficie de arriba se calcula dejando X3 fija en su promedio para poder graficarla en 3D. Para ver el efecto de las 3 variables a la vez, se compara además lo predicho contra lo observado:

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.

8.Test de Bondad

8.1 Correlación de Pearson

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)

data.frame(
  Variable = c("X1 (DISCOVERY_YEAR)", "X2 (FID_YEAR)", "X3 (STATUS_YEAR)"),
  Pearson  = c(round(pearson_x1, 4), round(pearson_x2, 4), round(pearson_x3, 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"))
) %>%
  gt() %>%
  tab_header(title = md("**Correlación de Pearson por variable**"))

Variable	Pearson	Resultado
Correlación de Pearson por variable
X1 (DISCOVERY_YEAR)	0.5950	No supera 0.7
X2 (FID_YEAR)	0.9783	Supera 0.7
X3 (STATUS_YEAR)	0.3990	No supera 0.7

Nota honesta: de las 3 variables, solo X2 (FID_YEAR) supera 0.7 de correlación con Y por sí sola; X1 y X3 no llegan a ese nivel individualmente.

8.2 Coeficiente de Determinación

r2       <- summary(modelo)$r.squared
f_pvalue <- pf(summary(modelo)$fstatistic[1],
               summary(modelo)$fstatistic[2],
               summary(modelo)$fstatistic[3], lower.tail = FALSE)

cat("R² del modelo:", round(r2, 4), "\n")

## R² del modelo: 0.9584

cat("F-test (p-valor):", format(f_pvalue, scientific = TRUE, digits = 4), "\n")

## F-test (p-valor): 9.661e-176

cat("¿F-test significativo (p < 0.05)?:", f_pvalue < 0.05, "\n")

## ¿F-test significativo (p < 0.05)?: TRUE

Aun cuando X1 y X3 no superan 0.7 de correlación por separado, 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.

9.Restricciones

Dominios:

resumen_dominios <- data.frame(
  Variable  = c("X1 (DISCOVERY_YEAR)", "X2 (FID_YEAR)", "X3 (STATUS_YEAR)",
                "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",
    "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("[", 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

Variable	Dominio_teorico	Rango_observado
Dominio teórico y rango observado de cada variable
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]
Y (PRODUCTION_START_YEAR)	Enteros positivos (Z+): año calendario de arranque de producción; no existen registros de años negativos	[1896, 2040]

Dominio [X1] (DISCOVERY_YEAR): D = {Z+ ∩ [1938, 2022]}

Dominio [X2] (FID_YEAR): D = {Z+ ∩ [1975, 2029]}

Dominio [X3] (STATUS_YEAR): D = {Z+ ∩ [2018, 2024]}

Dominio [Y] (PRODUCTION_START_YEAR): D = {Z+ ∩ [1975, 2034]}

¿Existe una combinación de X1, X2 y X3 que, reemplazada en el plano, genere un valor de Y fuera de su dominio (años negativos, cero o no enteros)?

Aunque b1 es negativo, esto no implica automáticamente una restricción: el modelo solo se interpreta dentro del dominio muestral. Se verifica evaluando la ecuación del modelo en los 8 vértices (combinaciones mín/máx de X1, X2 y X3):

combinaciones_extremas <- expand.grid(
  DISCOVERY_YEAR = c(min(tabla$DISCOVERY_YEAR), max(tabla$DISCOVERY_YEAR)),
  FID_YEAR        = c(min(tabla$FID_YEAR), max(tabla$FID_YEAR)),
  STATUS_YEAR     = c(min(tabla$STATUS_YEAR), max(tabla$STATUS_YEAR))
)
combinaciones_extremas$Y_estimado <- predict(modelo, newdata = combinaciones_extremas)
combinaciones_extremas <- combinaciones_extremas %>%
  mutate(Y_estimado = round(Y_estimado, 2),
         Fuera_de_dominio = Y_estimado <= 0)

combinaciones_extremas %>%
  gt() %>%
  tab_header(title = md("**Verificación de restricciones en los vértices del dominio**"))

DISCOVERY_YEAR	FID_YEAR	STATUS_YEAR	Y_estimado	Fuera_de_dominio
Verificación de restricciones en los vértices del dominio
1938	1975	2018	1977.31	FALSE
2022	1975	2018	1976.24	FALSE
1938	2029	2018	2031.08	FALSE
2022	2029	2018	2030.00	FALSE
1938	1975	2024	1979.53	FALSE
2022	1975	2024	1978.45	FALSE
1938	2029	2024	2033.29	FALSE
2022	2029	2024	2032.22	FALSE

hay_restriccion <- any(combinaciones_extremas$Fuera_de_dominio)
cat("¿Alguna combinación genera un Y fuera de dominio?:", hay_restriccion, "\n")

## ¿Alguna combinación genera un Y fuera de dominio?: FALSE

No hay restricciones: ninguna combinación de X1, X2 y X3 dentro del dominio muestral genera un valor de Y fuera de su dominio (Y siempre se mantiene positivo).

10.Estimación

escenarios <- data.frame(
  Escenario = c("Típico (medianas)", "Descubrimiento antiguo", "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)),
  FID_YEAR        = c(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), min(tabla$STATUS_YEAR))
)

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**"))

Escenario	DISCOVERY_YEAR	FID_YEAR	STATUS_YEAR	PRODUCTION_START_YEAR_ESTIMADA
Estimaciones del Modelo
Típico (medianas)	2005	2019	2022	2021.74
Descubrimiento antiguo	1938	2019	2022	2022.60
FID reciente	2005	2029	2022	2031.70
Estado más antiguo	2005	2019	2018	2020.26

11.Conclusión

Entre PRODUCTION_START_YEAR (Y) y DISCOVERY_YEAR (X1), FID_YEAR (X2) y STATUS_YEAR (X3) (medidos en años) existe una relación de tipo lineal definida por la ecuación del plano