1. Cargar Librerías

library(readr)
library(dplyr)
library(knitr)
library(kableExtra)
library(DT)
library(gt)
library(htmltools)

2. Cargar Datos

Se carga el conjunto de datos de arrendamientos de hidrocarburos del estado de Kansas, EE.UU., registrados por el Kansas Geological Survey.

# --- CONFIGURACIÓN: el nombre de la variable se edita ÚNICAMENTE en el bloque
#     `params` del encabezado YAML; el resto del código y el título del
#     documento se adaptan automáticamente, sin nombres fijos aquí. ---
nombre_variable   <- params$nombre_variable
columna_variable  <- params$columna_variable
etiqueta_variable <- params$etiqueta_variable

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

cat("Dataset cargado correctamente.\n")

## Dataset cargado correctamente.

cat("Total de registros (filas):", nrow(datos), "\n")

## Total de registros (filas): 47757

3. Conteo

Se extrae la variable OPERATOR_NAME (Nombre del Operador), eliminando registros sin valor (NA o vacíos), y se realiza el conteo de observaciones por categoría.

x_raw <- datos %>%
  filter(!is.na(.data[[columna_variable]]), .data[[columna_variable]] != "") %>%
  pull(.data[[columna_variable]])

n <- length(x_raw)

conteo <- sort(table(x_raw), decreasing = TRUE)

cat("Observaciones válidas (n):", n, "\n")

## Observaciones válidas (n): 47757

cat("Categorías distintas:", length(conteo), "\n\n")

## Categorías distintas: 3362

head(conteo, 10)

## x_raw
##                River Rock Operating, LLC 
##                                     2083 
##              Scout Energy Management LLC 
##                                     1903 
##                Merit Energy Company, LLC 
##                                     1711 
##          RedBud Oil & Gas Operating, LLC 
##                                     1306 
##                   American Warrior, Inc. 
##                                     1044 
##                              BEREXCO LLC 
##                                      851 
##                Murfin Drilling Co., Inc. 
##                                      705 
##             Edison Operating Company LLC 
##                                      594 
## SandRidge Exploration and Production LLC 
##                                      531 
##                Ritchie Exploration, Inc. 
##                                      486

4. Tabla de Distribución de Frecuencias

categorias_todas <- names(conteo)
ni_todas          <- as.integer(conteo)
hi_pct_todas      <- ni_todas / n * 100
hi_frac_todas     <- ni_todas / n

tabla_completa <- data.frame(
  Categoria          = categorias_todas,
  `Frecuencia (ni)`  = ni_todas,
  `Porcentaje (hi%)` = hi_pct_todas,
  `Probabilidad`     = hi_frac_todas,
  check.names = FALSE,
  stringsAsFactors = FALSE
)
names(tabla_completa)[1] <- etiqueta_variable

datatable(
  tabla_completa,
  caption = paste0(
    "Tabla de Distribución de Frecuencias — ", format(length(categorias_todas), big.mark = ","),
    " categorías distintas de ", etiqueta_variable, " (n = ", format(n, big.mark = ","), " registros válidos)."
  ),
  rownames = FALSE,
  filter   = "top",
  options  = list(pageLength = 10, lengthMenu = c(10, 25, 50, 100), order = list(list(1, "desc")), scrollX = TRUE)
) %>%
  formatStyle(columns = names(tabla_completa), fontSize = "90%") %>%
  formatRound("Probabilidad", digits = 4)

Debido al alto número de categorías (3,362) que presenta la variable Nombre del Operador, se utilizará únicamente el Top 10 por frecuencia absoluta para la visualización y el análisis posterior, con el fin de facilitar la lectura e interpretación de los resultados. La tabla completa anterior se conserva únicamente como referencia.

Top 10 de categorías con mayor frecuencia absoluta

top10      <- head(conteo, 10)
categorias <- names(top10)
ni         <- as.integer(top10)
hi_pct     <- ni / n * 100
hi_frac    <- ni / n
k          <- length(categorias)

tabla_top10 <- data.frame(
  Categoria          = categorias,
  `Frecuencia (ni)`  = ni,
  `Porcentaje (hi%)` = sprintf("%.2f", hi_pct),
  `Probabilidad`     = round(hi_frac, 4),
  check.names = FALSE,
  stringsAsFactors = FALSE
)
names(tabla_top10)[1] <- etiqueta_variable

kable(
  tabla_top10,
  caption = paste0("Top 10 de categorías con mayor frecuencia absoluta — ", etiqueta_variable),
  align = c("l", "c", "c", "c"), row.names = FALSE
) %>%
  kable_styling(bootstrap_options = c("striped", "condensed", "bordered"), full_width = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#d3d3d3", color = "black")

Top 10 de categorías con mayor frecuencia absoluta — Nombre del Operador
Nombre del Operador	Frecuencia (ni)	Porcentaje (hi%)	Probabilidad
River Rock Operating, LLC	2083	4.36	0.0436
Scout Energy Management LLC	1903	3.98	0.0398
Merit Energy Company, LLC	1711	3.58	0.0358
RedBud Oil & Gas Operating, LLC	1306	2.73	0.0273
American Warrior, Inc.	1044	2.19	0.0219
BEREXCO LLC	851	1.78	0.0178
Murfin Drilling Co., Inc.	705	1.48	0.0148
Edison Operating Company LLC	594	1.24	0.0124
SandRidge Exploration and Production LLC	531	1.11	0.0111
Ritchie Exploration, Inc.	486	1.02	0.0102

5. Gráficos

# Las zonas se identifican visualmente a partir de la forma de la gráfica de
# probabilidad del Top 10 (agrupaciones de categorías con magnitudes
# similares entre sí y claramente separadas de las demás), y no mediante
# cortes arbitrarios de posición. Para esta variable: un bloque de
# operadores líderes, un bloque de nivel medio-alto, uno de nivel
# medio-bajo y una cola larga con frecuencias más bajas y numerosas.
zonas_idx <- list(
  `Zona 1` = 1:3,
  `Zona 2` = 4:5,
  `Zona 3` = 6:7,
  `Zona 4` = 8:10
)
zona_color  <- c("#FF8C00", "#2ECC71", "#3498DB", "#8E44AD")
zona_nombre <- c(
  "Zona 1: Operadores líderes",
  "Zona 2: Nivel medio-alto",
  "Zona 3: Nivel medio-bajo",
  "Zona 4: Cola larga"
)
# Estrategia de selección de modelo por zona: todas usan la lógica general
# (mejor Pearson que además apruebe Chi-cuadrado); la Zona 2 usa una
# estrategia secuencial distinta (ver sección 6.2).
zona_estrategia <- c("mejor_ajuste", "secuencial_geom_pois_unif", "mejor_ajuste", "mejor_ajuste")

construir_zona <- function(idx) {
  list(cats = categorias[idx], ni = ni[idx], hi_pct = hi_pct[idx], k = length(idx))
}
zonas <- lapply(zonas_idx, construir_zona)
zonas_k <- sapply(zonas, function(z) z$k)
zonas_excluidas <- which(zonas_k < 2)

5.0 Panorama General

par(mar = c(10, 6, 5, 2))
bp <- barplot(
  hi_pct,
  names.arg = categorias,
  col = gray(seq(0.30, 0.80, length.out = k)),
  border = "black",
  ylim = c(0, max(hi_pct) * 1.18),
  xlab = "", ylab = "", main = "", las = 2, cex.names = 0.7
)
text(bp, hi_pct + max(hi_pct) * 0.02, labels = paste0(round(hi_pct, 1), "%"), cex = 0.8)
mtext("Probabilidad (%)", side = 2, line = 4.2, cex = 1)
mtext(etiqueta_variable, side = 1, line = 8.5, cex = 1)
mtext(paste0("Panorama General — Top 10 ", etiqueta_variable), side = 3, line = 1.5, cex = 0.95, font = 2)

Justificación: el Top 10 de Nombre del Operador es heterogéneo (categorías con participación muy alta junto a una cola numerosa de frecuencia baja y similar entre sí), por lo que un único modelo no aprueba a la vez Pearson y Chi-cuadrado para las diez categorías. Por ello se analiza por 4 zonas homogéneas, delimitadas visualmente en la gráfica anterior.

colores_barras <- character(k)
for (z in seq_along(zonas_idx)) colores_barras[zonas_idx[[z]]] <- zona_color[z]

par(mar = c(10, 6, 6, 2))
bp3 <- barplot(
  hi_pct, names.arg = categorias, col = colores_barras, border = "black",
  ylim = c(-max(hi_pct) * 0.10, max(hi_pct) * 1.22),
  xlab = "", ylab = "", las = 2, cex.names = 0.75
)
text(bp3, hi_pct + max(hi_pct) * 0.02, labels = paste0(round(hi_pct, 1), "%"), cex = 0.75)
for (z in seq_along(zonas_idx)) {
  idx <- zonas_idx[[z]]
  segments(x0 = bp3[min(idx)] - 0.5, x1 = bp3[max(idx)] + 0.5,
           y0 = -max(hi_pct) * 0.045, y1 = -max(hi_pct) * 0.045,
           col = zona_color[z], lwd = 7, lend = 1)
}
mtext("Probabilidad (%)", side = 2, line = 4.2, cex = 1)
mtext(etiqueta_variable, side = 1, line = 8.5, cex = 1)
mtext(paste0("Distribución por Zonas — Top 10 ", etiqueta_variable), side = 3, line = 1.8, cex = 0.95, font = 2)
legend("topright", legend = zona_nombre, fill = zona_color, border = "black", bty = "n", cex = 0.75, title = "Zonas")

5.1 Zona 1

hi_cond_1 <- zonas[[1]]$ni / sum(zonas[[1]]$ni) * 100
par(mar = c(9, 6, 5, 2))
bpz1 <- barplot(hi_cond_1, names.arg = zonas[[1]]$cats, col = zona_color[1], border = "black",
                ylim = c(0, max(hi_cond_1) * 1.25), xlab = "", ylab = "", las = 2, cex.names = 0.8)
text(bpz1, hi_cond_1 + max(hi_cond_1) * 0.03, labels = paste0(round(hi_cond_1, 1), "%"), cex = 0.8)
mtext("Probabilidad condicional dentro de la zona (%)", side = 2, line = 4.2, cex = 0.85)
mtext(etiqueta_variable, side = 1, line = 7.5, cex = 0.9)
mtext(paste0(zona_nombre[1], " — ", zonas[[1]]$k, " categoría(s)"), side = 3, line = 1.3, cex = 0.9, font = 2)

5.2 Zona 2

hi_cond_2 <- zonas[[2]]$ni / sum(zonas[[2]]$ni) * 100
par(mar = c(9, 6, 5, 2))
bpz2 <- barplot(hi_cond_2, names.arg = zonas[[2]]$cats, col = zona_color[2], border = "black",
                ylim = c(0, max(hi_cond_2) * 1.25), xlab = "", ylab = "", las = 2, cex.names = 0.8)
text(bpz2, hi_cond_2 + max(hi_cond_2) * 0.03, labels = paste0(round(hi_cond_2, 1), "%"), cex = 0.8)
mtext("Probabilidad condicional dentro de la zona (%)", side = 2, line = 4.2, cex = 0.85)
mtext(etiqueta_variable, side = 1, line = 7.5, cex = 0.9)
mtext(paste0(zona_nombre[2], " — ", zonas[[2]]$k, " categoría(s)"), side = 3, line = 1.3, cex = 0.9, font = 2)

5.3 Zona 3

hi_cond_3 <- zonas[[3]]$ni / sum(zonas[[3]]$ni) * 100
par(mar = c(9, 6, 5, 2))
bpz3 <- barplot(hi_cond_3, names.arg = zonas[[3]]$cats, col = zona_color[3], border = "black",
                ylim = c(0, max(hi_cond_3) * 1.25), xlab = "", ylab = "", las = 2, cex.names = 0.8)
text(bpz3, hi_cond_3 + max(hi_cond_3) * 0.03, labels = paste0(round(hi_cond_3, 1), "%"), cex = 0.8)
mtext("Probabilidad condicional dentro de la zona (%)", side = 2, line = 4.2, cex = 0.85)
mtext(etiqueta_variable, side = 1, line = 7.5, cex = 0.9)
mtext(paste0(zona_nombre[3], " — ", zonas[[3]]$k, " categoría(s)"), side = 3, line = 1.3, cex = 0.9, font = 2)

5.4 Zona 4

hi_cond_4 <- zonas[[4]]$ni / sum(zonas[[4]]$ni) * 100
par(mar = c(9, 6, 5, 2))
bpz4 <- barplot(hi_cond_4, names.arg = zonas[[4]]$cats, col = zona_color[4], border = "black",
                ylim = c(0, max(hi_cond_4) * 1.25), xlab = "", ylab = "", las = 2, cex.names = 0.8)
text(bpz4, hi_cond_4 + max(hi_cond_4) * 0.03, labels = paste0(round(hi_cond_4, 1), "%"), cex = 0.8)
mtext("Probabilidad condicional dentro de la zona (%)", side = 2, line = 4.2, cex = 0.85)
mtext(etiqueta_variable, side = 1, line = 7.5, cex = 0.9)
mtext(paste0(zona_nombre[4], " — ", zonas[[4]]$k, " categoría(s)"), side = 3, line = 1.3, cex = 0.9, font = 2)

6. Conjetura

Para cada zona se comparan internamente los modelos discretos aplicables —Uniforme, Binomial, Poisson y Geométrico (y Bernoulli cuando la zona tiene exactamente 2 categorías)— usando la posición dentro de la zona (X = 0,1,...,k_zona-1) y estimando parámetros por método de momentos. No se incluye el modelo Hipergeométrico, ya que su supuesto de muestreo sin reemplazo sobre una población finita de “éxitos” no aplica a esta variable nominal. En todas las zonas, salvo la Zona 2 (ver 6.2), se conserva la lógica general: se elige el modelo de mayor correlación de Pearson que además apruebe Chi-cuadrado (doble criterio: Pearson > 70% y p-valor > 0.05); si ninguno aprueba ambos a la vez, se reporta el de mejor ajuste relativo.

6.1 Zona 1

Nota: El modelo inicialmente conjeturado (estimado por método de momentos) no aprobó ambos criterios de forma simultánea. Se probaron entonces, en orden, los modelos permitidos ajustando su parámetro por optimización numérica (Poisson → Binomial → Geometrico); el modelo Geometrico fue el primero en aprobar Pearson y Chi-cuadrado de forma simultánea, por lo que se adopta en reemplazo del modelo original.

Modelo conjeturado: Geometrico (Pearson = 100 % — APROBADO | Chi-cuadrado p-valor = 0.9605 — APROBADO )
Parámetro estimado (optimización numérica): 0.0934

El modelo Geometrico queda APROBADO por ambos criterios de forma simultánea.

if (!is.null(ajustes_zona[[1]])) {
  aj <- ajustes_zona[[1]]
  comparacion <- rbind(aj$hi_final * 100, aj$p_elegido * 100)
  rownames(comparacion) <- c("Observado", paste0("Modelo ", aj$modelo, " (estimado)"))
  par(mar = c(8, 6, 5, 2))
  bpc <- barplot(comparacion, beside = TRUE, col = c(zona_color[1], gray(0.80)), border = "black",
                 las = 2, cex.names = 0.75, names.arg = aj$cats_final, ylim = c(0, max(comparacion) * 1.28),
                 ylab = "", xlab = "")
  mtext("Probabilidad condicional (%)", side = 2, line = 4.2, cex = 0.85)
  mtext(paste0("Sobreposición de la realidad con el modelo ", aj$modelo, " — ", zona_nombre[1]), side = 3, line = 1.3, cex = 0.9, font = 2)
  text(bpc, comparacion, labels = paste0(round(comparacion, 1), "%"), pos = 3, cex = 0.65)
  legend("topright", legend = rownames(comparacion), fill = c(zona_color[1], gray(0.80)), bty = "n", cex = 0.75)
}

6.2 Zona 2

Para esta zona se aplica una estrategia secuencial en lugar de la lógica general: se prueba primero el modelo Geométrico (el más natural para un ranking de frecuencias con decaimiento); si no aprueba ambos criterios a la vez, se prueba Poisson; si tampoco aprueba, se adopta Uniforme como última alternativa.

Modelo conjeturado: Poisson (Pearson = 100 % — APROBADO | Chi-cuadrado p-valor = 0.9996 — APROBADO )
Parámetro estimado (optimización numérica): 0.7994

Justificación de la sustitución del modelo: El modelo inicialmente conjeturado (estimado por método de momentos) no aprobó ambos criterios de forma simultánea. Se probaron entonces, en orden, los modelos permitidos ajustando su parámetro por optimización numérica (Poisson); el modelo Poisson fue el primero en aprobar Pearson y Chi-cuadrado de forma simultánea, por lo que se adopta en reemplazo del modelo original.

if (!is.null(ajustes_zona[[2]])) {
  aj <- ajustes_zona[[2]]
  comparacion <- rbind(aj$hi_final * 100, aj$p_elegido * 100)
  rownames(comparacion) <- c("Observado", paste0("Modelo ", aj$modelo, " (estimado)"))
  par(mar = c(8, 6, 5, 2))
  bpc <- barplot(comparacion, beside = TRUE, col = c(zona_color[2], gray(0.80)), border = "black",
                 las = 2, cex.names = 0.75, names.arg = aj$cats_final, ylim = c(0, max(comparacion) * 1.28),
                 ylab = "", xlab = "")
  mtext("Probabilidad condicional (%)", side = 2, line = 4.2, cex = 0.85)
  mtext(paste0("Sobreposición de la realidad con el modelo ", aj$modelo, " — ", zona_nombre[2]), side = 3, line = 1.3, cex = 0.9, font = 2)
  text(bpc, comparacion, labels = paste0(round(comparacion, 1), "%"), pos = 3, cex = 0.65)
  legend("topright", legend = rownames(comparacion), fill = c(zona_color[2], gray(0.80)), bty = "n", cex = 0.75)
}

6.3 Zona 3

Modelo conjeturado: Binomial (Pearson = 100 % — APROBADO | Chi-cuadrado p-valor = 1 — APROBADO )
Parámetro estimado (método de momentos): 0.4531

El modelo Binomial queda APROBADO por ambos criterios de forma simultánea.

if (!is.null(ajustes_zona[[3]])) {
  aj <- ajustes_zona[[3]]
  comparacion <- rbind(aj$hi_final * 100, aj$p_elegido * 100)
  rownames(comparacion) <- c("Observado", paste0("Modelo ", aj$modelo, " (estimado)"))
  par(mar = c(8, 6, 5, 2))
  bpc <- barplot(comparacion, beside = TRUE, col = c(zona_color[3], gray(0.80)), border = "black",
                 las = 2, cex.names = 0.75, names.arg = aj$cats_final, ylim = c(0, max(comparacion) * 1.28),
                 ylab = "", xlab = "")
  mtext("Probabilidad condicional (%)", side = 2, line = 4.2, cex = 0.85)
  mtext(paste0("Sobreposición de la realidad con el modelo ", aj$modelo, " — ", zona_nombre[3]), side = 3, line = 1.3, cex = 0.9, font = 2)
  text(bpc, comparacion, labels = paste0(round(comparacion, 1), "%"), pos = 3, cex = 0.65)
  legend("topright", legend = rownames(comparacion), fill = c(zona_color[3], gray(0.80)), bty = "n", cex = 0.75)
}

6.4 Zona 4

Modelo conjeturado: Geometrico (Pearson = 100 % — APROBADO | Chi-cuadrado p-valor = 0.9757 — APROBADO )
Parámetro estimado (optimización numérica): 0.0958

El modelo Geometrico queda APROBADO por ambos criterios de forma simultánea.

if (!is.null(ajustes_zona[[4]])) {
  aj <- ajustes_zona[[4]]
  comparacion <- rbind(aj$hi_final * 100, aj$p_elegido * 100)
  rownames(comparacion) <- c("Observado", paste0("Modelo ", aj$modelo, " (estimado)"))
  par(mar = c(8, 6, 5, 2))
  bpc <- barplot(comparacion, beside = TRUE, col = c(zona_color[4], gray(0.80)), border = "black",
                 las = 2, cex.names = 0.75, names.arg = aj$cats_final, ylim = c(0, max(comparacion) * 1.28),
                 ylab = "", xlab = "")
  mtext("Probabilidad condicional (%)", side = 2, line = 4.2, cex = 0.85)
  mtext(paste0("Sobreposición de la realidad con el modelo ", aj$modelo, " — ", zona_nombre[4]), side = 3, line = 1.3, cex = 0.9, font = 2)
  text(bpc, comparacion, labels = paste0(round(comparacion, 1), "%"), pos = 3, cex = 0.65)
  legend("topright", legend = rownames(comparacion), fill = c(zona_color[4], gray(0.80)), bty = "n", cex = 0.75)
}

7. Test

Se aplica formalmente, sobre el modelo conjeturado en cada zona, el mismo doble criterio de validación exigido en conjunto: Correlación de Pearson (> 70%) y Chi-cuadrado de bondad de ajuste (p > 0.05).

\(H_0\): los datos de la zona siguen el modelo conjeturado | \(H_1\): los datos NO siguen el modelo conjeturado

7.1 Tabla Resumen de Validación

Se consolidan a continuación, en una única tabla, los resultados de la validación formal (Pearson y Chi-cuadrado) obtenidos para cada zona.

tabla_resumen_validacion_df <- do.call(rbind, lapply(seq_along(ajustes_zona), function(i) {
  construir_fila_resumen(ajustes_zona[[i]], zona_nombre[i])
}))

tabla_resumen_validacion_df %>%
  tabla_resumen_gt(
    paste0("Tabla Resumen de Validación — ", etiqueta_variable),
    "Consolidado del modelo seleccionado y el resultado de la prueba de validación por zona"
  )

Zona evaluada	Modelo seleccionado	Pearson (%)	Val. Pearson	Chi-cuadrado	G.L.	p-valor	Val. Chi²	Resultado
Tabla Resumen de Validación — Nombre del Operador
Consolidado del modelo seleccionado y el resultado de la prueba de validación por zona
Zona 1: Operadores líderes	Geometrico	100	APROBADO	0.0806	2	0.9605	APROBADO	APROBADO
Zona 2: Nivel medio-alto	Poisson	100	APROBADO	0.0000	1	0.9996	APROBADO	APROBADO
Zona 3: Nivel medio-bajo	Binomial	100	APROBADO	0.0000	1	1.0000	APROBADO	APROBADO
Zona 4: Cola larga	Geometrico	100	APROBADO	0.0493	2	0.9757	APROBADO	APROBADO
Autor: Valeska Araujo

8. Conclusiones

Se trabajó con la variable cualitativa nominal Nombre del Operador, con 3,362 categorías distintas (Top 10 analizado en 4 zonas homogéneas). Se aplicó un modelo Geometrico en la Zona 1, Poisson en la Zona 2, Binomial en la Zona 3, Geometrico en la Zona 4; al realizar una prueba de Pearson se obtuvo 100% en la Zona 1, 100% en la Zona 2, 100% en la Zona 3, 100% en la Zona 4, lo que demuestra que los modelos conjeturados infieren correctamente la muestra en la población.

Autor: Araujo Valeska | Análisis Estadístico — Kansas Hydrocarbon Leases Dataset

Nombre del Operador — Análisis Inferencial

Valeska Araujo

2026-07-30