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:", nrow(datos), "\n")
## Número de registros: 8334
cat("Número de variables:", ncol(datos), "\n")
## Número de variables: 23

3.Selección de variables

Y = PRODUCTION_START_YEAR (efecto) X1 = DISCOVERY_YEAR, X2 = FID_YEAR (causas)

Un yacimiento comienza a producir después de haber sido descubierto (X1) y después de que se tome la Decisión Final de Inversión / FID (X2). Por eso X1 y X2 son las causas, y el año de inicio de producción (Y) es el efecto.

FID Year y Production start year vienen como texto porque incluyen valores como “2025 (expected)”; se extrae el año numérico de esos casos.

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

X1 <- as.numeric(datos$`Discovery year`)
X2 <- extraer_anio(datos$`FID Year`)
Y  <- extraer_anio(datos$`Production start year`)

tabla_original <- data.frame(DISCOVERY_YEAR = X1, FID_YEAR = X2, PRODUCTION_START_YEAR = Y)

4.Tabla de pares de valores

4.1 Tripleta de valores originales

nrow(tabla_original)
## [1] 8334
datatable(head(tabla_original, 50),
          options = list(pageLength = 10),
          colnames = c("Año de Descubrimiento (X1)", "Año de FID (X2)", "Año de Inicio de Producción (Y)"))

4.2 Tripleta depurada (solo registros completos)

Se conservan únicamente los registros donde X1, X2 y Y están presentes y son positivos, ya que faltantes o valores inventados distorsionarían la relación real entre las variables.

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

cat("Tripletas completas y válidas:", nrow(tabla), "\n")
## Tripletas completas y válidas: 259
datatable(tabla,
          options = list(pageLength = 10),
          colnames = c("Año de Descubrimiento (X1)", "Año de FID (X2)", "Año de Inicio de Producción (Y)"))

5.Gráfica de dispersión

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

6.Conjetura

Los puntos se ven creciendo de forma conjunta con X1 y X2, formando una nube más o menos plana (no caótica). Se prueba con regresión lineal múltiple (un plano):

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

7.Cálculo de parámetros

modelo <- lm(PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR, data = tabla)
summary(modelo)
## 
## Call:
## lm(formula = PRODUCTION_START_YEAR ~ DISCOVERY_YEAR + FID_YEAR, 
##     data = tabla)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.9703  -1.5127   0.1503   1.3026   8.2334 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    11.13115   26.98415   0.413    0.680    
## DISCOVERY_YEAR -0.01555    0.01302  -1.195    0.233    
## FID_YEAR        1.01132    0.01680  60.211   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.885 on 256 degrees of freedom
## Multiple R-squared:  0.9574, Adjusted R-squared:  0.9571 
## F-statistic:  2876 on 2 and 256 DF,  p-value: < 2.2e-16
b0 <- coef(modelo)[1]
b1 <- coef(modelo)[2]
b2 <- coef(modelo)[3]

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

8.Realidad y 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)
malla  <- outer(x1_seq, x2_seq, function(a, b) b0 + b1 * a + b2 * b)

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 = "Plano ajustado") %>%
  layout(scene = list(
    xaxis = list(title = "X1 DISCOVERY_YEAR"),
    yaxis = list(title = "X2 FID_YEAR"),
    zaxis = list(title = "Y PRODUCTION_START_YEAR")
  ))

Los puntos quedan cerca del plano; se evalúa el ajuste formalmente en la sección siguiente.

9.Test

pearson_x1 <- cor(tabla$DISCOVERY_YEAR, tabla$PRODUCTION_START_YEAR)
pearson_x2 <- cor(tabla$FID_YEAR, 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)", "F-test modelo"),
  Valor     = c(round(pearson_x1, 4), round(pearson_x2, 4), format(f_pvalue, scientific = TRUE, digits = 4)),
  Resultado = c(ifelse(abs(pearson_x1) > 0.7, "Aceptado", "No supera 0.7"),
                ifelse(abs(pearson_x2) > 0.7, "Aceptado", "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.595 No supera 0.7
Pearson X2 (FID_YEAR) 0.9783 Aceptado
F-test modelo 3.726e-176 Significativo
 
cat("R² del modelo:", round(r2, 4), "\n")
## R² del modelo: 0.9574

Nota honesta: X2 (FID_YEAR) pasa ampliamente el test de Pearson (0.9783), pero X1 (DISCOVERY_YEAR) individualmente no supera 0.7 (da 0.595). Aun así, el modelo conjunto explica 95.7% de la variabilidad de Y (R² = 0.9574), porque X2 aporta la mayor parte de la explicación y X1 suma información adicional. Si tu profesor exige que ambas variables superen 0.7 individualmente, este par no cumple esa condición estricta — avisame y probamos otra combinación.

10.Restricciones

  • Válido para:
    • X1 (DISCOVERY_YEAR) desde 1938 hasta infinito positivo
    • X2 (FID_YEAR) desde 1975 hasta infinito positivo
  • La inferencia de la muestra a la población (yacimientos reales) se respalda en el R² dentro del rango observado, pero entre más se aleje la predicción de ese rango, menor la confiabilidad de la estimación, ya que la relación podría dejar de comportarse de forma lineal a medida que cambian las condiciones (por ejemplo, cambios regulatorios o tecnológicos entre décadas).
  • Se usaron únicamente registros completos y reales (sin relleno artificial); no sirve para predecir yacimientos con datos faltantes en X1 o X2.

11.Estimación

escenarios <- expand.grid(
  DISCOVERY_YEAR = c(min(tabla$DISCOVERY_YEAR), 2000, 2015),
  FID_YEAR        = c(min(tabla$FID_YEAR), 2005, 2020)
)
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
DISCOVERY_YEAR FID_YEAR PRODUCTION_START_YEAR_ESTIMADA
1938 1975 1978.34
2000 1975 1977.38
2015 1975 1977.14
1938 2005 2008.68
2000 2005 2007.72
2015 2005 2007.48
1938 2020 2023.85
2000 2020 2022.89
2015 2020 2022.65

12.Conclusión

Con DISCOVERY_YEAR y FID_YEAR se obtiene el mejor ajuste con regresión lineal múltiple (un plano en 3D):

\[y = 11.13 + -0.0156x_1 + 1.0113x_2\]

Por cada unidad que aumenta el año de descubrimiento, Y cambia en -0.0156 (con FID_YEAR constante); por cada unidad que aumenta el año de FID, Y cambia en 1.0113 (con DISCOVERY_YEAR constante).

El modelo explica 95.7% de la variabilidad de Y (R²). Pearson: r = 0.595 para X1 (no supera 0.7), r = 0.978 para X2 (aceptado). F-test: significativo (p = 3.73e-176).