1 Librerías

library(readr)
library(dplyr)
library(ggplot2)
library(gt)
library(stringr)
library(DT)
library(knitr)
library(kableExtra)

cat("Librerías cargadas: readr, dplyr, ggplot2, gt, stringr, DT, knitr, kableExtra")
Librerías cargadas: readr, dplyr, ggplot2, gt, stringr, DT, knitr, kableExtra

2 Carga de Datos

ruta_archivo <- file.choose()
datos <- read_csv(ruta_archivo, show_col_types = FALSE)

# Rellenamiento de celdas faltantes con la mediana de cada variable
if (any(is.na(datos$TOWNSHIP))) {
  datos$TOWNSHIP[is.na(datos$TOWNSHIP)] <-
    median(datos$TOWNSHIP, na.rm = TRUE)
}
if (any(is.na(datos$LATITUDE))) {
  datos$LATITUDE[is.na(datos$LATITUDE)] <-
    median(datos$LATITUDE, na.rm = TRUE)
}

# Se descartan registros con valores no positivos
datos <- datos %>%
  filter(TOWNSHIP > 0, LATITUDE > 0)

str(datos[, c("TOWNSHIP", "LATITUDE")])
tibble [47,757 × 2] (S3: tbl_df/tbl/data.frame)
 $ TOWNSHIP: num [1:47757] 33 15 29 26 33 17 30 26 23 16 ...
 $ LATITUDE: num [1:47757] 37.1 38.8 37.5 37.8 37.1 ...

3 Selección de Variables

Variable Dependiente (Y): Latitud Geográfica. Se seleccionó como variable dependiente porque representa la ubicación geográfica absoluta de cada arrendamiento medida en grados decimales. Es el resultado que se desea analizar y explicar a partir de otra variable.

Variable Independiente (X): Township. Se eligió como variable independiente porque representa la posición norte-sur de cada arrendamiento dentro del sistema de agrimensura rectangular (PLSS) de Kansas. En general, a mayor número de Township, mayor es el desplazamiento hacia el sur del arrendamiento — cada unidad equivale aproximadamente a 6 millas — por lo que se espera que esta variable influya directamente sobre la latitud geográfica.

4 Tablas

4.1 Tabla de valores

La siguiente tabla contiene todos los 47757 registros del conjunto, mostrando los valores de Township y Latitud de cada arrendamiento.

datatable(
  datos %>% select(TOWNSHIP, LATITUDE) %>% rename("Township (X)" = TOWNSHIP, "Latitud (Y)" = LATITUDE),
  caption = htmltools::tags$caption(
    style = "caption-side: top; text-align: left; font-size: 16px;
             font-weight: 700; color:#1F2A33;",
    "Tabla N\u00b01: Todos los valores \u2014 Township y Latitud Geogr\u00e1fica"
  ),
  rownames = FALSE,
  class    = "display compact stripe hover",
  options  = list(
    pageLength = 10,
    dom        = "ltip",
    columnDefs = list(list(className = "dt-center", targets = "_all"))
  )
) %>%
  formatRound(columns = c("Latitud (Y)"), digits = 4)

La nube de puntos visualizada a continuación no permite conjeturar un modelo con claridad, por esta razón se optó por realizar una estrategia para depurar los datos y obtener un gráfico más limpio.

par(mar = c(5, 5, 4, 2))
plot(datos$TOWNSHIP, datos$LATITUDE,
     pch = 16, cex = 0.35,
     col = adjustcolor("#2E86AB", alpha.f = 0.15),
     xlab = "X (Township)",
     ylab = "Y (Latitud, grados decimales)",
     main = "Nube de puntos (todos los valores, sin agrupar)",
     cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()

Estrategia:

  1. Generación de pares únicos
  2. Llenar espacios donde no haya valores
  3. Donde se encuentren valores diversos usar indicadores de posición (mediana)
  4. Presentar nueva tabla con la mediana de los valores

4.2 Tabla de pares (X, Y) usando la mediana

tabla_xy <- datos %>%
  mutate(percentil = ntile(TOWNSHIP, 25)) %>%
  group_by(percentil) %>%
  summarise(
    n_pozos   = n(),
    x_mediana = median(TOWNSHIP),
    y_mediana = median(LATITUDE),
    .groups   = "drop"
  ) %>%
  arrange(percentil)

cat("Total de pares (x\u0303, y\u0303) obtenidos, uno por percentil:", nrow(tabla_xy), "\n")
Total de pares (x̃, ỹ) obtenidos, uno por percentil: 25 
datatable(
  tabla_xy %>% select(x_mediana, y_mediana) %>%
    rename("Township (X)" = x_mediana, "Mediana de Latitud (Y)" = y_mediana),
  caption = htmltools::tags$caption(
    style = "caption-side: top; text-align: left; font-size: 16px;
             font-weight: 700; color:#1F2A33;",
    "Tabla N\u00b02: Pares (x\u0303, y\u0303) \u2014 Mediana de Latitud Geogr\u00e1fica por Percentil de Township"
  ),
  rownames = FALSE,
  class    = "display compact stripe hover",
  options  = list(
    pageLength = 10,
    dom        = "ltip",
    columnDefs = list(list(className = "dt-center", targets = "_all"))
  )
) %>%
  formatRound(columns = c("Mediana de Latitud (Y)"), digits = 4)

4.3 Máximos, mínimos y resumen general

max_x <- max(tabla_xy$x_mediana)
min_x <- min(tabla_xy$x_mediana)
max_y <- max(tabla_xy$y_mediana)
min_y <- min(tabla_xy$y_mediana)

data.frame(
  Variable = c("TOWNSHIP (X)", "LATITUDE (Y)"),
  Minimo   = round(c(min_x, min_y), 4),
  Maximo   = round(c(max_x, max_y), 4),
  Rango    = round(c(max_x - min_x, max_y - min_y), 4),
  Mediana  = round(c(median(tabla_xy$x_mediana),
                     median(tabla_xy$y_mediana)), 4)
) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N\u00b02: Resumen General de las Variables (pares x\u0303, y\u0303)**")
  ) %>%
  tab_source_note(source_note = "Autor: Leslye Quinchiguango") %>%
  cols_align(align = "center", everything())
Tabla N°2: Resumen General de las Variables (pares x̃, ỹ)
Variable Minimo Maximo Rango Mediana
TOWNSHIP (X) 4.0000 35.000 31.0000 25.0000
LATITUDE (Y) 37.0236 39.735 2.7114 37.8802
Autor: Leslye Quinchiguango

5 Gráfica de Dispersión

Con los datos resumidos por Township mediante la mediana, la tendencia logarítmica decreciente se aprecia con mayor claridad.

par(mar = c(5, 5, 4, 2))
plot(tabla_xy$x_mediana, tabla_xy$y_mediana,
     pch = 19, col = "#2E86AB",
     xlab = "X\u0303 (Township)",
     ylab = "Y\u0303 (Mediana Latitud, grados decimales)",
     main = "Producci\u00f3n acumulada en funci\u00f3n de los a\u00f1os activos",
     cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()

6 Conjetura

La nube de puntos muestra una tendencia decreciente: al aumentar el número de Township la Latitud disminuye progresivamente. Aunque la curvatura no es evidente a simple vista — debido a que el rango de Township en Kansas (1–35) es estrecho — la tasa de descenso es mayor en los primeros townships y se aplana hacia los valores más altos, comportamiento inherente al sistema PLSS por la curvatura de la Tierra. Por lo anterior, se conjetura el modelo:

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

con pendiente \(b < 0\) (relación inversamente proporcional a tasa decreciente).

7 Cálculo de Parámetros

El modelo se ajusta sobre los pares (x̃, ỹ) por percentil (Tabla N°1). El uso de la mediana por percentil en lugar de los 47757 registros crudos reduce la influencia de coordenadas con errores de georreferenciación y permite capturar la tendencia central de la relación geográfica.

x <- tabla_xy$x_mediana
y <- tabla_xy$y_mediana

modelo <- lm(y ~ log(x))
a <- coef(modelo)[1]
b <- coef(modelo)[2]

cat("Modelo ajustado: y =", round(a, 4), "+", round(b, 4), "* ln(x)\n")
Modelo ajustado: y = 42.3318 + -1.4171 * ln(x)
summary(modelo)

Call:
lm(formula = y ~ log(x))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.63229 -0.14343 -0.00274  0.19019  0.27773 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 42.33180    0.27590  153.43 < 0.0000000000000002 ***
log(x)      -1.41708    0.08932  -15.87   0.0000000000000702 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2268 on 23 degrees of freedom
Multiple R-squared:  0.9163,    Adjusted R-squared:  0.9126 
F-statistic: 251.7 on 1 and 23 DF,  p-value: 0.00000000000007021

8 Sobreponer Modelo con la Realidad

x_seq <- seq(min(x), max(x), length.out = 500)
pred  <- predict(modelo,
                 newdata  = data.frame(x = x_seq),
                 interval = "confidence",
                 level    = 0.95)

par(mar = c(5, 5, 4, 2))
plot(x, y,
     pch = 19, col = "#2E86AB",
     xlab = "X\u0303 (Mediana Township)",
     ylab = "Y\u0303 (Mediana Latitud, grados decimales)",
     main = "Gr\u00e1fica N\u00b03: Sobreponer Modelo con la Realidad \u2014 Pares (x\u0303, y\u0303)",
     cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")

# Banda de confianza al 95%
polygon(c(x_seq, rev(x_seq)),
        c(pred[, "lwr"], rev(pred[, "upr"])),
        col    = rgb(0.5, 0.5, 0.5, 0.2),
        border = NA)

# Curva ajustada
lines(x_seq, pred[, "fit"], col = "#E74C3C", lwd = 3)

legend("topright",
       legend = c("Datos (x\u0303, y\u0303)",
                  "Modelo Logar\u00edtmico",
                  "I.C. 95%"),
       col    = c("#2E86AB", "#E74C3C", "gray"),
       pch    = c(16, NA, 15),
       lwd    = c(NA, 3, NA),
       pt.cex = c(1, NA, 2),
       bty    = "n")
box()

La curva logarítmica sigue de cerca la dirección de los puntos, ajustándose a los pares (x̃, ỹ) y confirmando la relación logarítmica decreciente entre el número de Township y la Latitud Geográfica.

9 Test de Bondad (Pearson)

El coeficiente de correlación de Pearson, calculado sobre el modelo linealizado \(\ln(x)\), debe superar 0.7 en valor absoluto para aceptar el modelo logarítmico. Se evalúa sobre los pares (x̃, ỹ) por percentil.

r_medianas <- cor.test(log(tabla_xy$x_mediana), tabla_xy$y_mediana)

data.frame(
  Conjunto   = c("Pares (x\u0303, y\u0303) por percentil"),
  r          = round(c(r_medianas$estimate), 4),
  Supera_0.7 = c(abs(r_medianas$estimate) > 0.7)
) %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla N\u00b03: Test de Bondad de Ajuste**"),
    subtitle = "Umbral |r| > 0.7"
  ) %>%
  tab_source_note(source_note = "Autor: Leslye Quinchiguango") %>%
  cols_align(align = "center", everything())
Tabla N°3: Test de Bondad de Ajuste
Umbral |r| > 0.7
Conjunto r Supera_0.7
Pares (x̃, ỹ) por percentil -0.9572 TRUE
Autor: Leslye Quinchiguango

Sobre los pares (x̃, ỹ), \(r =\) -0.957, por encima del umbral de 0.7, por lo que el modelo logarítmico se acepta para describir la tendencia central de la relación entre el Township y la Latitud Geográfica.

10 Restricciones

El modelo \(y = a + b \cdot \ln(x)\) con \(b < 0\) impone restricciones matemáticas y geográficas que se sintetizan en la siguiente tabla:

x_max_valido <- exp((-90 - a) / b)   # Township máximo para y = -90°

tabla_restricciones <- data.frame(
  Variable       = c("X (TOWNSHIP)", "Y (LATITUDE)"),
  Dominio        = c("x ∈ ℤ⁺",
                     "-90° ≤ y ≤ 90°"),
  Rango_valido   = c("1 \u2264 x \u2264 3.60e+40",
                     "-90° ≤ y ≤ 90°"),
  Rango_invalido = c("0 \u2264 x < 3.60e+40",
                     "-90° < y < 90°")
)

kable(tabla_restricciones,
      col.names = c("Variable", "Dominio", "Rango válido", "Rango inválido"),
      align     = "lccc", escape = TRUE) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE, position = "center", font_size = 14) %>%
  column_spec(1, bold = TRUE, width = "5cm") %>%
  column_spec(2, width = "5cm") %>%
  column_spec(3, width = "5cm") %>%
  column_spec(4, width = "5cm") %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2")
Variable Dominio Rango válido Rango inválido
X (TOWNSHIP) x ∈ ℤ⁺ 1 ≤ x ≤ 3.60e+40 0 ≤ x < 3.60e+40
Y (LATITUDE) -90° ≤ y ≤ 90° -90° ≤ y ≤ 90° -90° < y < 90°

La condición \(x > 0\) es matemáticamente necesaria para que \(\ln(x)\) esté definido; como Township toma valores enteros positivos (\(x \geq 1\)), esta condición siempre se cumple. El rango válido de X se extiende hasta \(3.59 \times 10^{40}\), valor a partir del cual la latitud predicha caería por debajo de \(-90°\), superando el límite geográfico terrestre. La latitud Y se mantiene dentro de su dominio \([-90°, 90°]\) para cualquier Township dentro del rango válido.

11 Estimación

x_est <- round(median(tabla_xy$x_mediana), 1)
y_est <- a + b * log(x_est)

if (b >= 0) {
  ecuacion <- paste0("y = ", round(a, 4), " + ", round(b, 4), " \u00b7 ln(x)")
} else {
  ecuacion <- paste0("y = ", round(a, 4), " - ", abs(round(b, 4)), " \u00b7 ln(x)")
}

data.frame(
  Ecuacion_del_Modelo      = ecuacion,
  X_estimado_township      = x_est,
  Y_estimado_latitud_grados = round(y_est, 4)
) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N\u00b04: Estimaci\u00f3n Puntual de Latitud Geogr\u00e1fica**")
  ) %>%
  tab_source_note(source_note = "Autor: Leslye Quinchiguango") %>%
  cols_align(align = "center", everything())
Tabla N°4: Estimación Puntual de Latitud Geográfica
Ecuacion_del_Modelo X_estimado_township Y_estimado_latitud_grados
y = 42.3318 - 1.4171 · ln(x) 25 37.7704
Autor: Leslye Quinchiguango

Para un arrendamiento ubicado en el Township 25 (mediana del conjunto analizado), el modelo estima una latitud geográfica de aproximadamente 37.7704° decimales.

12 Conclusiones

Entre el Township (X) y la Latitud Geográfica (Y) existe una relación de tipo logarítmica decreciente cuya ecuación matemática es:

\[y = 42.3318 - 1.4171 · ln(x)\]

Siendo X el número de Township dentro del sistema de agrimensura rectangular PLSS de Kansas, e Y la latitud geográfica del arrendamiento medida en grados decimales. Se evidencia una relación logarítmica negativa de alta intensidad entre las variables, consistente con el desplazamiento progresivo hacia el sur que representa cada unidad adicional de Township dentro del sistema PLSS.

Autor: Leslye Quinchiguango — Análisis Estadístico, Kansas Hydrocarbon Leases Dataset