library(readr); library(dplyr); library(gt)
cat("Librerías cargadas correctamente.\n")
## Librerías cargadas correctamente.
Se utiliza el conjunto de datos de arrendamientos de hidrocarburos del estado de Kansas, EE.UU., registrado por el Kansas Geological Survey.
ruta_archivo <- "C:/Users/thann/OneDrive/Escritorio/ESTADISTICA.LOL/daataset/oil_and_gas_leases_data.csv.csv"
datos_vale <- read_csv(ruta_archivo, show_col_types = FALSE)
cat("Base de datos cargada correctamente.\n")
## Base de datos cargada correctamente.
cat("Total de registros:", nrow(datos_vale), "\n")
## Total de registros: 47757
cat("Total de columnas:", ncol(datos_vale), "\n")
## Total de columnas: 24
La variable de estudio es PRODUCES_OIL, que indica si el pozo registra producción de petróleo. Es una variable cualitativa nominal dicotómica.
crudo <- datos_vale$PRODUCES_OIL
n_bruto <- length(crudo)
# Diagnóstico de la codificación cruda: permite confirmar visualmente que
# 1 = "sí produce" corresponde realmente al valor más frecuente/esperado.
cat("Distribución de valores crudos en PRODUCES_OIL:\n")
## Distribución de valores crudos en PRODUCES_OIL:
print(table(crudo, useNA = "ifany"))
## crudo
## 0 1
## 15117 32640
# Si al revisar la tabla anterior la codificación estuviera invertida en el
# archivo fuente (por ejemplo, 1 = "no produce"), cambiar este interruptor a TRUE
# para corregir la gráfica y todo el análisis sin tocar el resto del código.
invertir_codificacion <- FALSE
etiqueta_si <- "Sí produce"
etiqueta_no <- "No produce"
valor_exito <- if (invertir_codificacion) 0 else 1
valor_fracaso <- if (invertir_codificacion) 1 else 0
exito <- ifelse(crudo == valor_exito, 1, ifelse(crudo == valor_fracaso, 0, NA))
exito <- exito[!is.na(exito)]
n_total <- length(exito)
p_hat <- sum(exito) / n_total # proporción muestral de "Sí produce"
cat("Registros brutos: ", n_bruto, "\n")
## Registros brutos: 47757
cat("Registros válidos depurados: ", n_total, "\n")
## Registros válidos depurados: 47757
cat("Registros descartados (NA): ", n_bruto - n_total, "\n")
## Registros descartados (NA): 0
cat("Categoría 'éxito': ", etiqueta_si, "\n")
## Categoría 'éxito': Sí produce
cat("Categoría 'fracaso': ", etiqueta_no, "\n")
## Categoría 'fracaso': No produce
La categoría de éxito es Sí produce y la de fracaso es No produce.
tdf <- data.frame(Categoria = factor(ifelse(exito == 1, etiqueta_si, etiqueta_no),
levels = c(etiqueta_si, etiqueta_no))) %>%
count(Categoria, name = "ni") %>%
mutate(hi = ni / sum(ni), Probabilidad = round(100 * hi, 2))
total <- data.frame(Categoria = "TOTAL", ni = sum(tdf$ni), hi = 1, Probabilidad = 100)
tabla_tdf <- bind_rows(tdf, total)
n_filas_tdf <- nrow(tabla_tdf)
tabla_tdf %>%
gt() %>%
tab_header(
title = md("**Tabla N\u00b01: Distribución de Frecuencias**"),
subtitle = md("*Produce petróleo, Kansas*")
) %>%
cols_label(Categoria = md("**Produce petróleo**"), ni = md("**Frecuencia (ni)**"),
hi = md("**Frecuencia relativa (hi)**"), Probabilidad = md("**Probabilidad (%)**")) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_title(groups = "title")) %>%
tab_style(style = cell_fill(color = "#F5F5F5"), locations = cells_body(rows = seq(1, n_filas_tdf, by = 2))) %>%
tab_style(style = list(cell_fill(color = "#D6D6D6"), cell_text(weight = "bold")),
locations = cells_body(rows = Categoria == "TOTAL")) %>%
fmt_number(columns = hi, decimals = 4) %>%
cols_align(align = "center", columns = c(ni, hi, Probabilidad)) %>%
cols_align(align = "left", columns = Categoria) %>%
tab_source_note(source_note = md("*Autor: GRUPO*")) %>%
tab_options(table.width = pct(90), table.font.size = px(13),
heading.title.font.size = px(16), heading.subtitle.font.size = px(12),
data_row.padding = px(6),
column_labels.border.top.width = px(2),
column_labels.border.bottom.width = px(2),
table_body.border.bottom.width = px(2),
table.border.top.style = "hidden", table.border.bottom.style = "hidden")
| Tabla N°1: Distribución de Frecuencias | |||
| Produce petróleo, Kansas | |||
| Produce petróleo | Frecuencia (ni) | Frecuencia relativa (hi) | Probabilidad (%) |
|---|---|---|---|
| Sí produce | 32640 | 0.6835 | 68.35 |
| No produce | 15117 | 0.3165 | 31.65 |
| TOTAL | 47757 | 1.0000 | 100.00 |
| Autor: GRUPO | |||
bp <- barplot(tdf$Probabilidad, names.arg = tdf$Categoria,
col = c("#1B9E77", "#D95F02"), border = NA,
main = "Gráfica 1. Distribución empírica de Produce petróleo",
ylab = "Probabilidad (%)",
ylim = c(0, max(tdf$Probabilidad) * 1.18))
text(bp, tdf$Probabilidad, labels = paste0(round(tdf$Probabilidad, 2), "%"), pos = 3)
Al ser PRODUCES_OIL una variable dicotómica (solo dos categorías), no tiene sentido agrupar en bloques artificiales para generar más categorías: el modelo natural es tratar cada registro como un ensayo de Bernoulli independiente,
\[X_i \sim \text{Bernoulli}(p), \quad i = 1, \dots, n\]
donde \(p\) es la proporción poblacional de “Sí produce”. A partir de la gráfica empírica (sección anterior) se conjetura un valor concreto de \(p\) —en vez de suponer arbitrariamente que las dos categorías son igualmente probables—, para luego comprobar estadísticamente si el modelo conjeturado es aceptable:
H0: p = 0.683 (valor conjeturado a partir de la
proporción observada).
H1: p ≠ 0.683.
p0 <- round(p_hat, 3) # proporción conjeturada (sección "Conjetura")
error_estandar <- sqrt(p_hat * (1 - p_hat) / n_total)
tabla_parametros <- data.frame(
Parametro = c("n (registros válidos)", "p\u0302 (proporción muestral)",
"p0 (proporción bajo H0)", "Error estándar de p\u0302"),
Valor = c(n_total, round(p_hat, 4), p0, round(error_estandar, 4))
)
tabla_parametros %>%
gt() %>%
tab_header(title = md("**Tabla N\u00b02: Parámetros del modelo Bernoulli(p) conjeturado**")) %>%
cols_label(Parametro = md("**Parámetro**"), Valor = md("**Valor**")) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_title(groups = "title")) %>%
tab_style(style = cell_fill(color = "#F5F5F5"), locations = cells_body(rows = seq(1, nrow(tabla_parametros), by = 2))) %>%
cols_align(align = "center", columns = Valor) %>%
cols_align(align = "left", columns = Parametro) %>%
tab_source_note(source_note = md("*Autor: GRUPO*")) %>%
tab_options(table.width = pct(90), table.font.size = px(13),
heading.title.font.size = px(16),
data_row.padding = px(6),
column_labels.border.top.width = px(2),
column_labels.border.bottom.width = px(2),
table_body.border.bottom.width = px(2),
table.border.top.style = "hidden", table.border.bottom.style = "hidden")
| Tabla N°2: Parámetros del modelo Bernoulli(p) conjeturado | |
| Parámetro | Valor |
|---|---|
| n (registros válidos) | 47757.0000 |
| p̂ (proporción muestral) | 0.6835 |
| p0 (proporción bajo H0) | 0.6830 |
| Error estándar de p̂ | 0.0021 |
| Autor: GRUPO | |
Solo existen dos categorías, así que la comparación entre lo observado y lo esperado bajo H0 se hace directamente sobre ellas (sin inventar categorías intermedias).
tabla_modelo <- data.frame(Categoria = c(etiqueta_si, etiqueta_no),
ni = c(sum(exito == 1), sum(exito == 0)))
tabla_modelo$hi_obs <- tabla_modelo$ni / n_total
tabla_modelo$p_teorica <- c(p0, 1 - p0)
tabla_modelo$ni_esperada <- tabla_modelo$p_teorica * n_total
tabla_modelo %>%
gt() %>%
tab_header(title = md(paste0("**Tabla N\u00b03: Frecuencias observadas frente al modelo Bernoulli(p0 = ", round(p0, 3), ")**"))) %>%
cols_label(Categoria = md("**Categoría**"), ni = md("**ni observada**"),
hi_obs = md("**hi observada**"), p_teorica = md("**P bajo H0**"),
ni_esperada = md("**ni esperada bajo H0**")) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_title(groups = "title")) %>%
tab_style(style = cell_fill(color = "#F5F5F5"), locations = cells_body(rows = seq(1, nrow(tabla_modelo), by = 2))) %>%
fmt_number(columns = c(hi_obs, p_teorica), decimals = 4) %>%
fmt_number(columns = ni_esperada, decimals = 2) %>%
cols_align(align = "center", columns = c(ni, hi_obs, p_teorica, ni_esperada)) %>%
cols_align(align = "left", columns = Categoria) %>%
tab_source_note(source_note = md("*Autor: GRUPO*")) %>%
tab_options(table.width = pct(90), table.font.size = px(13),
heading.title.font.size = px(16),
data_row.padding = px(6),
column_labels.border.top.width = px(2),
column_labels.border.bottom.width = px(2),
table_body.border.bottom.width = px(2),
table.border.top.style = "hidden", table.border.bottom.style = "hidden")
| Tabla N°3: Frecuencias observadas frente al modelo Bernoulli(p0 = 0.683) | ||||
| Categoría | ni observada | hi observada | P bajo H0 | ni esperada bajo H0 |
|---|---|---|---|---|
| Sí produce | 32640 | 0.6835 | 0.6830 | 32,618.03 |
| No produce | 15117 | 0.3165 | 0.3170 | 15,138.97 |
| Autor: GRUPO | ||||
comparacion <- rbind(Observado = tabla_modelo$hi_obs * 100, `Esperado bajo H0` = tabla_modelo$p_teorica * 100)
colnames(comparacion) <- tabla_modelo$Categoria
bp <- barplot(comparacion, beside = TRUE, col = c("#1B9E77", "#D95F02"), border = NA,
main = paste0("Gráfica 2. Realidad observada vs. modelo Bernoulli(p0 = ", round(p0, 3), ")"),
ylab = "Probabilidad (%)", ylim = c(0, 100), legend.text = TRUE,
args.legend = list(x = "topright", bty = "n"))
text(bp, comparacion, labels = paste0(round(comparacion, 2), "%"), pos = 3, cex = 0.8)
Con \(\hat p = 0.6835\), se calculan probabilidades para una muestra futura de 30 pozos, \(Y \sim \text{Binomial}(30, \hat p)\) (aquí sí es apropiado un modelo de conteo, porque se trata de una muestra nueva de varios pozos, no de la misma variable dicotómica original):
n_futuro <- 30
media_futura <- n_futuro * p_hat
p_al_menos_mitad <- 1 - pbinom(floor(n_futuro / 2) - 1, size = n_futuro, prob = p_hat)
p_maximo_25pct <- pbinom(floor(0.25 * n_futuro), size = n_futuro, prob = p_hat)
tabla_probabilidades <- data.frame(
Evento = c(paste0("E[Y] esperado en n = ", n_futuro, " pozos"),
"P(Y \u2265 mitad de la muestra)",
"P(Y \u2264 25% de la muestra)"),
Valor = round(c(media_futura, p_al_menos_mitad, p_maximo_25pct), 4)
)
tabla_probabilidades %>%
gt() %>%
tab_header(title = md(paste0("**Tabla N\u00b05: Probabilidades bajo el modelo Binomial(30, p\u0302 = ", round(p_hat, 4), ")**")),
subtitle = md("*Para una muestra futura*")) %>%
cols_label(Evento = md("**Evento**"), Valor = md("**Valor**")) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_title(groups = "title")) %>%
tab_style(style = cell_fill(color = "#F5F5F5"), locations = cells_body(rows = seq(1, nrow(tabla_probabilidades), by = 2))) %>%
cols_align(align = "center", columns = Valor) %>%
cols_align(align = "left", columns = Evento) %>%
tab_source_note(source_note = md("*Autor: GRUPO*")) %>%
tab_options(table.width = pct(90), table.font.size = px(13),
heading.title.font.size = px(16), heading.subtitle.font.size = px(12),
data_row.padding = px(6),
column_labels.border.top.width = px(2),
column_labels.border.bottom.width = px(2),
table_body.border.bottom.width = px(2),
table.border.top.style = "hidden", table.border.bottom.style = "hidden")
| Tabla N°5: Probabilidades bajo el modelo Binomial(30, p̂ = 0.6835) | |
| Para una muestra futura | |
| Evento | Valor |
|---|---|
| E[Y] esperado en n = 30 pozos | 20.5038 |
| P(Y ≥ mitad de la muestra) | 0.9888 |
| P(Y ≤ 25% de la muestra) | 0.0000 |
| Autor: GRUPO | |
Intervalo de confianza al 95% para la proporción poblacional \(p\), calculado por los métodos de Wald y de Wilson, a partir de \(\hat p = 0.6835\) y \(n = 47757\).
z <- qnorm(0.975)
wald_li <- p_hat - z * sqrt(p_hat * (1 - p_hat) / n_total)
wald_ls <- p_hat + z * sqrt(p_hat * (1 - p_hat) / n_total)
denom <- 1 + z^2 / n_total
centro <- (p_hat + z^2 / (2 * n_total)) / denom
margen <- (z / denom) * sqrt((p_hat * (1 - p_hat) / n_total) + (z^2 / (4 * n_total^2)))
wilson_li <- centro - margen
wilson_ls <- centro + margen
tabla_ic <- data.frame(Metodo = c("Wald", "Wilson"),
Lim_Inferior = round(c(wald_li, wilson_li), 4),
Lim_Superior = round(c(wald_ls, wilson_ls), 4))
tabla_ic %>%
gt() %>%
tab_header(title = md("**Tabla N\u00b06: Intervalo de Confianza al 95% para p**")) %>%
cols_label(Metodo = md("**Método**"), Lim_Inferior = md("**Límite inferior**"),
Lim_Superior = md("**Límite superior**")) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()) %>%
tab_style(style = list(cell_fill(color = "#2C2C2C"), cell_text(color = "white", weight = "bold")),
locations = cells_title(groups = "title")) %>%
tab_style(style = cell_fill(color = "#F5F5F5"), locations = cells_body(rows = seq(1, nrow(tabla_ic), by = 2))) %>%
cols_align(align = "center", columns = c(Lim_Inferior, Lim_Superior)) %>%
cols_align(align = "left", columns = Metodo) %>%
tab_source_note(source_note = md("*Autor: GRUPO*")) %>%
tab_options(table.width = pct(90), table.font.size = px(13),
heading.title.font.size = px(16),
data_row.padding = px(6),
column_labels.border.top.width = px(2),
column_labels.border.bottom.width = px(2),
table_body.border.bottom.width = px(2),
table.border.top.style = "hidden", table.border.bottom.style = "hidden")
| Tabla N°6: Intervalo de Confianza al 95% para p | ||
| Método | Límite inferior | Límite superior |
|---|---|---|
| Wald | 0.6793 | 0.6876 |
| Wilson | 0.6793 | 0.6876 |
| Autor: GRUPO | ||
De los 47.757 registros válidos, la proporción estimada de Sí produce es \(\hat p = 0.6835\) (68.35%). La prueba Z para una proporción, contrastada contra el valor conjeturado \(p_0 = 0.683\), arrojó un estadístico de 0.216 con valor p de 0.8289, por lo que al 5% de significancia se acepta el modelo Bernoulli(p0) conjeturado. El intervalo de confianza al 95% para \(p\) (Wilson) es [0.6793, 0.6876].
Autor: GRUPO — Análisis Estadístico, Kansas Hydrocarbon Leases Dataset