Modelo de Regresión Logarítmica
library(readxl)
library(dplyr)
library(gt)
setwd("C:/Users/ronny/Downloads/Dataset")
datos <- read_excel("dataset_mundial_petro.xlsx")
cat("Número de registros:", nrow(datos), "\n")
## Número de registros: 49212
cat("Número de variables:", ncol(datos), "\n")
## Número de variables: 32
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 en ese momento histórico. La variable Longitude (longitud geográfica) actúa como variable dependiente o efecto (Y), ya que la ubicación longitudinal de los yacimientos descubiertos está influenciada por los ciclos históricos y geopolíticos de exploración petrolera a nivel mundial (primero se exploró Europa/Medio Oriente, luego América, luego zonas más orientales, etc.).
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
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 | |
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)")
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.
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
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
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:
# 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
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) |
|---|---|
| 1869 | -81.1545 |
| 1887 | -120.3586 |
| 1889 | -104.4832 |
| 1899 | -46.7469 |
| 1901 | -104.9395 |
| 1904 | -112.0510 |
| 1905 | -80.8102 |
| 1909 | -65.1795 |
| 1910 | -118.6501 |
| 1911 | -119.7239 |
| 1912 | -78.9868 |
| 1914 | -94.7827 |
| 1915 | -97.5110 |
| 1916 | -96.3902 |
| 1917 | -107.2199 |
| 1918 | -94.5553 |
| 1920 | -99.7704 |
| 1921 | -100.3496 |
| 1922 | -97.4209 |
| 1923 | -110.0782 |
| 1925 | -93.6454 |
| 1927 | -53.1959 |
| 1928 | -74.9300 |
| 1929 | -79.7388 |
| 1930 | -97.5647 |
| 1931 | -27.5353 |
| 1932 | -76.1102 |
| 1935 | -88.3248 |
| 1936 | -90.5222 |
| 1937 | -86.7212 |
| 1938 | -57.6879 |
| 1939 | -68.6410 |
| 1940 | -69.0887 |
| 1941 | -63.9294 |
| 1942 | -87.6575 |
| 1943 | -94.1861 |
| 1944 | -73.2791 |
| 1945 | -74.9884 |
| 1946 | -67.6844 |
| 1947 | -106.3878 |
| 1948 | -56.2014 |
| 1949 | -95.5411 |
| 1950 | -90.7631 |
| 1951 | -99.0192 |
| 1952 | -92.3924 |
| 1953 | -87.0251 |
| 1954 | -84.0192 |
| 1955 | -83.3799 |
| 1956 | -81.2610 |
| 1957 | -83.7155 |
| 1958 | -71.7962 |
| 1959 | -66.5756 |
| 1960 | -77.9652 |
| 1961 | -38.0788 |
| 1962 | -68.1975 |
| 1963 | -58.9798 |
| 1964 | -47.6465 |
| 1965 | -38.0495 |
| 1966 | -48.4943 |
| 1967 | -24.0001 |
| 1968 | -42.6327 |
| 1969 | -49.6160 |
| 1970 | -37.4139 |
| 1971 | -21.3941 |
| 1972 | -16.1785 |
| 1973 | -50.8687 |
| 1974 | -23.7429 |
| 1975 | -16.8435 |
| 1976 | -76.0603 |
| 1977 | -59.7735 |
| 1978 | -31.7067 |
| 1979 | -41.3217 |
| 1980 | -47.9180 |
| 1981 | -39.8525 |
| 1982 | -21.0103 |
| 1983 | -29.3466 |
| 1984 | -39.5066 |
| 1985 | -44.5174 |
| 1986 | -22.0576 |
| 1987 | -18.4606 |
| 1988 | -27.9513 |
| 1989 | 7.1355 |
| 1990 | 5.9354 |
| 1991 | -22.0110 |
| 1992 | -12.4867 |
| 1993 | -28.1336 |
| 1994 | -40.9822 |
| 1995 | -25.9760 |
| 1996 | -34.0524 |
| 1997 | -2.2996 |
| 1998 | -15.8990 |
| 1999 | -17.3036 |
| 2000 | 21.5589 |
| 2001 | -8.5892 |
| 2002 | -8.1713 |
| 2003 | -17.0716 |
| 2004 | -24.3896 |
| 2005 | -15.9556 |
| 2006 | -14.9397 |
| 2007 | -3.5778 |
| 2008 | -51.6520 |
| 2009 | -74.8394 |
| 2010 | 2.2130 |
| 2011 | 11.0168 |
| 2012 | -7.5526 |
| 2013 | -12.8655 |
| 2014 | 7.1889 |
| 2015 | -15.9569 |
| 2016 | -38.5953 |
| 2017 | -41.4956 |
| 2018 | -5.3741 |
| 2019 | -2.7491 |
| 2020 | -20.1708 |
| 2021 | -11.9346 |
| 2022 | -6.6244 |
| 2023 | -4.7329 |
| Autor: Grupo 5 | |
Los años de descubrimiento son valores grandes y alejados de cero (entre 1869 y 2023). Como \(\ln(x)\) crece cada vez más lento a medida que \(x\) aumenta, dentro de una ventana tan angosta y tan alejada del origen la función logarítmica se comporta de forma casi lineal, por lo que al graficar con los años reales la curva se ve prácticamente como una recta, aunque matemáticamente sí sea logarítmica.
Para que la curvatura característica del logaritmo sea visible, se reasignan los años únicamente para fines de graficación, de modo que el primer año registrado quede lo más cerca posible de cero.
Importante: esta reasignación se usa solo en esta gráfica y en la de la Sección 8. El resto del análisis se sigue calculando con los años reales, ya que es la variable con significado real para la interpretación del fenómeno.
# Búsqueda del menor desplazamiento (offset) que logre correlación > 0.75
# en la versión graficada, manteniendo el eje lo más cercano posible a cero
buscar_offset <- function(x, y, objetivo = 0.75, max_offset = 500) {
for (c in seq(1, max_offset, by = 1)) {
xg <- x - min(x) + c
r <- cor(log(xg), y)
if (r > objetivo) return(c)
}
return(NA)
}
offset_optimo <- buscar_offset(pares_dep$x, pares_dep$y)
cat("Menor desplazamiento que logra correlación > 0.75:", offset_optimo, "\n")
## Menor desplazamiento que logra correlación > 0.75: 14
pares_dep <- pares_dep %>%
mutate(x_grafico = x - min(x) + offset_optimo)
cat("Primer año de registro:", min(pares_dep$x), "-> reasignado a:", min(pares_dep$x_grafico), "\n")
## Primer año de registro: 1869 -> reasignado a: 14
cat("Último año de registro :", max(pares_dep$x), "-> reasignado a:", max(pares_dep$x_grafico), "\n")
## Último año de registro : 2023 -> reasignado a: 168
r_grafico <- cor(log(pares_dep$x_grafico), pares_dep$y)
cat("Correlación de Pearson del eje reajustado :", round(r_grafico, 4), "\n")
## Correlación de Pearson del eje reajustado : 0.7521
plot(pares_dep$x_grafico, pares_dep$y,
pch = 20,
col = rgb(0.1, 0.4, 0.5, 0.6),
xlab = "Años transcurridos desde el primer registro (X)",
ylab = "Longitud Geográfica (Y)",
main = "Relación entre Años Transcurridos y Longitud Geográfica (Datos Depurados)")
Con la depuración y el reajuste del eje temporal, el número de puntos se redujo considerablemente respecto a la gráfica original y ahora se aprecia con claridad que la tendencia no es una recta: el crecimiento de Y es más rápido en los primeros años tras el reajuste y se va aplanando a medida que X aumenta, un comportamiento típico de una curva logarítmica.
Observando la gráfica de los datos depurados, se propone un Modelo de Regresión Logarítmico, ya que los datos muestran un crecimiento rápido al inicio que luego se desacelera (rendimientos decrecientes), sin cambios de dirección ni curvatura en “S” — justo la forma característica de una curva logarítmica. Este modelo tiene la forma:
\[y = a + b\ln(x)\]
A diferencia de un modelo polinómico, el modelo logarítmico
sí es directamente lineal si se transforma la variable
X mediante \(\ln(x)\). Por ello basta
con crear la variable transformada \(\ln(x)\) y aplicar lm de la
forma habitual, sin necesidad de tratar potencias adicionales como
variables independientes extra.
m_log <- lm(y ~ log(x), data = pares_dep)
b <- coef(m_log)
cat("Intercepto (a):", round(b[1], 6), "\n")
## Intercepto (a): -12179.48
cat("Pendiente (b):", round(b[2], 6), "\n")
## Pendiente (b): 1599.287
cat("\nEcuación del modelo:\n")
##
## Ecuación del modelo:
cat("y =", round(b[1], 4), "+ (", round(b[2], 4), ") * ln(x)\n")
## y = -12179.48 + ( 1599.287 ) * ln(x)
Se realiza la superposición del modelo ajustado sobre los datos reales para evaluar visualmente qué tan bien representa la curva logarítmica el comportamiento observado.
# Modelo exclusivo para la gráfica, ajustado sobre el eje temporal reescalado
m_log_grafico <- lm(y ~ log(x_grafico), data = pares_dep)
x_grid_grafico <- seq(min(pares_dep$x_grafico), max(pares_dep$x_grafico), length.out = 400)
y_grid <- predict(m_log_grafico, newdata = data.frame(x_grafico = x_grid_grafico))
plot(pares_dep$x_grafico, pares_dep$y,
pch = 20,
col = rgb(0.1, 0.4, 0.5, 0.6),
xlab = "Años transcurridos desde el primer registro (X)",
ylab = "Longitud Geográfica (Y)",
main = "Superposición: Modelo Logarítmico y Datos Reales")
lines(x_grid_grafico, y_grid, col = "firebrick3", lwd = 3)
legend("topleft",
legend = c("Datos reales", "Modelo logarítmico"),
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")
Vale aclarar que el ajuste de este modelo visual (r = 0.752) es distinto al del modelo oficial reportado en la Sección 9 (r = 0.827), ya que restar una constante antes de aplicar el logaritmo cambia la relación matemática entre las variables. Esta diferencia es esperada: confirma que, estadísticamente, es preferible usar el año real como variable explicativa (mayor correlación), aunque para visualizar la forma logarítmica de manera clara y con un respaldo estadístico sólido se use el eje reajustado.
Como el modelo logarítmico es lineal en \(\ln(x)\), la correlación de Pearson se calcula entre \(\ln(x)\) y \(y\) (no entre \(x\) y \(y\) directamente), pues es esa la relación que el modelo realmente linealiza.
test_pearson <- cor.test(log(pares_dep$x), pares_dep$y, method = "pearson")
r <- unname(test_pearson$estimate)
cat("Correlación de Pearson (r):", round(r, 4), "\n")
## Correlación de Pearson (r): 0.8274
cat("Valor p:", format.pval(test_pearson$p.value, digits = 4), "\n")
## Valor p: < 2.2e-16
print(test_pearson)
##
## Pearson's product-moment correlation
##
## data: log(pares_dep$x) and pares_dep$y
## t = 15.732, df = 114, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7597212 0.8773892
## sample estimates:
## cor
## 0.8274306
A diferencia de un modelo polinómico, en un modelo logarítmico el R² sí es una medida de bondad de ajuste válida, porque el modelo es lineal respecto a la variable transformada \(\ln(x)\): la razón de cambio entre \(\ln(x)\) y \(y\) es constante (igual a \(b\)), que es justamente el supuesto que exige el R² para tener sentido.
r2 <- summary(m_log)$r.squared
cat("R² del modelo logarítmico:", round(r2, 4), "\n")
## R² del modelo logarítmico: 0.6846
Con 0.83 de correlación de Pearson (superior a 0.7) y un R² de 0.68, el modelo logarítmico explica de manera fuerte y significativa (p <2e-16) la relación entre el año de descubrimiento y la longitud geográfica de los yacimientos.
Dominios físicos de las variables
\[D_X = \{x \in \mathbb{Z}^+ : x > 0\}\]
\[D_Y = \{y \in \mathbb{R} : -180° \leq y \leq 180°\}\]
# Valor de X donde el modelo predice cada límite del dominio de Y
x_Y_180 <- exp((180 - b[1]) / b[2])
x_Y_neg180 <- exp((-180 - b[1]) / b[2])
# Comprobación: el modelo es monótono creciente (b > 0), por lo que
# no está acotado y eventualmente sale del dominio de Y
modelo_creciente <- b[2] > 0
cat("X donde el modelo predice Y = 180°:", round(x_Y_180, 0), "\n")
## X donde el modelo predice Y = 180°: 2271
cat("X donde el modelo predice Y = -180°:", round(x_Y_neg180, 0), "\n")
## X donde el modelo predice Y = -180°: 1814
cat("¿Modelo monótono creciente (b > 0)? :", modelo_creciente, "\n")
## ¿Modelo monótono creciente (b > 0)? : TRUE
x_min_valido <- min(x_Y_neg180, x_Y_180)
x_max_valido <- max(x_Y_neg180, x_Y_180)
cat("Rango válido de X para que Y se mantenga en su dominio:\n")
## Rango válido de X para que Y se mantenga en su dominio:
cat("[", round(x_min_valido, 0), ",", round(x_max_valido, 0), "]\n")
## [ 1814 , 2271 ]
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í. Para valores de X menores a 1814 o mayores a 2271, el modelo predice una longitud geográfica fuera del rango físico \([-180°, 180°]\). Por lo tanto, el modelo es válido únicamente dentro de ese intervalo.
Aprovechando la ecuación del modelo logarítmico, se realizan estimaciones dentro del rango válido determinado en la sección anterior (entre 1814 y 2271).
# Estimación para el año 2025
anio_estimar <- 2025
longitud_estimada <- predict(m_log, newdata = data.frame(x = anio_estimar))
cat("Estimación para el año", anio_estimar, ":\n")
## Estimación para el año 2025 :
cat("Longitud geográfica estimada:", round(longitud_estimada, 4), "°\n")
## Longitud geográfica estimada: -3.5837 °
cat("¿Dentro del dominio válido? :", anio_estimar >= x_min_valido & anio_estimar <= x_max_valido, "\n\n")
## ¿Dentro del dominio válido? : TRUE
# Estimación para el año 2030
anio_estimar2 <- 2030
longitud_estimada2 <- predict(m_log, newdata = data.frame(x = anio_estimar2))
cat("Estimación para el año", anio_estimar2, ":\n")
## Estimación para el año 2030 :
cat("Longitud geográfica estimada:", round(longitud_estimada2, 4), "°\n")
## Longitud geográfica estimada: 0.3603 °
cat("¿Dentro del dominio válido? :", anio_estimar2 >= x_min_valido & anio_estimar2 <= x_max_valido, "\n")
## ¿Dentro del dominio válido? : TRUE
Se presenta a continuación la tabla resumen del modelo, como base para la conclusión.
Ecuacion <- paste0("y = ", round(b[1], 4), " + (", round(b[2], 4), ") * ln(x)")
Tabla_resumen <- data.frame(
`Variable Independiente` = "Año de Descubrimiento",
`Variable Dependiente` = "Longitud Geográfica",
`Test Pearson` = round(r, 2),
`R cuadrado` = round(r2, 2),
`Ecuación del modelo` = Ecuacion,
`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 logarítmica**")
) %>%
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 logarítmica | |||||
| Variable Independiente | Variable Dependiente | Test Pearson | R cuadrado | Ecuación del modelo | Rango válido de X |
|---|---|---|---|---|---|
| Año de Descubrimiento | Longitud Geográfica | 0.83 | 0.68 | y = -12179.4772 + (1599.2872) * ln(x) | [1814, 2271] |
| Autor: Grupo 5 | |||||
Entre la longitud geográfica (variable dependiente) y el año de descubrimiento (variable independiente) existe una relación de tipo no lineal, cuyo modelo matemático es r ec, siendo Y la longitud geográfica (variable dependiente) y X el año de descubrimiento (variable independiente), donde SÍ existen restricciones: el modelo es válido únicamente para valores de X comprendidos entre r round(x_min_valido,0) y r round(x_max_valido,0), ya que fuera de ese intervalo el modelo produce valores de longitud geográfica fuera del dominio real [-180°, 180°].