Modelo de Regresión Potencial


1.Carga de librerías

library(readxl)
library(dplyr)
library(gt)

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

La variable Discovery year (año de descubrimiento) actúa como variable independiente o causa (X), ya que el año en que se descubre un yacimiento antecede necesariamente a su puesta en producción. La variable Production start year (año de inicio de producción) actúa como variable dependiente o efecto (Y), ya que depende del tiempo de desarrollo posterior al descubrimiento.

Ambas variables son siempre positivas (años calendario), condición necesaria para un modelo potencial \(y = a \cdot x^{b}\), ya que requiere calcular \(\ln(x)\) y \(\ln(y)\).

Production start year viene como texto porque incluye 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)))

x_raw <- as.numeric(datos$`Discovery year`)
y_raw <- extraer_anio(datos$`Production start year`)

cat("Registros con Discovery year:", sum(!is.na(x_raw)), "\n")
## Registros con Discovery year: 4935
cat("Registros con Production start year:", sum(!is.na(y_raw)), "\n")
## Registros con Production start year: 2233
cat("Pares completos (ambos con dato):", sum(!is.na(x_raw) & !is.na(y_raw)), "\n")
## Pares completos (ambos con dato): 1889
cat("X sin Y:", sum(!is.na(x_raw) & is.na(y_raw)), "\n")
## X sin Y: 3046
cat("Y sin X:", sum(is.na(x_raw) & !is.na(y_raw)), "\n")
## Y sin X: 344

4.Tabla de pares de valores

Se presenta un extracto (primeros 20 registros) tal como fueron extraídos del dataset, antes de cualquier depuración.

df_pares <- data.frame(x = x_raw, y = y_raw)

df_pares %>%
  head(20) %>%
  rename(`Año de Descubrimiento (X)`       = x,
         `Año de Inicio de Producción (Y)` = y) %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla de Pares de Valores**"),
    subtitle = md("Valores originales sin depurar")
  ) %>%
  tab_source_note(source_note = "Autor: Grupo 5") %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.border.top.color            = "black",
    table.border.bottom.color         = "black",
    table.border.top.style            = "solid",
    table.border.bottom.style         = "solid",
    column_labels.font.weight         = "bold",
    column_labels.border.top.color    = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color       = "black",
    heading.border.bottom.width       = px(2),
    table_body.hlines.color           = "grey",
    table_body.border.bottom.color    = "black"
  )
Tabla de Pares de Valores
Valores originales sin depurar
Año de Descubrimiento (X) Año de Inicio de Producción (Y)
1949 1951
2001 2009
1966 1969
1975 1979
1984 1987
1986 1998
1981 1981
2004 2005
NA NA
NA NA
1981 1985
1986 2002
1982 1987
2007 NA
1965 1965
2000 2010
2013 2020
2001 NA
1979 1983
1999 2021
Autor: Grupo 5

5.Gráfica

Gráfica original

plot(df_pares$x, df_pares$y,
     pch  = 20,
     col  = rgb(0.1, 0.4, 0.5, 0.4),
     xlab = "Año de Descubrimiento (X)",
     ylab = "Año de Inicio de Producción (Y)",
     main = "Relación entre Año de Descubrimiento y Año de Inicio de Producción (Datos Originales)")

Debido a la dispersión y a los múltiples valores de Y repetidos para un mismo X, se procede a aplicar una estrategia de tratamiento de datos antes de proponer un modelo.

5.1 Descarte de pares incompletos

A diferencia de un modelo polinómico, un modelo potencial trabaja en escala logarítmica (\(\ln x\), \(\ln y\)), donde cualquier valor artificial (como rellenar con la media global) distorsiona fuertemente la relación: si muchos años X distintos no tienen dato real de Y y todos se completan con el mismo valor, se genera una banda plana de puntos falsos que rompe la tendencia y hunde la correlación. Por eso, en este modelo no se rellenan los faltantes: se descartan directamente los registros que no tengan ambos valores.

df_pares <- df_pares %>% filter(!is.na(x), !is.na(y))
cat("Pares completos reales (ambos datos presentes):", nrow(df_pares), "\n")
## Pares completos reales (ambos datos presentes): 1889

5.2 Agrupación de múltiples Y por X

Cuando un mismo año de descubrimiento (X) tiene múltiples años de inicio de producción (Y), se calcula la media aritmética de todos esos Y para obtener un único par representativo (un único X, un único Y). Se conserva además el número de registros originales (n) que dieron origen a cada promedio, dato necesario para la depuración posterior.

pares <- df_pares %>%
  filter(!is.na(x), !is.na(y)) %>%
  group_by(x) %>%
  summarise(y = mean(y, na.rm = TRUE), n = n(), .groups = "drop") %>%
  arrange(x)

cat("Pares únicos (un X, un Y) para el modelo:", nrow(pares), "\n")
## Pares únicos (un X, un Y) para el modelo: 98
cat("Rango de años:", min(pares$x), "-", max(pares$x), "\n")
## Rango de años: 1905 - 2023

5.3 Estrategia de depuración

Al analizar los datos agrupados se identificaron valores atípicos (pares con años de producción muy alejados de la tendencia central) y años con un solo registro original (que no aportan representatividad, pues su valor de Y no proviene de un promedio). La estrategia adoptada consiste en:

  • Eliminar pares cuya Y supere 2 desviaciones estándar de la media, considerados valores atípicos.
  • Conservar únicamente los años con más de un registro original (n > 1), para garantizar representatividad.
media_y  <- mean(pares$y)
sd_y     <- sd(pares$y)
lim_sup  <- media_y + 2 * sd_y
lim_inf  <- media_y - 2 * sd_y

cat("Media de Y:", round(media_y, 4), "\n")
## Media de Y: 1986.281
cat("Desv. estándar de Y:", round(sd_y, 4), "\n")
## Desv. estándar de Y: 30.3537
cat("Límite superior:", round(lim_sup, 4), "\n")
## Límite superior: 2046.988
cat("Límite inferior:", round(lim_inf, 4), "\n")
## Límite inferior: 1925.574
pares_dep <- pares %>%
  filter(y >= lim_inf & y <= lim_sup) %>%
  filter(n > 1) %>%
  filter(x > 0, y > 0)   # el modelo potencial exige X, Y positivos

cat("\nPares antes de depuración:", nrow(pares), "\n")
## 
## Pares antes de depuración: 98
cat("Pares después de depuración:", nrow(pares_dep), "\n")
## Pares después de depuración: 84
cat("Pares eliminados (atípicos, un único registro, o no positivos):",
    nrow(pares) - nrow(pares_dep), "\n")
## Pares eliminados (atípicos, un único registro, o no positivos): 14

5.4 Tabla de pares depurados

pares_dep %>%
  select(x, y) %>%
  rename(`Año de Descubrimiento (X)`       = x,
         `Año de Inicio de Producción (Y)` = y) %>%
  mutate(`Año de Inicio de Producción (Y)` = round(`Año de Inicio de Producción (Y)`, 4)) %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla de Pares Depurados**"),
    subtitle = md("Año de Descubrimiento y Año de Inicio de Producción")
  ) %>%
  tab_source_note(source_note = "Autor: Grupo 5") %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.border.top.color            = "black",
    table.border.bottom.color         = "black",
    table.border.top.style            = "solid",
    table.border.bottom.style         = "solid",
    column_labels.font.weight         = "bold",
    column_labels.border.top.color    = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color       = "black",
    heading.border.bottom.width       = px(2),
    table_body.hlines.color           = "grey",
    table_body.border.bottom.color    = "black"
  )
Tabla de Pares Depurados
Año de Descubrimiento y Año de Inicio de Producción
Año de Descubrimiento (X) Año de Inicio de Producción (Y)
1927 1940.500
1928 1944.500
1932 1935.500
1937 1963.500
1938 1966.000
1941 1949.500
1944 1944.000
1945 1947.000
1946 1973.333
1948 1951.250
1949 1955.143
1951 1966.571
1952 1958.857
1953 1964.000
1954 1966.444
1955 1977.077
1956 1973.143
1957 1968.833
1958 1975.188
1959 1972.625
1960 1980.923
1961 1982.133
1962 1974.850
1963 1972.909
1964 1977.529
1965 1984.879
1966 1980.667
1967 1989.281
1968 1986.303
1969 1990.536
1970 1986.762
1971 1992.086
1972 1988.590
1973 1996.821
1974 1996.902
1975 1993.756
1976 1998.350
1977 1997.767
1978 1995.826
1979 1998.600
1980 1995.500
1981 1993.963
1982 1998.605
1983 1999.939
1984 1998.903
1985 2002.500
1986 2000.735
1987 2000.710
1988 2001.148
1989 2005.100
1990 2005.739
1991 2002.036
1992 2007.576
1993 2004.667
1994 2005.000
1995 2003.846
1996 2004.482
1997 2005.500
1998 2004.697
1999 2007.708
2000 2010.471
2001 2008.811
2002 2012.040
2003 2013.704
2004 2013.963
2005 2015.524
2006 2015.708
2007 2016.160
2008 2015.512
2009 2017.667
2010 2019.686
2011 2020.103
2012 2023.000
2013 2022.808
2014 2021.023
2015 2020.250
2016 2021.929
2017 2023.095
2018 2023.786
2019 2022.680
2020 2025.067
2021 2022.500
2022 2027.600
2023 2023.667
Autor: Grupo 5

5.5 Gráfica simplificada

plot(pares_dep$x, pares_dep$y,
     pch  = 20,
     col  = rgb(0.1, 0.4, 0.5, 0.6),
     xlab = "Año de Descubrimiento (X)",
     ylab = "Año de Inicio de Producción (Y)",
     main = "Relación entre Año de Descubrimiento y Año de Inicio de Producción (Datos Depurados)")

6.Conjetura de Modelo

Observando la gráfica de los datos depurados, se propone un Modelo de Regresión Potencial, ya que ambas variables son estrictamente positivas y la tendencia muestra un crecimiento sostenido de Y a medida que X aumenta, sin cambios de dirección. Este modelo tiene la forma:

\[y = a \cdot x^{b}\]

7.Cálculo de parámetros

El modelo potencial no es lineal en X, pero se linealiza aplicando logaritmo natural a ambos lados:

\[\ln(y) = \ln(a) + b \cdot \ln(x)\]

Esta expresión sí es lineal en \(\ln(x)\), por lo que R puede estimar \(\ln(a)\) y \(b\) mediante mínimos cuadrados con lm, y luego se recupera \(a = e^{\ln(a)}\).

m_potencial <- lm(log(y) ~ log(x), data = pares_dep)

coefs <- coef(m_potencial)
ln_a  <- coefs[1]
b     <- coefs[2]
a     <- exp(ln_a)

cat("ln(a) :", round(ln_a, 6), "\n")
## ln(a) : 0.871664
cat("a     :", round(a, 8), "\n")
## a     : 2.390887
cat("b     :", round(b, 6), "\n")
## b     : 0.886026
cat("\nEcuación del modelo:\n")
## 
## Ecuación del modelo:
cat("y =", format(a, scientific = TRUE, digits = 6), "* x^(", round(b, 6), ")\n")
## y = 2.39089e+00 * x^( 0.886026 )

8.Comparación del modelo con la realidad

Se realiza la superposición del modelo ajustado sobre los datos reales para evaluar visualmente qué tan bien representa la curva potencial el comportamiento observado.

x_grid <- seq(min(pares_dep$x), max(pares_dep$x), length.out = 400)
y_grid <- a * x_grid^b

plot(pares_dep$x, pares_dep$y,
     pch  = 20,
     col  = rgb(0.1, 0.4, 0.5, 0.6),
     xlab = "Año de Descubrimiento (X)",
     ylab = "Año de Inicio de Producción (Y)",
     main = "Superposición: Modelo Potencial y Datos Reales")

lines(x_grid, y_grid, col = "firebrick3", lwd = 3)

legend("topleft",
       legend = c("Datos reales", "Modelo potencial"),
       col    = c(rgb(0.1, 0.4, 0.5, 0.6), "firebrick3"),
       pch    = c(20, NA),
       lty    = c(NA, 1),
       lwd    = c(NA, 3),
       bty    = "n")

9.Test de Pearson

Correlación lineal (escala original)

r <- cor(pares_dep$x, pares_dep$y)
cat("Correlación de Pearson (r):", round(r, 4), "\n")
## Correlación de Pearson (r): 0.9699

Correlación en escala logarítmica

El modelo potencial se ajusta realmente sobre \(\ln(x)\) y \(\ln(y)\), por lo que también se reporta la correlación en esa escala, que es la que sustenta directamente el ajuste por mínimos cuadrados.

r_log <- cor(log(pares_dep$x), log(pares_dep$y))
cat("Correlación de Pearson en escala log-log (r):", round(r_log, 4), "\n")
## Correlación de Pearson en escala log-log (r): 0.9696

Coeficiente de determinación

Como el modelo se linealiza en escala logarítmica, el R² obtenido de m_potencial corresponde a la bondad de ajuste sobre \(\ln(x)\) y \(\ln(y)\), no sobre la escala original:

r2 <- summary(m_potencial)$r.squared
cat("R² (escala log-log):", round(r2, 4), "\n")
## R² (escala log-log): 0.9401

10.Restricciones

El dominio de la variable independiente X (año de descubrimiento) corresponde a números enteros positivos: \(X \in \mathbb{Z}^{+}\), condición que exige el modelo potencial (no está definido para \(X \le 0\)). El dominio de la variable dependiente Y (año de inicio de producción) corresponde también a años calendario positivos, y por realismo histórico se acota entre 1859 (primer pozo petrolero comercial documentado) y 2100.

Como el modelo es potencial con b > 0 (crecimiento monótono), Y aumenta de forma continua conforme X aumenta, sin cambios de dirección. Por lo tanto, basta con despejar X a partir de los límites de Y para conocer el rango válido:

\[x = \left(\frac{y}{a}\right)^{1/b}\]

y_min_dominio <- 1859
y_max_dominio <- 2100

x_min_valido <- (y_min_dominio / a)^(1 / b)
x_max_valido <- (y_max_dominio / a)^(1 / b)

cat("El modelo se mantiene dentro del dominio de Y [", y_min_dominio, ",", y_max_dominio, "]\n")
## El modelo se mantiene dentro del dominio de Y [ 1859 , 2100 ]
cat("para valores de X comprendidos entre:", round(x_min_valido, 0), "y", round(x_max_valido, 0), "\n")
## para valores de X comprendidos entre: 1830 y 2100

Conclusión de la sección: el modelo sí presenta restricciones. Es válido únicamente para años de descubrimiento comprendidos entre 1830 y 2100. Fuera de este intervalo, la curva potencial produce años de inicio de producción fuera del rango histórico/realista considerado.

11.Estimación del modelo

Aprovechando la ecuación del modelo potencial, se realizan estimaciones dentro del rango válido determinado en la sección anterior.

anio_estimar <- 2015
produccion_estimada <- a * anio_estimar^b

cat("Estimación para un yacimiento descubierto en", anio_estimar, ":\n")
## Estimación para un yacimiento descubierto en 2015 :
cat("Año estimado de inicio de producción:", round(produccion_estimada, 1), "\n")
## Año estimado de inicio de producción: 2024.1
cat("¿Dentro del dominio válido? :", anio_estimar >= x_min_valido & anio_estimar <= x_max_valido, "\n\n")
## ¿Dentro del dominio válido? : TRUE
anio_estimar2 <- 2020
produccion_estimada2 <- a * anio_estimar2^b

cat("Estimación para un yacimiento descubierto en", anio_estimar2, ":\n")
## Estimación para un yacimiento descubierto en 2020 :
cat("Año estimado de inicio de producción:", round(produccion_estimada2, 1), "\n")
## Año estimado de inicio de producción: 2028.6
cat("¿Dentro del dominio válido? :", anio_estimar2 >= x_min_valido & anio_estimar2 <= x_max_valido, "\n")
## ¿Dentro del dominio válido? : TRUE

12.Conclusión

Se presenta a continuación la tabla resumen del modelo, como base para la conclusión.

Tabla_resumen <- data.frame(
  `Variable Independiente` = "Año de Descubrimiento",
  `Variable Dependiente`   = "Año de Inicio de Producción",
  `Test Pearson`           = round(r, 4),
  `Test Pearson (log-log)` = round(r_log, 4),
  `Ecuación del modelo`    = ec,
  `Rango válido de X`      = paste0("[", round(x_min_valido, 0), ", ", round(x_max_valido, 0), "]"),
  check.names = FALSE
)

Tabla_resumen %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N°1**"),
    subtitle = md("**Resumen del modelo de regresión potencial**")
  ) %>%
  tab_source_note(source_note = md("Autor: Grupo 5")) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_options(
    table.border.top.color            = "black",
    table.border.bottom.color         = "black",
    table.border.top.style            = "solid",
    table.border.bottom.style         = "solid",
    column_labels.font.weight         = "bold",
    column_labels.border.top.color    = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color       = "black",
    heading.border.bottom.width       = px(2),
    table_body.hlines.color           = "grey",
    table_body.border.bottom.color    = "black"
  )
Tabla N°1
Resumen del modelo de regresión potencial
Variable Independiente Variable Dependiente Test Pearson Test Pearson (log-log) Ecuación del modelo Rango válido de X
Año de Descubrimiento Año de Inicio de Producción 0.9699 0.9696 y = 2.39089e+00 * x^(0.886026) [1830, 2100]
Autor: Grupo 5

Entre el año de descubrimiento (X) y el año de inicio de producción (Y) existe una relación potencial, cuya ecuación es:

y = 2.39089e+00 * x^(0.886026)

Siendo X el año de descubrimiento del yacimiento y Y el año estimado de inicio de producción. La correlación de Pearson en escala log-log es de 0.9696, lo que indica un ajuste fuerte del modelo. El modelo presenta restricciones: es válido únicamente para años de descubrimiento comprendidos entre 1830 y 2100, ya que fuera de ese intervalo el modelo produce años de inicio de producción fuera del rango histórico/realista considerado.