Modelo de Regresión Polinómica


0.- Carga de librerías

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

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

2.- Definición de las 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 determina en qué zona geográfica se realizó la exploración. La variable Longitude (longitud geográfica) actúa como variable dependiente o efecto (Y), ya que la ubicación longitudinal de los yacimientos está influenciada por los ciclos históricos y geopolíticos de exploración petrolera a nivel mundial.

x_raw <- as.numeric(datos$`Discovery year`)
y_raw <- as.numeric(datos$Longitude)

cat("Registros con Discovery year:", sum(!is.na(x_raw)), "\n")
## Registros con Discovery year: 4935
cat("Registros con Longitude:", sum(!is.na(y_raw)), "\n")
## Registros con Longitude: 7537
cat("Pares completos (ambos con dato):", sum(!is.na(x_raw) & !is.na(y_raw)), "\n")
## Pares completos (ambos con dato): 4641
cat("X sin Y:", sum(!is.na(x_raw) & is.na(y_raw)), "\n")
## X sin Y: 294
cat("Y sin X:", sum(is.na(x_raw) & !is.na(y_raw)), "\n")
## Y sin X: 2896

3.- Tabla 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,
         `Longitud Geográfica (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) Longitud Geográfica (Y)
1949 16.71667
2001 -39.61200
1966 -36.88900
1975 -36.26200
1984 -39.96100
1986 -39.74600
1981 -36.73400
2004 -36.07300
NA NA
NA -43.47230
1981 -40.52400
1986 -36.73300
1982 -36.56500
2007 -38.10900
1965 -38.17100
2000 -39.85300
2013 -42.46900
2001 -41.89300
1979 -38.97100
1999 -58.17800
Autor: Grupo 5

4.- Gráfica de dispersión

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 = "Longitud Geográfica (Y)",
     main = "Relación entre el Año de Descubrimiento y la Longitud Geográfica (Datos Originales)")

5.- Tratamiento de datos

Debido a la complejidad y a la fuerte dispersión de los puntos observados en la gráfica anterior, se procede a aplicar una estrategia de tratamiento de datos antes de proponer un modelo.

Relleno de valores faltantes. Los registros que tienen X pero no tienen Y se rellenan con la media aritmética global de Y, que es la mejor estimación cuando no se dispone del dato real.

media_y_global <- mean(df_pares$y, na.rm = TRUE)
cat("Media global de Longitude (Y):", round(media_y_global, 4), "\n")
## Media global de Longitude (Y): -54.6526
df_pares$y[is.na(df_pares$y) & !is.na(df_pares$x)] <- media_y_global
cat("Pares disponibles tras relleno:", sum(!is.na(df_pares$x) & !is.na(df_pares$y)), "\n")
## Pares disponibles tras relleno: 4935

Agrupación de múltiples Y por X. Cuando un mismo año (X) tiene múltiples longitudes (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: 125
cat("Rango de años:", min(pares$x), "-", max(pares$x), "\n")
## Rango de años: 1869 - 2023

Estrategia de depuración. Al analizar los datos agrupados se identificaron valores atípicos (pares con longitudes muy alejadas 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 longitud (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.
# Calcular límites para valores atípicos (±2 desviaciones estándar)
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: -53.2393
cat("Desv. estándar de Y:", round(sd_y, 4), "\n")
## Desv. estándar de Y: 38.7936
cat("Límite superior:", round(lim_sup, 4), "\n")
## Límite superior: 24.3479
cat("Límite inferior:", round(lim_inf, 4), "\n")
## Límite inferior: -130.8266
# Filtrar valores atípicos y conservar solo años con más de un registro original
pares_dep <- pares %>%
  filter(y >= lim_inf & y <= lim_sup) %>%
  filter(n > 1)

cat("\nPares antes de depuración:", nrow(pares), "\n")
## 
## Pares antes de depuración: 125
cat("Pares después de depuración:", nrow(pares_dep), "\n")
## Pares después de depuración: 116
cat("Pares eliminados (atípicos o con un único registro):", nrow(pares) - nrow(pares_dep), "\n")
## Pares eliminados (atípicos o con un único registro): 9

Transformación de X desde cero. Dado que los datos van desde el año 1869 hasta 2023, la curvatura del modelo polinómico no es visible en la escala original ya que trabajamos con valores muy grandes de X. Para evidenciar correctamente la curvatura, se transforma X restando el año mínimo, de modo que X empiece desde 0. Para las predicciones finales se vuelve a la variable original sumando el año mínimo.

x_min <- min(pares_dep$x)
pares_dep$x_orig <- pares_dep$x          # guardar X original
pares_dep$x      <- pares_dep$x - x_min  # X transformada desde 0

cat("X mínimo original:", x_min, "\n")
## X mínimo original: 1869
cat("Rango de X transformada:", min(pares_dep$x), "-", max(pares_dep$x), "\n")
## Rango de X transformada: 0 - 154

5.1.- Tabla pares de valores simplificada

pares_dep %>%
  select(x, y) %>%
  rename(`Año de Descubrimiento (X)` = x,
         `Longitud Geográfica (Y)`   = y) %>%
  mutate(`Longitud Geográfica (Y)` = round(`Longitud Geográfica (Y)`, 4)) %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla de Pares Depurados**"),
    subtitle = md("Año de Descubrimiento y Longitud Geográfica")
  ) %>%
  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 Longitud Geográfica
Año de Descubrimiento (X) Longitud Geográfica (Y)
0 -81.1545
18 -120.3586
20 -104.4832
30 -46.7469
32 -104.9395
35 -112.0510
36 -80.8102
40 -65.1795
41 -118.6501
42 -119.7239
43 -78.9868
45 -94.7827
46 -97.5110
47 -96.3902
48 -107.2199
49 -94.5553
51 -99.7704
52 -100.3496
53 -97.4209
54 -110.0782
56 -93.6454
58 -53.1959
59 -74.9300
60 -79.7388
61 -97.5647
62 -27.5353
63 -76.1102
66 -88.3248
67 -90.5222
68 -86.7212
69 -57.6879
70 -68.6410
71 -69.0887
72 -63.9294
73 -87.6575
74 -94.1861
75 -73.2791
76 -74.9884
77 -67.6844
78 -106.3878
79 -56.2014
80 -95.5411
81 -90.7631
82 -99.0192
83 -92.3924
84 -87.0251
85 -84.0192
86 -83.3799
87 -81.2610
88 -83.7155
89 -71.7962
90 -66.5756
91 -77.9652
92 -38.0788
93 -68.1975
94 -58.9798
95 -47.6465
96 -38.0495
97 -48.4943
98 -24.0001
99 -42.6327
100 -49.6160
101 -37.4139
102 -21.3941
103 -16.1785
104 -50.8687
105 -23.7429
106 -16.8435
107 -76.0603
108 -59.7735
109 -31.7067
110 -41.3217
111 -47.9180
112 -39.8525
113 -21.0103
114 -29.3466
115 -39.5066
116 -44.5174
117 -22.0576
118 -18.4606
119 -27.9513
120 7.1355
121 5.9354
122 -22.0110
123 -12.4867
124 -28.1336
125 -40.9822
126 -25.9760
127 -34.0524
128 -2.2996
129 -15.8990
130 -17.3036
131 21.5589
132 -8.5892
133 -8.1713
134 -17.0716
135 -24.3896
136 -15.9556
137 -14.9397
138 -3.5778
139 -51.6520
140 -74.8394
141 2.2130
142 11.0168
143 -7.5526
144 -12.8655
145 7.1889
146 -15.9569
147 -38.5953
148 -41.4956
149 -5.3741
150 -2.7491
151 -20.1708
152 -11.9346
153 -6.6244
154 -4.7329
Autor: Grupo 5

5.2.- Gráfica de dispersión simplificada

plot(pares_dep$x, pares_dep$y,
     pch  = 20,
     col  = rgb(0.1, 0.4, 0.5, 0.6),
     xlab = paste0("Años desde ", x_min, " (X transformada)"),
     ylab = "Longitud Geográfica (Y)",
     main = "Relación entre el Año de Descubrimiento y la Longitud Geográfica")

Con la transformación y depuración, la curvatura del modelo polinómico se evidencia con mayor claridad.

6.- Conjetura

Observando la gráfica de los datos depurados, se propone un Modelo de Regresión Polinómica de grado 3, ya que los datos presentan una tendencia curvilínea con cambios de dirección. Este modelo tiene la forma:

\[y = a + b_1x + b_2x^2 + b_3x^3\]

7.- Parámetros

Como se trata de una curva, no se puede aplicar lm directamente sobre X. Sin embargo, tratando cada potencia de X (\(x\), \(x^2\), \(x^3\)) como una variable independiente adicional, la ecuación se vuelve lineal en esos términos, por lo que R puede calcular automáticamente los coeficientes mediante mínimos cuadrados con lm.

m_poli3 <- lm(y ~ poly(x, 3, raw = TRUE), data = pares_dep)

Pendiente e intercepto

b <- coef(m_poli3)

cat("Intercepto (a)  :", round(b[1], 6), "\n")
## Intercepto (a)  : -74.44101
cat("Pendiente b1    :", round(b[2], 8), "\n")
## Pendiente b1    : -1.614955
cat("Pendiente b2    :", round(b[3], 10), "\n")
## Pendiente b2    : 0.02927269
cat("Pendiente b3    :", round(b[4], 12), "\n")
## Pendiente b3    : -0.0001046021
cat("\nEcuación del modelo:\n")
## 
## Ecuación del modelo:
cat("y =", round(b[1], 2),
    "+ (", round(b[2], 6), ")x",
    "+ (", round(b[3], 8), ")x²",
    "+ (", round(b[4], 10), ")x³\n")
## y = -74.44 + ( -1.614955 )x + ( 0.02927269 )x² + ( -0.0001046021 )x³

8.- Comparación de la realidad con el modelo

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

x_grid <- seq(min(pares_dep$x), max(pares_dep$x), length.out = 400)
y_grid <- predict(m_poli3, newdata = data.frame(x = x_grid))

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 = "Longitud Geográfica (Y)",
     main = "Superposición: Modelo Polinómico y Datos Reales")

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

legend("topright",
       legend = c("Datos reales", "Modelo polinómico"),
       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 Bondad

Correlación lineal

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

Coeficiente de determinación

En modelos no lineales como el polinómico, el coeficiente de determinación R² no aplica directamente como medida de bondad de ajuste, ya que este indicador fue diseñado para modelos de regresión lineal, donde la razón de cambio (pendiente) entre X y Y es constante. En una curva, esa razón de cambio varía en cada punto, por lo que R² pierde su interpretación. En su lugar, la calidad del ajuste se evalúa visualmente mediante la superposición del modelo con los datos reales, y mediante la correlación de Pearson entre las variables.

10.- Restricciones

1. Dominio de X

El dominio de la variable independiente X (año de descubrimiento) corresponde al conjunto de los números enteros: \(X \in \mathbb{Z}\).

2. Dominio de Y

El dominio de la variable dependiente Y (longitud geográfica) corresponde al intervalo: \(Y \in [-180°, 180°]\).

3. Pregunta

¿Existe algún valor del dominio de X que, al sustituirse en el modelo matemático, genere un valor de Y fuera de su dominio?

Respuesta: Sí

4. Cálculo de la restricción utilizando la ecuación del modelo

Dado que el modelo es una curva polinómica de grado 3, existe la posibilidad de que, para valores de X fuera del rango de los datos observados, el valor estimado de Y se salga del dominio \([-180°, 180°]\). Para comprobarlo, se resuelve el polinomio igualándolo a los límites del dominio de Y, y se despeja X:

\[a + b_1x + b_2x^2 + b_3x^3 = 180 \qquad \text{y} \qquad a + b_1x + b_2x^2 + b_3x^3 = -180\]

# Raíces reales donde el modelo predice Y = 180°
raices_sup <- polyroot(c(b[1] - 180, b[2], b[3], b[4]))
raices_sup_reales <- Re(raices_sup)[abs(Im(raices_sup)) < 1e-6]

# Raíces reales donde el modelo predice Y = -180°
raices_inf <- polyroot(c(b[1] + 180, b[2], b[3], b[4]))
raices_inf_reales <- Re(raices_inf)[abs(Im(raices_inf)) < 1e-6]

cat("Valores de X donde el modelo predice Y = 180°:\n")
## Valores de X donde el modelo predice Y = 180°:
print(round(sort(raices_sup_reales), 0))
## [1] -65
cat("\nValores de X donde el modelo predice Y = -180°:\n")
## 
## Valores de X donde el modelo predice Y = -180°:
print(round(sort(raices_inf_reales), 0))
## [1] 232

A partir de las raíces anteriores se determina el intervalo de X para el cual la curva se mantiene dentro del dominio válido de Y:

x_amplio <- seq(min(pares_dep$x) - 100, max(pares_dep$x) + 100, by = 1)
y_amplio <- predict(m_poli3, newdata = data.frame(x = x_amplio))

x_validos <- x_amplio[y_amplio >= -180 & y_amplio <= 180]

x_min_valido <- min(x_validos)
x_max_valido <- max(x_validos)

cat("El modelo se mantiene dentro del dominio de Y [-180°, 180°]\n")
## El modelo se mantiene dentro del dominio de Y [-180°, 180°]
cat("para valores de X comprendidos entre:", x_min_valido, "y", x_max_valido, "\n")
## para valores de X comprendidos entre: -64 y 232

Conclusión de la sección: el modelo sí presenta restricciones. Es válido únicamente para años comprendidos entre -64 y 232. Fuera de este intervalo, el polinomio produce valores de longitud que exceden el dominio geográfico real \([-180°, 180°]\), por lo que el modelo no debe usarse para estimar fuera de ese rango.

11.- Estimación

Aprovechando la ecuación del modelo polinómico, se realizan estimaciones dentro del rango válido determinado en la sección anterior (entre -64 y 232).

# Nota: X fue transformada restando x_min, por lo que hay que convertir
# el año original a X transformada antes de predecir

# Estimación para el año 2025
anio_estimar      <- 2025
x_estimar_transf  <- anio_estimar - x_min   # volver a X transformada
longitud_estimada <- predict(m_poli3, newdata = data.frame(x = x_estimar_transf))

cat("Estimación para el año", anio_estimar, "(X transformada =", x_estimar_transf, "):\n")
## Estimación para el año 2025 (X transformada = 156 ):
cat("Longitud geográfica estimada:", round(longitud_estimada, 4), "°\n")
## Longitud geográfica estimada: -11.107 °
cat("¿Dentro del dominio válido? :", anio_estimar >= x_min_valido & anio_estimar <= x_max_valido, "\n\n")
## ¿Dentro del dominio válido? : FALSE
# Estimación para el año 2030
anio_estimar2      <- 2030
x_estimar_transf2  <- anio_estimar2 - x_min
longitud_estimada2 <- predict(m_poli3, newdata = data.frame(x = x_estimar_transf2))

cat("Estimación para el año", anio_estimar2, "(X transformada =", x_estimar_transf2, "):\n")
## Estimación para el año 2030 (X transformada = 161 ):
cat("Longitud geográfica estimada:", round(longitud_estimada2, 4), "°\n")
## Longitud geográfica estimada: -12.2055 °
cat("¿Dentro del dominio válido? :", anio_estimar2 >= x_min_valido & anio_estimar2 <= x_max_valido, "\n")
## ¿Dentro del dominio válido? : FALSE

12.- Conclusión

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

Ecuacion <- paste0(
  "y = ", round(b[1], 2),
  " + (", round(b[2], 6), ")x",
  " + (", round(b[3], 8), ")x²",
  " + (", round(b[4], 10), ")x³"
)

Tabla_resumen <- data.frame(
  `Variable Independiente` = "Año de Descubrimiento",
  `Variable Dependiente`   = "Longitud Geográfica",
  `Test Pearson`           = round(r, 2),
  `Ecuación del modelo`    = Ecuacion,
  `Rango válido de X`      = paste0("[", x_min_valido, ", ", x_max_valido, "]"),
  check.names = FALSE
)

Tabla_resumen %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N°1**"),
    subtitle = md("**Resumen del modelo de regresión polinómica**")
  ) %>%
  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 polinómica
Variable Independiente Variable Dependiente Test Pearson Ecuación del modelo Rango válido de X
Año de Descubrimiento Longitud Geográfica 0.83 y = -74.44 + (-1.614955)x + (0.02927269)x² + (-0.0001046021)x³ [-64, 232]
Autor: Grupo 5

Entre la longitud geográfica (Y) y el año de descubrimiento (X) existe una relación de tipo no lineal (polinómica de tercer grado), cuyo modelo matemático es:

y = -74.44 + (-1.614955)x + (0.02927269)x² + (-0.0001046021)x³

Siendo Y la longitud geográfica donde se ubica el yacimiento y X el año de descubrimiento, donde SÍ existen restricciones: el modelo es válido únicamente para valores de X comprendidos entre -64 y 232, ya que fuera de ese intervalo el modelo produce valores de longitud fuera del dominio real \([-180°, 180°]\).