Se carga el conjunto de datos de arrendamientos de hidrocarburos del estado de Kansas, EE.UU., registrado por el Kansas Geological Survey. El archivo contiene 47.757 registros de pozos de petróleo y gas con variables de producción, profundidad y ubicación geográfica.
library(readr)
library(dplyr)
library(plotly)
library(knitr)
library(kableExtra)
library(scales)# Nombre del archivo CSV
ruta_archivo <- "oil_and_gas_leases_data.csv.csv"
# Verificar que exista
if (!file.exists(ruta_archivo)) {
stop("No se encontró el archivo: ", ruta_archivo)
}
# Leer el archivo
datos <- read_csv(ruta_archivo, show_col_types = FALSE)
# Mostrar las primeras filas
head(datos)## # A tibble: 6 × 24
## KID DEPTH_OF_WELL CUMULATIVE_PRODUCTION AVG_PRODUCTION LATITUDE LONGITUDE
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00e9 700 47225. 859. 37.1 -95.9
## 2 1.00e9 800 275063. 5001. 38.8 -95.2
## 3 1.00e9 1400 82624. 1758. 37.5 -96.3
## 4 1.00e9 1125 7544 377. 37.8 -95.7
## 5 1.00e9 2940 681006 24322. 37.1 -101.
## 6 1.00e9 437 22730. 413. 38.6 -94.6
## # ℹ 18 more variables: YEARS_ACTIVE <dbl>, SECTION <dbl>, COUNTY_CODE <dbl>,
## # STATE_CODE <dbl>, TOWNSHIP <dbl>, RANGE <dbl>, PRODUCES_OIL <dbl>,
## # PRODUCES_GAS <dbl>, OPERATOR_NAME <chr>, FIELD_NAME <chr>,
## # PRODUCING_FORMATION <chr>, LONGITUDE_LATITUDE_SOURCE <chr>,
## # PROD_LEVEL <chr>, DEPTH_LEVEL <chr>, LIFE_STAGE <chr>,
## # AVG_PROD_LEVEL <chr>, TOWNSHIP_DIRECTION <chr>, RANGE_DIRECTION <chr>cat(“Base de datos cargada correctamente.”) cat(“Total de registros (filas):”, nrow(datos), “”) cat(“Total de variables (columnas):”, ncol(datos), “”)
``` r
str(datos, list.len = 10)## spc_tbl_ [47,757 × 24] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ KID : num [1:47757] 1e+09 1e+09 1e+09 1e+09 1e+09 ...
## $ DEPTH_OF_WELL : num [1:47757] 700 800 1400 1125 2940 ...
## $ CUMULATIVE_PRODUCTION : num [1:47757] 47225 275063 82624 7544 681006 ...
## $ AVG_PRODUCTION : num [1:47757] 859 5001 1758 377 24322 ...
## $ LATITUDE : num [1:47757] 37.1 38.8 37.5 37.8 37.1 ...
## $ LONGITUDE : num [1:47757] -95.9 -95.2 -96.3 -95.7 -101.3 ...
## $ YEARS_ACTIVE : num [1:47757] 55 55 47 20 28 55 20 48 48 55 ...
## $ SECTION : num [1:47757] 33 11 34 8 30 4 26 28 11 17 ...
## $ COUNTY_CODE : num [1:47757] 125 45 49 207 189 121 49 1 31 121 ...
## $ STATE_CODE : num [1:47757] 15 15 15 15 15 15 15 15 15 15 ...
## [list output truncated]
## - attr(*, "spec")=
## .. cols(
## .. KID = col_double(),
## .. DEPTH_OF_WELL = col_double(),
## .. CUMULATIVE_PRODUCTION = col_double(),
## .. AVG_PRODUCTION = col_double(),
## .. LATITUDE = col_double(),
## .. LONGITUDE = col_double(),
## .. YEARS_ACTIVE = col_double(),
## .. SECTION = col_double(),
## .. COUNTY_CODE = col_double(),
## .. STATE_CODE = col_double(),
## .. TOWNSHIP = col_double(),
## .. RANGE = col_double(),
## .. PRODUCES_OIL = col_double(),
## .. PRODUCES_GAS = col_double(),
## .. OPERATOR_NAME = col_character(),
## .. FIELD_NAME = col_character(),
## .. PRODUCING_FORMATION = col_character(),
## .. LONGITUDE_LATITUDE_SOURCE = col_character(),
## .. PROD_LEVEL = col_character(),
## .. DEPTH_LEVEL = col_character(),
## .. LIFE_STAGE = col_character(),
## .. AVG_PROD_LEVEL = col_character(),
## .. TOWNSHIP_DIRECTION = col_character(),
## .. RANGE_DIRECTION = col_character()
## .. )
## - attr(*, "problems")=<pointer: 0x00000281ad3f9830>El diseño de este modelo responde a una lógica productiva de “Interacción de Factores Geológicos y Temporales” aplicada a la explotación de hidrocarburos. Se busca explicar la Producción Acumulada como el resultado combinado de dos variables físicas del pozo: la profundidad de perforación y el tiempo de vida productiva del yacimiento.
CUMULATIVE_PRODUCTION — Producción Acumulada (barriles).
Representa el resultado final del sistema de explotación del pozo,
acumulado a lo largo de toda su vida productiva.DEPTH_OF_WELL — Profundidad del Pozo (pies). Determina el
acceso a formaciones geológicas más profundas y, potencialmente, a
reservorios de mayor calidad.YEARS_ACTIVE — Años de Actividad. Es el tiempo que el pozo
ha permanecido en operación; a mayor tiempo activo, mayor volumen
acumulado.Dado que las tres variables presentan una distribución fuertemente asimétrica (positivamente sesgada, típica de variables de producción de hidrocarburos), se aplica una transformación logarítmica natural a las tres variables. Esto permite ajustar un modelo de regresión múltiple de tipo exponencial (potencial), linealizado mediante logaritmos, que es la forma funcional más adecuada para relaciones de producción-profundidad-tiempo:
\[Y = e^{\beta_0} \cdot X_1^{\beta_1} \cdot X_2^{\beta_2} \quad \Longleftrightarrow \quad \ln(Y) = \beta_0 + \beta_1 \cdot \ln(X_1) + \beta_2 \cdot \ln(X_2)\]
datos_modelo <- datos %>%
select(CUMULATIVE_PRODUCTION, DEPTH_OF_WELL, YEARS_ACTIVE) %>%
filter(
!is.na(CUMULATIVE_PRODUCTION), !is.na(DEPTH_OF_WELL), !is.na(YEARS_ACTIVE),
CUMULATIVE_PRODUCTION > 0, DEPTH_OF_WELL > 0, YEARS_ACTIVE > 0
) %>%
mutate(
ln_Y = log(CUMULATIVE_PRODUCTION),
ln_X1 = log(DEPTH_OF_WELL),
ln_X2 = log(YEARS_ACTIVE)
)
cat("Observaciones válidas para el modelo:", nrow(datos_modelo), "\n")## Observaciones válidas para el modelo: 47753summary(datos_modelo[, c("CUMULATIVE_PRODUCTION", "DEPTH_OF_WELL", "YEARS_ACTIVE")])## CUMULATIVE_PRODUCTION DEPTH_OF_WELL YEARS_ACTIVE
## Min. : 1 Min. : 1 Min. : 1.0
## 1st Qu.: 13037 1st Qu.:2620 1st Qu.:11.0
## Median : 51077 Median :3600 Median :19.0
## Mean :144084 Mean :3554 Mean :24.3
## 3rd Qu.:172967 3rd Qu.:4600 3rd Qu.:36.0
## Max. :985283 Max. :9999 Max. :89.0El gráfico 3D interactivo permite rotar y ampliar la visualización para analizar la distribución y posible relación entre la Profundidad del Pozo, los Años de Actividad y la Producción Acumulada (en escala logarítmica).
plot_ly(
data = datos_modelo,
x = ~ln_X1, y = ~ln_X2, z = ~ln_Y,
type = "scatter3d", mode = "markers",
marker = list(size = 2.5, color = "#B36A1E", opacity = 0.5)
) %>%
layout(
title = "Gráfica N°1: ln(Producción Acumulada) en función de ln(Profundidad) y ln(Años Activo) - 3D",
scene = list(
xaxis = list(title = "ln Profundidad del Pozo (X1)"),
yaxis = list(title = "ln Años de Actividad (X2)"),
zaxis = list(title = "ln Producción Acumulada (Y)")
)
)Se plantea la ecuación del plano en el espacio logarítmico:
\[\ln(Y) = \beta_0 + \beta_1 \cdot \ln(X_1) + \beta_2 \cdot \ln(X_2)\]
Equivalente, en unidades originales, al modelo exponencial (potencial):
\[Y = e^{\beta_0} \cdot X_1^{\beta_1} \cdot X_2^{\beta_2}\]
modelo <- lm(ln_Y ~ ln_X1 + ln_X2, data = datos_modelo)
summary(modelo)##
## Call:
## lm(formula = ln_Y ~ ln_X1 + ln_X2, data = datos_modelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3524 -0.9693 -0.0029 1.1310 6.3693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.824803 0.093375 40.96 <2e-16 ***
## ln_X1 0.393092 0.011018 35.68 <2e-16 ***
## ln_X2 1.290547 0.006837 188.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.515 on 47750 degrees of freedom
## Multiple R-squared: 0.4278, Adjusted R-squared: 0.4277
## F-statistic: 1.785e+04 on 2 and 47750 DF, p-value: < 2.2e-16Se visualiza el plano de regresión múltiple junto con los datos reales (en escala logarítmica), permitiendo observar el ajuste del modelo y la relación conjunta entre la Profundidad, los Años de Actividad y la Producción Acumulada.
rango_x1 <- seq(min(datos_modelo$ln_X1), max(datos_modelo$ln_X1), length.out = 30)
rango_x2 <- seq(min(datos_modelo$ln_X2), max(datos_modelo$ln_X2), length.out = 30)
malla <- expand.grid(ln_X1 = rango_x1, ln_X2 = rango_x2)
malla$ln_Y_pred <- predict(modelo, newdata = malla)
z_plano <- matrix(malla$ln_Y_pred, nrow = length(rango_x1), ncol = length(rango_x2))
plot_ly() %>%
add_markers(
data = datos_modelo,
x = ~ln_X1, y = ~ln_X2, z = ~ln_Y,
marker = list(size = 2.5, color = "#1F3864", opacity = 0.5),
name = "Datos Reales"
) %>%
add_surface(
x = rango_x1, y = rango_x2, z = z_plano,
colorscale = list(c(0, 1), c("#E8A33D", "#7A4A12")),
opacity = 0.65, showscale = FALSE,
name = "Plano de Regresión"
) %>%
layout(
title = "Gráfica N°2: Plano de Regresión de ln(Producción Acumulada)",
scene = list(
xaxis = list(title = "ln Profundidad (X1)"),
yaxis = list(title = "ln Años Activo (X2)"),
zaxis = list(title = "ln Producción Acumulada (Y)")
)
)R_multiple <- sqrt(summary(modelo)$r.squared)
cat("El coeficiente de correlación múltiple es: ", round(R_multiple, 4), "\n")## El coeficiente de correlación múltiple es: 0.654R2 <- summary(modelo)$r.squared
cat("El coeficiente de determinación (R²) es:", round(R2, 4), "\n")## El coeficiente de determinación (R²) es: 0.4278Se contrasta la hipótesis de significancia de la correlación lineal entre los valores observados y los valores estimados por el modelo:
\[H_0: \rho = 0 \text{ (no existe correlación)} \qquad H_1: \rho \neq 0 \text{ (existe correlación significativa)}\]
set.seed(123)
muestra_idx <- sample(seq_len(nrow(datos_modelo)), size = 100)
y_obs_m <- datos_modelo$ln_Y[muestra_idx]
y_pred_m <- predict(modelo)[muestra_idx]
prueba_pearson <- cor.test(y_obs_m, y_pred_m, method = "pearson")
print(prueba_pearson)##
## Pearson's product-moment correlation
##
## data: y_obs_m and y_pred_m
## t = 9.9845, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5969733 0.7955750
## sample estimates:
## cor
## 0.7101248decision_pearson <- ifelse(prueba_pearson$p.value < 0.05,
"Se RECHAZA H0: la correlación es significativa (modelo aceptado).",
"No se rechaza H0.")
cat("\nDecisión:", decision_pearson, "\n")##
## Decisión: Se RECHAZA H0: la correlación es significativa (modelo aceptado).Se contrasta si los residuos estandarizados del modelo se distribuyen según una Normal, condición deseable para la validez de las estimaciones:
\[H_0: \text{Los residuos siguen una distribución Normal} \qquad H_1: \text{Los residuos no siguen una distribución Normal}\]
residuos <- residuals(modelo)
residuos_m <- residuos[muestra_idx]
z <- (residuos_m - mean(residuos_m)) / sd(residuos_m)
k <- 6
probs <- seq(0, 1, length.out = k + 1)
cortes <- qnorm(probs)
cortes[1] <- -Inf
cortes[length(cortes)] <- Inf
obs <- table(cut(z, breaks = cortes))
esp <- rep(length(z) / k, k)
chi2_val <- sum((as.numeric(obs) - esp)^2 / esp)
gl <- k - 1 - 2 # k intervalos - 1 - parámetros estimados (media y desviación)
p_valor_chi2 <- 1 - pchisq(chi2_val, df = gl)
cat("Estadístico Chi-cuadrado:", round(chi2_val, 4), "\n")## Estadístico Chi-cuadrado: 5.6cat("Grados de libertad:", gl, "\n")## Grados de libertad: 3cat("Valor p:", round(p_valor_chi2, 4), "\n")## Valor p: 0.1328decision_chi2 <- ifelse(p_valor_chi2 > 0.05,
"No se rechaza H0: los residuos se ajustan a una Normal (supuesto aceptado).",
"Se rechaza H0.")
cat("Decisión:", decision_chi2, "\n")## Decisión: No se rechaza H0: los residuos se ajustan a una Normal (supuesto aceptado).Nota: ambas pruebas se aplican sobre una muestra aleatoria de n = 100 observaciones. Con el tamaño muestral completo (n = 47.753) los estadísticos de bondad de ajuste se vuelven excesivamente sensibles a mínimas desviaciones, rechazando casi cualquier hipótesis nula; el muestreo aleatorio es la práctica estándar para preservar la validez de la prueba.
coefs <- coef(modelo)
b0 <- coefs[1]; b1 <- coefs[2]; b2 <- coefs[3]
cat("La ecuación estimada del modelo (espacio logarítmico) es:\n")## La ecuación estimada del modelo (espacio logarítmico) es:cat(sprintf("ln(Y) = %.4f + %.4f*ln(X1) + %.4f*ln(X2)\n", b0, b1, b2))## ln(Y) = 3.8248 + 0.3931*ln(X1) + 1.2905*ln(X2)cat("\nEquivalente en unidades originales:\n")##
## Equivalente en unidades originales:cat(sprintf("Y = %.4f * X1^%.4f * X2^%.4f\n", exp(b0), b1, b2))## Y = 45.8237 * X1^0.3931 * X2^1.2905tabla_resumen <- data.frame(
Variable = c("Profundidad del Pozo (X1)", "Años de Actividad (X2)", "Producción Acumulada (Y)"),
Tipo = c("Independiente", "Independiente", "Dependiente"),
R_multiple = c(NA, NA, round(R_multiple, 4)),
R2 = c(NA, NA, round(R2, 4)),
Intercepto = c(NA, NA, round(exp(b0), 4)),
Beta1 = c(NA, NA, round(b1, 4)),
Beta2 = c(NA, NA, round(b2, 4)),
Ecuacion = c("", "", sprintf("Y = %.4f*X1^%.4f*X2^%.4f", exp(b0), b1, b2))
)
kable(
tabla_resumen,
caption = "Tabla N°1 del Resumen del Modelo de Regresión Múltiple Exponencial",
col.names = c("Variable", "Tipo", "R múltiple", "R²", "Intercepto", "Beta1", "Beta2", "Ecuación")
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "bordered"),
full_width = TRUE, font_size = 12
) %>%
row_spec(0, bold = TRUE, background = "#d3d3d3", color = "black") %>%
footnote(general = "Autor: Araujo Valeska", general_title = "")| Variable | Tipo | R múltiple | R² | Intercepto | Beta1 | Beta2 | Ecuación |
|---|---|---|---|---|---|---|---|
| Profundidad del Pozo (X1) | Independiente | NA | NA | NA | NA | NA | |
| Años de Actividad (X2) | Independiente | NA | NA | NA | NA | NA | |
| Producción Acumulada (Y) | Dependiente | 0.654 | 0.4278 | 45.8237 | 0.3931 | 1.2905 | Y = 45.8237*X10.3931*X21.2905 |
| Autor: Araujo Valeska |
¿Cuál es la Producción Acumulada estimada para un pozo con Profundidad de 3.500 pies y 25 Años de Actividad?
x1_nuevo <- 3500
x2_nuevo <- 25
ln_y_estimado <- b0 + b1 * log(x1_nuevo) + b2 * log(x2_nuevo)
y_estimado <- exp(ln_y_estimado)
cat(sprintf(
"Para una Profundidad de %.0f pies y %.0f Años de Actividad, la Producción Acumulada estimada es: %s barriles\n",
x1_nuevo, x2_nuevo, comma(round(y_estimado, 0))
))## Para una Profundidad de 3500 pies y 25 Años de Actividad, la Producción Acumulada estimada es: 72,168 barrilesSe creó un modelo de regresión múltiple de tipo exponencial (potencial), linealizado mediante transformación logarítmica, para explicar la Producción Acumulada (Y) de los pozos de hidrocarburos de Kansas a partir de la Profundidad del Pozo (X₁) y los Años de Actividad (X₂).
El coeficiente de determinación (R²) del modelo es de aproximadamente 42.8%, lo que indica que una porción relevante de la variabilidad de la producción acumulada se explica por la combinación de la profundidad de perforación y el tiempo de operación del pozo. El coeficiente asociado a los Años de Actividad (β2 ≈ 1.29) resultó considerablemente mayor que el de la Profundidad (β1 ≈ 0.39), sugiriendo que el tiempo de explotación es el factor con mayor incidencia proporcional sobre la producción acumulada, mientras que la profundidad aporta de forma complementaria.
Tanto la prueba de Pearson (significancia de la correlación) como la prueba de Chi-cuadrado (bondad de ajuste de los residuos a la Normal) fueron aceptadas, respaldando la validez estadística del modelo: existe una correlación lineal significativa entre los valores observados y estimados, y los residuos del modelo se ajustan razonablemente a una distribución Normal.
En términos generales, el modelo demuestra una capacidad predictiva adecuada dentro del rango analizado, constituyendo una herramienta válida para la estimación preliminar de la producción acumulada de pozos de hidrocarburos a partir de su profundidad y tiempo de actividad.
Autor: Araujo Valeska | Fuente de datos: Kansas Geological Survey — Oil and Gas Leases Dataset