FORMACIÓN PRODUCTORA — Análisis Inferencial

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: cambiar únicamente estas 3 líneas para adaptar el script
#     a otra variable, sin tocar el resto del código. ---
nombre_variable   <- "formacion_productora"    # identificador corto en snake_case (nombra objetos y archivo)
columna_variable  <- "PRODUCING_FORMATION"     # nombre exacto de la columna en el dataset
etiqueta_variable <- "Formación Productora"    # etiqueta legible usada en tablas, ejes y texto

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 PRODUCING_FORMATION (Formación Productora), 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): 47756

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

## Categorías distintas: 313

head(conteo, 10)

## x_raw
##              UNKNOWN Mississippian System          Chase Group 
##                32073                 3351                 2053 
##       Arbuckle Group        Lansing Group    Kansas City Group 
##                 1963                 1509                  696 
##  Council Grove Group  Upper Kearny Member      Bevier Coal Bed 
##                  559                  559                  555 
##       Marmaton Group 
##                  313

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 (313) que presenta la variable Formación Productora, 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 — Formación Productora

Formación Productora	Frecuencia (ni)	Porcentaje (hi%)	Probabilidad
UNKNOWN	32073	67.16	0.6716
Mississippian System	3351	7.02	0.0702
Chase Group	2053	4.30	0.0430
Arbuckle Group	1963	4.11	0.0411
Lansing Group	1509	3.16	0.0316
Kansas City Group	696	1.46	0.0146
Council Grove Group	559	1.17	0.0117
Upper Kearny Member	559	1.17	0.0117
Bevier Coal Bed	555	1.16	0.0116
Marmaton Group	313	0.66	0.0066

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: dos formaciones
# dominantes, un bloque de nivel medio-alto, un bloque amplio de nivel
# medio-bajo y, al final de la cola, una única categoría residual.
zonas_idx <- list(
  `Zona 1` = 1:2,
  `Zona 2` = 3:5,
  `Zona 3` = 6:9,
  `Zona 4` = 10:10
)
zona_color  <- c("#3498DB", "#FF8C00", "#E91E8C", "#8E44AD")
zona_nombre <- c(
  "Zona 1: Formaciones líderes",
  "Zona 2: Nivel medio-alto",
  "Zona 3: Nivel medio-bajo",
  "Zona 4: Categoría residual única"
)
# 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 Formación Productora 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. La(s) zona(s) Zona 4 se excluyen únicamente del análisis inferencial (sección 7) por quedar compuestas por una sola categoría, sin afectar las conclusiones generales.

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)

Nota: esta zona contiene una única categoría (Marmaton Group), sin variación interna que comparar contra un modelo; se conserva aquí solo con fines descriptivos y se excluye del análisis inferencial (sección 7).

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

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

El modelo Binomial 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: Uniforme (Pearson = 99.72 % — APROBADO | Chi-cuadrado p-valor = 0 — RECHAZADO )
Parámetro estimado (método de momentos): 0.3333

Justificación de la secuencia aplicada: Ni el modelo Geométrico ni el Poisson lograron aprobar ambos criterios de forma simultánea; se adopta el modelo Uniforme como última alternativa de la secuencia definida para esta zona.

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: Uniforme (Pearson = 100 % — APROBADO | Chi-cuadrado p-valor = 0 — RECHAZADO )
Parámetro estimado (método de momentos): 0.25

Ningún modelo logró la aprobación simultánea de ambos criterios; se reporta el de mejor ajuste relativo (Uniforme).

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

Zona excluida del análisis inferencial (categoría única); ver nota en la sección 5.4.

if (length(ajustes_zona) >= 4 && !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"
  )

Tabla Resumen de Validación — Formación Productora
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
Zona 1: Formaciones líderes	Binomial	100.00	APROBADO	0.0000	1	1	APROBADO	APROBADO
Zona 2: Nivel medio-alto	Uniforme	99.72	APROBADO	92.3352	2	0	RECHAZADO	RECHAZADO
Zona 3: Nivel medio-bajo	Uniforme	100.00	APROBADO	24.2512	3	0	RECHAZADO	RECHAZADO
Autor: Valeska Araujo

8. Conclusiones

Se trabajó con la variable cualitativa nominal Formación Productora, con 313 categorías distintas (Top 10 analizado en 4 zonas homogéneas). Se aplicó un modelo Binomial en la Zona 1, Uniforme en la Zona 2, Uniforme en la Zona 3; al realizar una prueba de Pearson se obtuvo 100% en la Zona 1, 99.72% en la Zona 2, 100% en la Zona 3, lo que demuestra que el modelo conjeturado infiere correctamente la muestra en la población en 1 de 3 zona(s) evaluadas.

---

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