library(dplyr)
library(ggplot2)
library(gt)
library(stringr)
ruta_csv <- "C:/Users/PATRICIA/Desktop/pr-estadistica/Oil__Gas____Other_Regulated_Wells__Beginning_1860 (3).csv"
Datos <- read.csv(ruta_csv,
header = TRUE,
sep = ";",
dec = ".",
fileEncoding = "Latin1",
stringsAsFactors = FALSE)
cat("Dimensiones del dataset:", nrow(Datos), "filas x", ncol(Datos), "columnas\n")
## Dimensiones del dataset: 47407 filas x 55 columnas
El TVD (True Vertical Depth, ft) es la variable independiente (x): una magnitud física definida previamente, sin depender del costo. La Tarifa de Permiso (Permit Fee, USD) es la variable dependiente (y), pues se cobra en función de la profundidad alcanzada. A mayor TVD, mayor tarifa, con un crecimiento que se ajusta a una tasa proporcional constante (elasticidad), de ahí el modelo potencial.
# Localizamos las columnas de forma robusta (los nombres pueden llegar
# alterados por R como "True.Vertical.Depth..ft" y "Permit.Fee")
col_tvd <- grep("True.?Vertical.?Depth", names(Datos), value = TRUE)[1]
col_fee <- grep("^Permit.?Fee$", names(Datos), value = TRUE)[1]
if (is.na(col_tvd)) stop("ERROR: No se encontró la columna de Profundidad Vertical Real (TVD).")
if (is.na(col_fee)) stop("ERROR: No se encontró la columna de Tarifa de Permiso (Permit Fee).")
# Selección de variables
datos_raw <- Datos %>%
select(all_of(c(col_tvd, col_fee))) %>%
setNames(c("tvd", "permit_fee")) %>%
mutate(
x_raw = abs(as.numeric(str_replace(as.character(tvd), ",", "."))),
y_raw = abs(as.numeric(str_replace(as.character(permit_fee), ",", ".")))
) %>%
filter(!is.na(x_raw) & !is.na(y_raw) & x_raw > 0 & y_raw > 0) %>%
filter(x_raw <= 27500 & y_raw <= 11000) # rango de dominio válido de cada variable
Se indica el tamaño muestral obtenido y se muestran únicamente las primeras filas de los pares de valores (x, y).
cat("Tamaño muestral: N =", nrow(datos_raw), "pares de valores (TVD, Tarifa de Permiso)\n")
## Tamaño muestral: N = 11045 pares de valores (TVD, Tarifa de Permiso)
datos_raw %>%
select(x_raw, y_raw) %>%
head(10) %>%
gt() %>%
tab_header(title = md("**Tabla N°1: Primeras filas de los pares de valores (x, y)**")) %>%
cols_label(x_raw = "TVD (ft)", y_raw = "Tarifa de Permiso (USD)") %>%
cols_align(align = "center", columns = everything())
| Tabla N°1: Primeras filas de los pares de valores (x, y) | |
| TVD (ft) | Tarifa de Permiso (USD) |
|---|---|
| 1800 | 860 |
| 1445 | 475 |
| 2070 | 725 |
| 6321 | 1725 |
| 1206 | 475 |
| 1590 | 475 |
| 9914 | 380 |
| 9674 | 5230 |
| 1529 | 1050 |
| 4281 | 1225 |
par(mar = c(5, 5, 4, 2))
color_trans <- rgb(0.2, 0.6, 0.86, 0.4)
plot(datos_raw$x_raw, datos_raw$y_raw,
main = "Gráfica N°1: Diagrama de Dispersión de la Tarifa de Permiso\nen función de la Profundidad Vertical Real (TVD)",
xlab = "TVD (ft)",
ylab = "Tarifa de Permiso (USD)",
col = color_trans,
pch = 16,
cex = 0.6,
cex.main = 0.9,
frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
La nube de puntos observada en la Gráfica N°1 presenta una alta variabilidad, con un comportamiento disperso y ruidoso que dificulta identificar visualmente una tendencia clara. Al tratarse de una nube compleja/caótica, se procede con el tratamiento de los datos descrito a continuación.
Para reducir el ruido y resaltar la tendencia, se segmenta x en bins de 500 ft (único x), se promedia y por bin (único y), y se omiten outliers: bins con menos de 3 datos y valores de y fuera del percentil 5%-95%.
# Agrupamiento cada 500 unidades de TVD
datos_model <- datos_raw %>%
mutate(x_bin = round(x_raw / 500) * 500) %>%
group_by(x_bin) %>%
summarise(
y = mean(y_raw, na.rm = TRUE),
conteo = n(),
.groups = "drop"
) %>%
rename(x = x_bin) %>%
filter(conteo >= 3, x > 0) # se excluye x = 0 (ln(0) no está definido)
# Limpieza de outliers
lim_y <- quantile(datos_model$y, probs = c(0.05, 0.95))
datos_model <- datos_model %>%
filter(y >= lim_y[1] & y <= lim_y[2])
x <- datos_model$x
y <- datos_model$y
par(mar = c(5, 5, 4, 2))
plot(x, y,
main = "Gráfica N°2: Dispersión de la Tarifa de Permiso promedio\npor intervalos de TVD (ft)",
xlab = "TVD (ft)",
ylab = "Tarifa de Permiso promedio (USD)",
col = "#3498DB",
pch = 16,
cex = 0.9,
cex.main = 0.9,
frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
Tras el tratamiento de los datos, la nube de puntos muestra una tendencia creciente compatible con una tasa de crecimiento proporcional (elasticidad constante) entre las variables, patrón compatible con un modelo potencial:
\[y = a \cdot x^{b}\]
Para linealizar el modelo se aplica logaritmo natural a ambos lados de la ecuación:
\[\ln(y) = \ln(a) + b \cdot \ln(x)\]
modelo_pot <- lm(log(y) ~ log(x))
ln_a_bin <- coef(modelo_pot)[1]
b_bin <- coef(modelo_pot)[2]
a_bin <- exp(ln_a_bin)
ecuacion <- paste0("y = ",
round(a_bin, 4),
" * x^",
round(b_bin, 4))
cat("La ecuación estimada del modelo es:\n\n", ecuacion)
## La ecuación estimada del modelo es:
##
## y = 1.4542 * x^0.8175
r <- cor(log(x), log(y), use = "complete.obs")
r2 <- summary(modelo_pot)$r.squared
tabla_resumen <- data.frame(
Variable = c("TVD (ft)", "Tarifa de Permiso (USD)"),
Tipo = c("Independiente (x)", "Dependiente (y)"),
R = c("", round(r, 2)),
R2 = c("", round(r2, 2)),
Coeficiente_a = c("", round(a_bin, 4)),
Exponente_b = c("", round(b_bin, 4)),
Ecuación = c("", ecuacion)
)
tabla_resumen %>%
gt() %>%
tab_header(title = md("**Tabla N°2: Resumen del Modelo de Regresión Potencial**")) %>%
tab_source_note(source_note = "Autor: JENNY") %>%
cols_align(align = "center", columns = everything())
| Tabla N°2: Resumen del Modelo de Regresión Potencial | ||||||
| Variable | Tipo | R | R2 | Coeficiente_a | Exponente_b | Ecuación |
|---|---|---|---|---|---|---|
| TVD (ft) | Independiente (x) | |||||
| Tarifa de Permiso (USD) | Dependiente (y) | 0.96 | 0.92 | 1.4542 | 0.8175 | y = 1.4542 * x^0.8175 |
| Autor: JENNY | ||||||
Se presenta el ajuste del modelo sobre los datos reales (agrupados), incluyendo la banda de incertidumbre estadística (Intervalo de Confianza del 95%). El intervalo se calcula sobre el modelo linealizado y se transforma de vuelta a la escala original con la función exponencial.
par(mar = c(5, 5, 4, 2))
plot(x, y,
main = "Gráfica N°3: Modelo Potencial de la Tarifa de Permiso\nen función del TVD (ft)",
xlab = "TVD (ft)",
ylab = "Tarifa de Permiso (USD)",
col = "#3498DB",
pch = 16,
cex = 1.0,
cex.main = 0.9,
frame.plot = FALSE)
grid(nx = NULL, ny = NULL, col = "#D7DBDD", lty = "dotted")
# Secuencia suave
x_seq <- seq(min(x), max(x), length.out = 500)
pred_pot <- predict(modelo_pot,
newdata = data.frame(x = x_seq),
interval = "confidence",
level = 0.95)
# Se transforma de la escala logarítmica a la escala original
pred_pot_original <- exp(pred_pot)
# Intervalo de confianza
polygon(c(x_seq, rev(x_seq)),
c(pred_pot_original[, "lwr"], rev(pred_pot_original[, "upr"])),
col = rgb(0.5, 0.5, 0.5, 0.2),
border = NA)
# Línea ajustada
lines(x_seq, pred_pot_original[, "fit"], col = "#E74C3C", lwd = 3)
legend("topleft",
legend = c("Datos promediados (binning)",
"Modelo Potencial",
"I.C. 95%"),
col = c("#3498DB", "#E74C3C", "gray"),
pch = c(16, NA, 15),
lwd = c(NA, 3, NA),
pt.cex = c(1, NA, 2),
bty = "n")
cat("El coeficiente de correlación es: ", round(r, 2))
## El coeficiente de correlación es: 0.96
cat(paste0("El coeficiente de determinación (R²) es: ", round(r2, 2)))
## El coeficiente de determinación (R²) es: 0.92
El modelo presenta como única condición matemática que x > 0 y y > 0, lo cual se cumple naturalmente en el contexto físico del TVD (ft) y de la Tarifa de Permiso (USD), por lo que no existen restricciones prácticas adicionales dentro del rango analizado.
¿Cuál es la Tarifa de Permiso estimada para un pozo con un TVD de 3000 ft?
tvd_test <- 3000
ln_y_est <- predict(modelo_pot,
newdata = data.frame(x = tvd_test))
y_est <- exp(ln_y_est)
cat("Para un TVD de", tvd_test,
"ft, la Tarifa de Permiso estimada es:",
round(y_est, 4), "USD")
## Para un TVD de 3000 ft, la Tarifa de Permiso estimada es: 1012.127 USD
Entre el TVD (ft) y la Tarifa de Permiso (USD) existe una relación de tipo potencial, con un coeficiente de determinación R² = 0.92, lo que indica un ajuste excelente del modelo.
La ecuación estimada es: y = 1.4542 * x^0.8175.
El modelo presenta como única condición matemática que x > 0 y y > 0, lo cual se cumple naturalmente en el contexto físico del TVD (ft) y de la Tarifa de Permiso (USD), por lo que no existen restricciones prácticas adicionales dentro del rango analizado.