Modelo de Regresión Potencial
Analisis estadísticos de Pozos Regulados del Estado de Nueva York
Eduardo Loayza
2026-07-23
1 Librerías
# Instalación automática de paquetes faltantes (solo si no están ya instalados)
paquetes_necesarios <- c("dplyr", "gt", "kableExtra")
paquetes_faltantes <- paquetes_necesarios[!(paquetes_necesarios %in% installed.packages()[, "Package"])]
if (length(paquetes_faltantes) > 0) {
install.packages(paquetes_faltantes, repos = "https://cloud.r-project.org")
}
library(dplyr)
library(gt)
library(kableExtra)2 Carga de datos
# IMPORTANTE: ajustar esta ruta a la ubicación local del archivo
ruta_carpeta <- "C:/Users/ASUS/Desktop/Estadistica/new_york_exel"
archivo_csv <- "Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv"
if (dir.exists(ruta_carpeta)) {
setwd(ruta_carpeta)
} else {
warning("La carpeta indicada en 'ruta_carpeta' no existe en este equipo. ",
"Se intentará leer 'archivo_csv' tal como está escrito (ruta ",
"relativa al directorio de trabajo actual o ruta absoluta). ",
"Ajusta 'ruta_carpeta' o escribe la ruta completa en 'archivo_csv'.")
}
lineas <- readLines(archivo_csv, encoding = "latin1", warn = FALSE)
Encoding(lineas) <- "latin1"
lineas <- iconv(lineas, from = "latin1", to = "UTF-8")
# Detección de separador (evita depender de un delimitador fijo)
separadores <- c(",", ";", "\t", "|")
mejor_sep <- NULL
mejor_ncol <- 1
for (s in separadores) {
n_campos <- tryCatch(length(strsplit(lineas[1], s, fixed = TRUE)[[1]]),
error = function(e) 1)
if (!is.na(n_campos) && n_campos > mejor_ncol) {
mejor_ncol <- n_campos
mejor_sep <- s
}
}
if (is.null(mejor_sep)) mejor_sep <- ";"
cat("Separador detectado:", ifelse(mejor_sep == "\t", "TAB", mejor_sep),
"| Columnas detectadas:", mejor_ncol, "\n")## Separador detectado: , | Columnas detectadas: 58datos <- read.csv(text = lineas, sep = mejor_sep, header = TRUE,
check.names = TRUE, stringsAsFactors = FALSE)
cat("Número de registros:", nrow(datos), "\n")## Número de registros: 47390cat("Número de variables:", ncol(datos), "\n")## Número de variables: 52Extracto del dataset (55 variables, primeras 5 filas):
datos %>%
head(5) %>%
kbl(caption = "Extracto del dataset: todas las columnas, primeras 5 filas") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
font_size = 11, full_width = FALSE) %>%
scroll_box(width = "100%", height = "420px")| API.Well.Number | County.Code | API.Hole.Number | Sidetrack | Completion | Well.Name | Company.Name | Operator.Number | Well.Type | Map.Symbol | Well.Status | Status.Date | Permit.Application.Date | Permit.Issued.Date | Date.Spudded | Date.of.Total.Depth | Date.Well.Completed | Date.Well.Plugged | Date.Well.Confidentiality.Ends | Confidentiality.Code | Town | Quad | Quad.Section | Producing.Field | Producing.Formation | Financial.Security | Slant | County | Region | State.Lease | Proposed.Depth..ft | Surface.Longitude | Surface.Latitude | Bottom.Hole.Longitude | Bottom.Hole.Latitude | True.Vertical.Depth..ft | Measured.Depth..ft | Kickoff..ft | Drilled.Depth..ft | Elevation..ft | Original.Well.Type | Permit.Fee | Objective.Formation | Depth.Fee | Spacing | Spacing.Acres | Integration | Hearing.Date | Date.Last.Modified | DEC.Database.Link | Location.1 | Georeference |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 31003026700000 | 3 | 2670 | 0 | 0 | Francisco 1 | Van Gilder | 9279 | DW | DP | PA | 03/11/1953 | 04/10/1953 | Pre-1989 Well (N/A) | Amity | Belmont | F | False | Vertical | Allegany | 9 | NA | 0 | -78.01913 | 42.19713 | -78.01913 | 42.19713 | 2006 | 2006 | 0 | 2006 | 1815 | NL | 0 | 0 | NA | 10/12/1995 12:00:00 AM | http://extapps.dec.ny.gov/cfmx/extapps/GasOil/search/wells/index.cfm?api=31003026700000 | (42.19713, -78.01913) | POINT (-78.01913 42.19713) | ||||||||||||
| 31003045990000 | 3 | 4599 | 0 | 0 | Francisco 1 | Christman Raymond L. | 9111 | OW | OW | UN | Pre-1989 Well (N/A) | Amity | Belmont | F | False | Vertical | Allegany | 9 | NA | 0 | -78.02615 | 42.20221 | -78.02615 | 42.20221 | 0 | 0 | 0 | 0 | 1520 | NL | 0 | 0 | NA | 01/24/2003 03:53:37 PM | http://extapps.dec.ny.gov/cfmx/extapps/GasOil/search/wells/index.cfm?api=31003045990000 | (42.20221, -78.02615) | POINT (-78.02615 42.20221) | ||||||||||||||
| 31003048420000 | 3 | 4842 | 0 | 0 | Guyer Devonian #10 | Pennzoil Products Co. | 29 | NL | O | VP | Pre-1989 Well (N/A) | False | Vertical | 9 | NA | NA | NA | NA | NA | NA | 0 | 0 | 0 | 0 | NA | NL | 0 | 0 | NA | 12/27/2013 03:00:05 PM | http://extapps.dec.ny.gov/cfmx/extapps/GasOil/search/wells/index.cfm?api=31003048420000 | ||||||||||||||||||||
| 31003054190000 | 3 | 5419 | 0 | 0 | Regan 2142 | Iroquois Gas Corp. | 16 | GD | GWP | PA | Pre-1989 Well (N/A) | Alma | Wellsville South | D | False | Vertical | Allegany | 9 | NA | 0 | -77.97909 | 42.07020 | -77.97909 | 42.07020 | 0 | 0 | 0 | 0 | 2100 | NL | 0 | 0 | NA | 02/28/1995 12:00:00 AM | http://extapps.dec.ny.gov/cfmx/extapps/GasOil/search/wells/index.cfm?api=31003054190000 | (42.0702, -77.97909) | POINT (-77.97909 42.0702) | ||||||||||||||
| 31003069290000 | 3 | 6929 | 0 | 0 | 2667 | Bradley Producing Company | 9673 | NL | O | VP | Pre-1989 Well (N/A) | False | Vertical | Statewide | 9 | NA | NA | NA | NA | NA | NA | 0 | 0 | 0 | 0 | NA | NL | 0 | 0 | NA | 11/28/2018 10:25:33 AM | http://extapps.dec.ny.gov/cfmx/extapps/GasOil/search/wells/index.cfm?api=31003069290000 |
3 Selección de variables
La tarifa por profundidad (Depth Fee)
se calcula por reglamento a partir de la profundidad
propuesta (Proposed Depth), según una tabla de
tarifas escalonada. Hay una relación causa (X) → efecto (Y) clara:
- X:
Proposed Depth, ft. - Y:
Depth Fee(USD).
Las columnas se localizan por patrón con grepl() en vez
de por nombre fijo, ya que check.names = TRUE transforma
los encabezados originales (espacios y comas por puntos).
col_proposed <- names(datos)[grepl("proposed.*depth", names(datos), ignore.case = TRUE)]
col_fee <- names(datos)[grepl("^depth.*fee", names(datos), ignore.case = TRUE)]
if (length(col_proposed) == 0 || length(col_fee) == 0) {
columnas_txt <- paste(names(datos), collapse = " | ")
stop("ERROR: No se encontraron las columnas esperadas ",
"(se buscaron los patrones 'proposed...depth' y 'depth...fee').\n",
"COLUMNAS ENCONTRADAS EN EL ARCHIVO:\n", columnas_txt)
}
col_proposed <- col_proposed[1]
col_fee <- col_fee[1]
cat("Columna Proposed Depth identificada:", col_proposed, "\n")## Columna Proposed Depth identificada: Proposed.Depth..ftcat("Columna Depth Fee identificada:", col_fee, "\n")## Columna Depth Fee identificada: Depth.Feex_raw <- suppressWarnings(as.numeric(as.character(datos[[col_proposed]])))
y_raw <- suppressWarnings(as.numeric(as.character(datos[[col_fee]])))
cat("Registros con Proposed Depth:", sum(!is.na(x_raw) & x_raw > 0), "\n")## Registros con Proposed Depth: 16507cat("Registros con Depth Fee:", sum(!is.na(y_raw) & y_raw > 0), "\n")## Registros con Depth Fee: 11696Solo se conservan los pares donde ambas variables son mayores que 0; un valor de 0 o vacío indica que el permiso no llegó a esa etapa administrativa.
4 Tabla de pares de valores
pares_iniciales <- data.frame(x = x_raw, y = y_raw) %>%
filter(!is.na(x), !is.na(y), x > 0, y > 0)
cat("Tamaño muestral (pares completos X > 0 y Y > 0):", nrow(pares_iniciales), "\n")## Tamaño muestral (pares completos X > 0 y Y > 0): 11175Vista previa (primeras filas):
tabla_ini_gt <- pares_iniciales %>%
head(10) %>%
rename(`Proposed Depth, ft (X)` = x,
`Depth Fee, USD (Y)` = y)
tabla_ini_gt %>%
gt() %>%
tab_header(
title = md("**Tabla de Pares de Valores (vista previa)**"),
subtitle = md(paste0("Tamaño muestral n = ", nrow(pares_iniciales)))
) %>%
tab_source_note(source_note = "Autor: Grupo 1") %>%
cols_align(align = "center", columns = everything()) %>%
tab_style(
style = list(cell_fill(color = col_encabezado), cell_text(color = "white", weight = "bold")),
locations = cells_title()
) %>%
tab_style(
style = list(cell_fill(color = col_encabezado), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()
) %>%
tab_style(
style = list(cell_fill(color = col_fila_alt)),
locations = cells_body(rows = seq(1, nrow(tabla_ini_gt), 2))
) %>%
opt_table_outline(style = "solid", width = px(3), color = col_borde) %>%
tab_options(
table.border.top.color = col_borde,
table.border.bottom.color = col_borde,
table.border.top.style = "solid",
table.border.bottom.style = "solid",
column_labels.border.top.color = col_borde,
column_labels.border.bottom.color = col_borde,
column_labels.border.bottom.width = px(2),
heading.border.bottom.color = col_borde,
heading.border.bottom.width = px(2),
table_body.hlines.color = "#CBD5DE",
table_body.border.bottom.color = col_borde,
table.font.size = px(14),
data_row.padding = px(6),
table_body.border.top.style = "solid",
column_labels.background.color = col_encabezado
)| Tabla de Pares de Valores (vista previa) | |
| Tamaño muestral n = 11175 | |
| Proposed Depth, ft (X) | Depth Fee, USD (Y) |
|---|---|
| 900 | 250 |
| 4200 | 1710 |
| 1555 | 500 |
| 1600 | 500 |
| 1783 | 760 |
| 1600 | 760 |
| 1750 | 760 |
| 1530 | 500 |
| 1250 | 375 |
| 2350 | 625 |
| Autor: Grupo 1 | |
5 Gráfica de dispersión
plot(pares_iniciales$x, pares_iniciales$y,
pch = 20,
col = rgb(0.04, 0.18, 0.29, 0.15),
xlab = "Proposed Depth, ft (X)",
ylab = "Depth Fee, USD (Y)",
main = "Gráfica original: Proposed Depth en comparación con Depth Fee")
6 Conjetura
La nube de puntos original es compleja/caótica: un mismo X puede tener varios Y (ajustes de tarifa por año) y hay tarifas desproporcionadas (probables errores de digitación). Por eso se trata el dato antes de modelar.
6.1 Tratamiento de los datos
Paso A — Omisión de outliers. Se calcula la razón
Depth Fee / Proposed Depth; si queda fuera de 0.05–1.0
USD/ft se descarta por ser casi con certeza un error de digitación.
pares_tratados <- pares_iniciales %>% mutate(razon = y / x)
n_antes <- nrow(pares_tratados)
pares_tratados <- pares_tratados %>% filter(razon >= 0.05, razon <= 1.0) %>% select(-razon)
n_despues <- nrow(pares_tratados)
cat("Pares antes del filtro de coherencia:", n_antes, "\n")## Pares antes del filtro de coherencia: 11175cat("Pares tras el filtro de coherencia:", n_despues, "\n")## Pares tras el filtro de coherencia: 11134cat("Descartados por razón no coherente (outliers):", n_antes - n_despues, "\n")## Descartados por razón no coherente (outliers): 41Paso B — Único X, único Y. Cuando una misma X tiene varios valores de Y, se calcula la media para obtener un único par por valor exacto de X (no es binning, no se crean rangos artificiales).
pares <- pares_tratados %>%
group_by(x) %>%
summarise(y = mean(y, na.rm = TRUE), .groups = "drop") %>%
arrange(x)
cat("Pares únicos tras el tratamiento (único X, único Y):", nrow(pares), "\n")## Pares únicos tras el tratamiento (único X, único Y): 2319Tras el tratamiento, la relación es clara y no requiere segmentar por partes.
6.2 Nueva gráfica de dispersión
plot(pares$x, pares$y,
pch = 20,
col = rgb(0.04, 0.18, 0.29, 0.35),
xlab = "Proposed Depth, ft (X)",
ylab = "Depth Fee, USD (Y)",
main = "Gráfica tras el tratamiento de datos")
6.3 Nueva conjetura
La nube muestra una tendencia creciente y cóncava: la tarifa aumenta con la profundidad pero a un ritmo cada vez menor, típico de un modelo potencial:
\[y = a \cdot x^{b}\]
Al no ser lineal en sus parámetros, se linealiza aplicando logaritmo natural a ambas variables:
\[\ln(y) = \ln(a) + b \cdot \ln(x)\]
Se ajusta lm() sobre \((\ln x,
\ln y)\) y luego se recupera \(a\) con exp().
pares$log_x <- log(pares$x)
pares$log_y <- log(pares$y)
plot(pares$log_x, pares$log_y,
pch = 20,
col = rgb(0.04, 0.18, 0.29, 0.35),
xlab = "ln(Proposed Depth, ft) — eje X reducido",
ylab = "ln(Depth Fee, USD) — eje Y reducido",
main = "Linealización: ln(Y) en comparación con ln(X)")
7 Cálculo de parámetros
m_pot <- lm(log_y ~ log_x, data = pares)
sum_reg <- summary(m_pot)coefs <- coef(m_pot)
intercepto <- coefs[1]
pendiente <- coefs[2]
a <- exp(intercepto)
b <- pendiente
cat("Intercepto en espacio log [ln(a)] :", round(intercepto, 6), "\n")## Intercepto en espacio log [ln(a)] : 0.080781cat("Pendiente (b) :", round(b, 6), "\n")## Pendiente (b) : 0.858846cat("Parámetro a = exp(intercepto) :", round(a, 4), "\n")## Parámetro a = exp(intercepto) : 1.0841cat("\nEcuación del modelo potencial:\n")##
## Ecuación del modelo potencial:cat("y =", round(a, 4), "* x^(", round(b, 6), ")\n")## y = 1.0841 * x^( 0.858846 )Como \(b \approx 0.86 < 1\), el crecimiento es sublineal: la tarifa crece más lentamente que la profundidad.
8 Realidad y modelo
x_grid <- seq(min(pares$x), max(pares$x), length.out = 400)
log_x_grid <- log(x_grid)
log_y_grid <- predict(m_pot, newdata = data.frame(log_x = log_x_grid))
y_grid <- exp(log_y_grid)
plot(pares$x, pares$y,
pch = 20,
col = rgb(0.04, 0.18, 0.29, 0.15),
xlab = "Proposed Depth, ft (X)",
ylab = "Depth Fee, USD (Y)",
main = "Superposición: Modelo Potencial y Datos Reales")
lines(x_grid, y_grid, col = col_curva, lwd = 3)
legend("topleft",
legend = c("Datos reales", "Modelo potencial"),
col = c(col_puntos, col_curva),
pch = c(20, NA),
lty = c(NA, 1),
lwd = c(NA, 3),
bty = "n")
9 Test
Como el modelo se ajustó linealizado, el test de Pearson se calcula entre \(\ln(X)\) y \(\ln(Y)\), donde la relación es realmente lineal.
r <- cor(pares$log_x, pares$log_y)
r2 <- sum_reg$r.squared * 100
cat("Correlación de Pearson (r), entre ln(X) y ln(Y):", round(r, 4), "\n")## Correlación de Pearson (r), entre ln(X) y ln(Y): 0.88cat("Coeficiente de determinación (R²%):", round(r2, 2), "%\n")## Coeficiente de determinación (R²%): 77.44 %Con 0.88 de correlación de Pearson, el modelo muestra una asociación fuerte, coherente con una tarifa que depende de la profundidad por reglamento pero no de forma perfectamente determinística.
10 Restricciones
x_min <- round(min(pares$x), 2)
x_max <- round(max(pares$x), 2)
cat("Dominio válido de X (Proposed Depth, ft):", x_min, "a", x_max, "\n")## Dominio válido de X (Proposed Depth, ft): 200 a 13450- Válido solo entre 200 ft y 13450 ft; no extrapolar fuera de este rango.
- Requiere \(x > 0\) (profundidad propuesta estrictamente positiva).
- Describe una tendencia promedio, no el monto exacto de un permiso individual, que puede variar por año de vigencia u otras condiciones administrativas.
11 Estimación
x_estimar1 <- x_min
y_estimada1 <- a * x_estimar1^b
cat("Estimación para X =", x_estimar1, "ft:\n")## Estimación para X = 200 ft:cat("Depth Fee estimado: $", round(y_estimada1, 2), "\n\n")## Depth Fee estimado: $ 102.64x_estimar2 <- x_max
y_estimada2 <- a * x_estimar2^b
cat("Estimación para X =", x_estimar2, "ft:\n")## Estimación para X = 13450 ft:cat("Depth Fee estimado: $", round(y_estimada2, 2), "\n\n")## Depth Fee estimado: $ 3810.83# Ejemplo puntual de referencia
x_estimar3 <- 1500
y_estimada3 <- a * x_estimar3^b
cat("Estimación para X =", x_estimar3, "ft:\n")## Estimación para X = 1500 ft:cat("Depth Fee estimado: $", round(y_estimada3, 2), "\n")## Depth Fee estimado: $ 579.24Tabla resumen del modelo:
Ecuacion <- paste0("y = ", round(a, 4), " * x^(", round(b, 6), ")")
Tabla_resumen <- data.frame(
`Variable Independiente` = "Proposed Depth, ft",
`Variable Dependiente` = "Depth Fee, USD",
`Test Pearson (ln X y ln Y)` = round(r, 4),
`Coeficiente de determinación` = round(r2, 2),
`Ecuación del modelo` = Ecuacion,
check.names = FALSE
)
Tabla_resumen %>%
gt() %>%
tab_header(
title = md("**Tabla N°1**"),
subtitle = md("**Resumen del modelo de regresión potencial**")
) %>%
tab_source_note(source_note = md("Autor: Grupo 1")) %>%
cols_align(align = "center", columns = everything()) %>%
tab_style(
style = list(cell_fill(color = col_encabezado), cell_text(color = "white", weight = "bold")),
locations = cells_title()
) %>%
tab_style(
style = list(cell_fill(color = col_encabezado), cell_text(color = "white", weight = "bold")),
locations = cells_column_labels()
) %>%
tab_style(
style = list(cell_fill(color = col_fila_alt)),
locations = cells_body(rows = 1)
) %>%
opt_table_outline(style = "solid", width = px(3), color = col_borde) %>%
tab_options(
table.border.top.color = col_borde,
table.border.bottom.color = col_borde,
table.border.top.style = "solid",
table.border.bottom.style = "solid",
column_labels.border.top.color = col_borde,
column_labels.border.bottom.color = col_borde,
column_labels.border.bottom.width = px(2),
heading.border.bottom.color = col_borde,
heading.border.bottom.width = px(2),
table_body.hlines.color = "#CBD5DE",
table_body.border.bottom.color = col_borde,
table.font.size = px(14),
data_row.padding = px(6)
)| Tabla N°1 | ||||
| Resumen del modelo de regresión potencial | ||||
| Variable Independiente | Variable Dependiente | Test Pearson (ln X y ln Y) | Coeficiente de determinación | Ecuación del modelo |
|---|---|---|---|---|
| Proposed Depth, ft | Depth Fee, USD | 0.88 | 77.44 | y = 1.0841 * x^(0.858846) |
| Autor: Grupo 1 | ||||
12 Conclusión
La relación entre Proposed Depth (X) y
Depth Fee (Y) es potencial y
sublineal:
\[y = 1.0841 \cdot x^{0.858846}\]
Con Pearson 0.88 y R² de 77.44%, el modelo es válido solo entre 200 y 13450 ft.