1 Carga de Librerías

library(tidyverse)
library(knitr)
library(kableExtra)
library(broom)

2 Carga de Datos

archivo <- "oil_and_gas_leases_data (2)(2).csv"
if(!file.exists(archivo)) archivo <- "oil_and_gas_leases_data.csv"
if(!file.exists(archivo)) archivo <- "Kansas_limpio.csv"

datos <- read.csv(archivo)
str(datos)

## 'data.frame':    47757 obs. of  24 variables:
##  $ KID                      : int  1001106903 1001106572 1001106590 1001107343 1001108234 1001106684 1001107377 1001107386 1001107740 1001106710 ...
##  $ DEPTH_OF_WELL            : num  700 800 1400 1125 2940 ...
##  $ CUMULATIVE_PRODUCTION    : num  47225 275063 82624 7544 681006 ...
##  $ AVG_PRODUCTION           : num  859 5001 1758 377 24322 ...
##  $ LATITUDE                 : num  37.1 38.8 37.5 37.8 37.1 ...
##  $ LONGITUDE                : num  -95.9 -95.2 -96.3 -95.7 -101.3 ...
##  $ YEARS_ACTIVE             : num  55 55 47 20 28 55 20 48 48 55 ...
##  $ SECTION                  : num  33 11 34 8 30 4 26 28 11 17 ...
##  $ COUNTY_CODE              : num  125 45 49 207 189 121 49 1 31 121 ...
##  $ STATE_CODE               : int  15 15 15 15 15 15 15 15 15 15 ...
##  $ TOWNSHIP                 : num  33 15 29 26 33 17 30 26 23 16 ...
##  $ RANGE                    : num  14 20 10 16 36 25 12 21 16 24 ...
##  $ PRODUCES_OIL             : num  1 1 1 1 0 1 1 1 1 1 ...
##  $ PRODUCES_GAS             : num  0 0 0 0 1 0 0 0 0 0 ...
##  $ OPERATOR_NAME            : chr  "Horton, John" "Whitlow Energy, Inc." "Suerte Oil Company" "Patterson-Blackford" ...
##  $ FIELD_NAME               : chr  "WAYSIDE-HAVANA" "BALDWIN" "DUNKLEBERGER" "ROSE EAST" ...
##  $ PRODUCING_FORMATION      : chr  "UNKNOWN" "UNKNOWN" "UNKNOWN" "UNKNOWN" ...
##  $ LONGITUDE_LATITUDE_SOURCE: chr  "CENTER_OF_SECTION" "CENTER_OF_SECTION" "CENTER_OF_SECTION" "CENTER_OF_SECTION" ...
##  $ PROD_LEVEL               : chr  "MEDIUM" "HIGH" "MEDIUM" "LOW" ...
##  $ DEPTH_LEVEL              : chr  "SHALLOW" "SHALLOW" "SHALLOW" "SHALLOW" ...
##  $ LIFE_STAGE               : chr  "OLD" "OLD" "OLD" "MATURE" ...
##  $ AVG_PROD_LEVEL           : chr  "LOW" "MEDIUM" "MEDIUM" "LOW" ...
##  $ TOWNSHIP_DIRECTION       : chr  "S" "S" "S" "S" ...
##  $ RANGE_DIRECTION          : chr  "E" "E" "E" "E" ...

3 Extracción y Depuración de la Variable

Se trabaja con la relación entre:

Variable independiente: YEARS_ACTIVE
Variable dependiente: AVG_PRODUCTION

Se analiza una sección de 1 a 20 años activos, porque en esta zona la producción promedio presenta un descenso más evidente.

datos_modelo <- datos %>%
  transmute(
    X = suppressWarnings(as.numeric(.data[["YEARS_ACTIVE"]])),
    Y = suppressWarnings(as.numeric(.data[["AVG_PRODUCTION"]]))
  ) %>%
  filter(!is.na(X), !is.na(Y), X > 0, Y > 0, X >= 1, X <= 20)

n_total <- nrow(datos_modelo)
n_total

## [1] 26246

4 Tabla de Distribución de Frecuencias

Se aplica la regla de Sturges con un máximo de 10 intervalos:

\[ k = 1 + 3.322\log_{10}(n) \]

k_sturges <- ceiling(1 + 3.322 * log10(n_total))
k <- min(k_sturges, 10)
k

## [1] 10

breaks <- pretty(range(datos_modelo$X), n = k)
if(length(breaks) - 1 > 10) breaks <- seq(min(datos_modelo$X), max(datos_modelo$X), length.out = 11)

resumen_modelo <- datos_modelo %>%
  mutate(intervalo = cut(X, breaks = breaks, include.lowest = TRUE, right = TRUE)) %>%
  group_by(intervalo) %>%
  summarise(
    MC = mean(X),
    ni = n(),
    Y_MEDIA = mean(Y),
    .groups = "drop"
  ) %>%
  mutate(
    hi = round(100 * ni / sum(ni), 2),
    Ni = cumsum(ni),
    Hi = round(cumsum(hi), 2)
  ) %>%
  filter(!is.na(intervalo))

TDF <- resumen_modelo %>%
  select(Intervalo = intervalo, MC, ni, hi, Ni, Hi)

TDF_total <- bind_rows(
  TDF,
  tibble(
    Intervalo = "TOTAL",
    MC = NA_real_,
    ni = sum(TDF$ni),
    hi = 100,
    Ni = NA_integer_,
    Hi = 100
  )
)

kable(TDF_total, caption = "Tabla N°1: Distribución de frecuencias") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Tabla N°1: Distribución de frecuencias
Intervalo	MC	ni	hi	Ni	Hi
[0,2]	1.330070	3269	12.46	3269	12.46
(2,4]	3.467194	2149	8.19	5418	20.65
(4,6]	5.521918	1939	7.39	7357	28.04
(6,8]	7.504696	2023	7.71	9380	35.75
(8,10]	9.570470	2235	8.52	11615	44.27
(10,12]	11.467408	3237	12.33	14852	56.60
(12,14]	13.500176	2847	10.85	17699	67.45
(14,16]	15.467043	2655	10.12	20354	77.57
(16,18]	17.470515	3222	12.28	23576	89.85
(18,20]	19.398876	2670	10.17	26246	100.02
TOTAL	NA	26246	100.00	NA	100.00

5 Gráfica

ggplot(datos_modelo, aes(x = X, y = Y)) +
  geom_point(color = "gray35", alpha = 0.35, size = 1) +
  labs(
    title = "Gráfica N°1: Relación observada entre YEARS_ACTIVE y AVG_PRODUCTION",
    x = "YEARS_ACTIVE",
    y = "AVG_PRODUCTION"
  ) +
  theme_bw()

6 Conjetura

La nube de puntos muestra que la producción promedio disminuye con rapidez durante los primeros años de actividad del pozo y luego tiende a estabilizarse. Por ello, se propone un modelo exponencial decreciente.

7 Cálculo de Parámetros

modelo <- lm(log(Y_MEDIA) ~ MC, data = resumen_modelo)
resumen_modelo$Y_MODELO <- exp(predict(modelo, newdata = resumen_modelo))
a <- exp(coef(modelo)[1])
b <- coef(modelo)[2]
ecuacion <- paste0("y = ", round(a, 4), " e^(", round(b, 4), "x)")
parametros <- data.frame(
  Parametro = c("a = exp(intercepto)", "b"),
  Valor = round(c(a, b), 6)
)


parametros %>%
  kable(caption = "Tabla N°2: Parámetros estimados del modelo") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Tabla N°2: Parámetros estimados del modelo
	Parametro	Valor
(Intercept)	a = exp(intercepto)	7768.382171
MC	b	-0.013235

ecuacion

## [1] "y = 7768.3822 e^(-0.0132x)"

Gráfica: realidad vs modelo

ggplot(resumen_modelo, aes(x = MC)) +
  geom_point(aes(y = Y_MEDIA), color = "gray20", size = 3) +
  geom_line(aes(y = Y_MEDIA), color = "gray45", linewidth = 0.8) +
  geom_point(aes(y = Y_MODELO), color = "black", shape = 17, size = 3) +
  geom_line(aes(y = Y_MODELO), color = "black", linewidth = 0.9, linetype = "dashed") +
  labs(
    title = "Gráfica N°2: Realidad vs Modelo Exponencial",
    x = "Marca de clase de YEARS_ACTIVE",
    y = "Media de AVG_PRODUCTION",
    caption = "Puntos grises: realidad observada | Línea negra punteada: modelo ajustado"
  ) +
  theme_bw()

Test de Pearson

obs_pct <- 100 * resumen_modelo$Y_MEDIA / sum(resumen_modelo$Y_MEDIA)
esp_pct <- 100 * resumen_modelo$Y_MODELO / sum(resumen_modelo$Y_MODELO)

pearson_r <- cor(obs_pct, esp_pct)
pearson_porcentaje <- pearson_r * 100
pearson_porcentaje

## [1] 61.16926

Test Chi-cuadrado

chi_calculado <- sum((obs_pct - esp_pct)^2 / esp_pct)
numero_parametros <- length(coef(modelo))
gl <- length(obs_pct) - numero_parametros - 1
chi_critico <- qchisq(0.95, df = gl)
p_valor <- 1 - pchisq(chi_calculado, df = gl)
decision <- ifelse(chi_calculado < chi_critico & p_valor > 0.05,
                   "No se rechaza H0: modelo aceptado",
                   "Se rechaza H0: modelo no aceptado")

resumen_bondad <- tibble(
  Modelo = "Exponencial",
  `Pearson (%)` = round(pearson_porcentaje, 2),
  `Chi-cuadrado calculado` = round(chi_calculado, 4),
  `gl` = gl,
  `Chi-cuadrado crítico` = round(chi_critico, 4),
  `p-valor` = round(p_valor, 4),
  Decisión = decision
)

kable(resumen_bondad, caption = "Tabla N°3: Resumen de bondad del ajuste") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Tabla N°3: Resumen de bondad del ajuste
Modelo	Pearson (%)	Chi-cuadrado calculado	gl	Chi-cuadrado crítico	p-valor	Decisión
Exponencial	61.17	1.1485	7	14.0671	0.9921	No se rechaza H0: modelo aceptado

8 Cálculo de Probabilidades

Se estima la probabilidad relativa de que la producción promedio modelada se ubique en cada intervalo de años activos.

probabilidades <- resumen_modelo %>%
  mutate(
    Prob_observada = round(obs_pct / 100, 4),
    Prob_modelo = round(esp_pct / 100, 4)
  ) %>%
  select(Intervalo = intervalo, MC, Prob_observada, Prob_modelo)

probabilidades_total <- bind_rows(
  probabilidades,
  tibble(
    Intervalo = "TOTAL",
    MC = NA_real_,
    Prob_observada = sum(probabilidades$Prob_observada),
    Prob_modelo = sum(probabilidades$Prob_modelo)
  )
)

kable(probabilidades_total, caption = "Tabla N°4: Probabilidades observadas y modeladas por intervalo") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Tabla N°4: Probabilidades observadas y modeladas por intervalo
Intervalo	MC	Prob_observada	Prob_modelo
[0,2]	1.330070	0.1169	0.1125
(2,4]	3.467194	0.1268	0.1094
(4,6]	5.521918	0.1047	0.1065
(6,8]	7.504696	0.0903	0.1037
(8,10]	9.570470	0.0952	0.1009
(10,12]	11.467408	0.1011	0.0984
(12,14]	13.500176	0.0796	0.0958
(14,16]	15.467043	0.0849	0.0933
(16,18]	17.470515	0.0973	0.0909
(18,20]	19.398876	0.1033	0.0886
TOTAL	NA	1.0001	1.0000

9 Intervalo de Confianza

Se calcula el intervalo de confianza del 95 % para los parámetros del modelo ajustado.

IC_parametros <- confint(modelo, level = 0.95) %>%
  as.data.frame() %>%
  rownames_to_column("Parametro") %>%
  rename(`Límite inferior 95%` = `2.5 %`, `Límite superior 95%` = `97.5 %`)

kable(IC_parametros, caption = "Tabla N°5: Intervalos de confianza del 95 % para los parámetros") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Tabla N°5: Intervalos de confianza del 95 % para los parámetros
Parametro	Límite inferior 95%	Límite superior 95%
(Intercept)	8.7738798	9.1417546
MC	-0.0286341	0.0021645

10 Conclusiones

Se construyó la relación entre YEARS_ACTIVE y AVG_PRODUCTION usando el modelo Exponencial.
La conjetura se planteó a partir de la forma observada en la nube de puntos y se trabajó por sección para evitar forzar el modelo sobre toda la dispersión original.
El ajuste se evaluó mediante el coeficiente de Pearson y el estadístico Chi-cuadrado, tomando como criterio la no contradicción entre los valores observados y los valores esperados por el modelo.
El modelo permite representar la tendencia principal de la variable dependiente dentro del intervalo de análisis seleccionado.

Years Active vs Avg Production

Leslye Quinchiguango

2026-06-22