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

Se definió el Township (TOWNSHIP, cuantitativa discreta, adimensional) como variable independiente / causa (x), pues representa la posición norte-sur de cada arrendamiento dentro del sistema de agrimensura rectangular (PLSS) de Kansas, donde cada unidad equivale aproximadamente a 6 millas de desplazamiento hacia el sur. La Latitud Geográfica (LATITUDE, cuantitativa continua, medida en grados decimales) actúa como variable dependiente / efecto (y), ya que refleja la ubicación geográfica absoluta del arrendamiento. Al aumentar x, y disminuye de forma logarítmica decreciente, lo que justifica el modelo propuesto.

4 Tabla de Pares de Valores

4.1 Mediana de la variable Y por percentil de X

Los 47757 registros se agrupan en 25 percentiles según Township y se calcula la mediana de ambas variables por grupo, obteniendo un único par (x̃, ỹ) por percentil. Se emplea la mediana por ser más robusta ante valores atípicos, dado que TOWNSHIP es discreta con pocos valores únicos y LATITUDE puede presentar irregularidades de georreferenciación.

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,
  caption = htmltools::tags$caption(
    style = "caption-side: top; text-align: left; font-size: 16px;
             font-weight: 700; color:#1F2A33;",
    "Tabla N\u00b01: 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("x_mediana", "y_mediana"), digits = 4)

4.2 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

5.1 Todos los valores (47757 registros)

Se grafican todos los registros del conjunto, sin agrupar ni resumir, para visualizar la nube completa de puntos.

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 = "Gr\u00e1fica N\u00b01: Nube de Puntos \u2014 Todos los Valores",
     cex.main = 0.9, frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
box()

5.2 Pares (x̃, ỹ) por percentil

Con los datos resumidos por percentil 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 (Mediana Township)",
     ylab = "Y\u0303 (Mediana Latitud, grados decimales)",
     main = "Gr\u00e1fica N\u00b02: Nube de Puntos \u2014 Pares (x\u0303, y\u0303)",
     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 con curvatura logarítmica: al aumentar el número de Township la Latitud disminuye progresivamente, con una tasa de descenso mayor en los primeros townships que se aplana hacia los valores más altos. Este comportamiento es inherente al sistema PLSS, donde cada township adicional representa un desplazamiento hacia el sur no perfectamente lineal debido a 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),
  R2         = round(c(r_medianas$estimate^2), 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 R2 Supera_0.7
Pares (x̃, ỹ) por percentil -0.9572 0.9163 TRUE
Autor: Leslye Quinchiguango

Sobre los pares (x̃, ỹ), \(r =\) -0.957 (\(R^2 =\) 0.916), 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. Con un coeficiente de correlación de Pearson \(r =\) -0.957, el modelo explica 91.6% de la variabilidad de Y (\(R^2 =\) 0.916); el resto corresponde a factores no capturados como imprecisiones de georreferenciación o variaciones internas dentro de cada township, evidenciando una relación logarítmica negativa de alta intensidad entre las variables.

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