Procedimiento de Regresión Clásica — Regresión Lineal Múltiple (Oil & Gas Leases: Producción Acumulada, Producción Promedio y Años Activos)

Usar Run All (Code → Run Region → Run All), no Knit, porque el chunk de carga usa file.choose().

Dataset: oil_and_gas_leases_data.csv (47,757 pozos). De las 13 variables numéricas disponibles, se evaluaron todas las combinaciones posibles de 2 variables independientes y se seleccionaron las 3 que dieron el mejor resultado en el test de Pearson (ver sección 9):

\[y = b_0 + b_1 \cdot x_1 + b_2 \cdot x_2\]

1. Librerías

library(dplyr)
library(knitr)
library(kableExtra)
library(DT)
library(plotly)

2. Carga de datos

ruta <- file.choose()
datos <- read.csv(ruta, stringsAsFactors = FALSE)

datos <- datos %>%
  mutate(
    CUMULATIVE_PRODUCTION = as.numeric(CUMULATIVE_PRODUCTION),
    AVG_PRODUCTION        = as.numeric(AVG_PRODUCTION),
    YEARS_ACTIVE          = as.numeric(YEARS_ACTIVE)
  ) %>%
  select(CUMULATIVE_PRODUCTION, AVG_PRODUCTION, YEARS_ACTIVE) %>%
  na.omit() %>%
  filter(CUMULATIVE_PRODUCTION > 0, AVG_PRODUCTION > 0, YEARS_ACTIVE > 0)

cat("Tamaño muestral original (n):", nrow(datos))

## Tamaño muestral original (n): 47757

3. Selección de variables (causa – efecto)

Efecto (Y): CUMULATIVE_PRODUCTION — Causas (X1, X2): AVG_PRODUCTION, YEARS_ACTIVE

Justificación de la relación causa-efecto:

• Mecanismo físico directo: la producción acumulada de un pozo es, por construcción, el resultado de acumular producción a lo largo del tiempo. Si un pozo produce en promedio AVG_PRODUCTION por año y ha estado activo YEARS_ACTIVE años, su producción acumulada será aproximadamente el producto de ambas: CUMULATIVE_PRODUCTION ≈ AVG_PRODUCTION × YEARS_ACTIVE. Esto hace que la relación causa-efecto sea, en este caso, prácticamente una identidad física.
• Precedencia temporal: tanto la tasa de producción promedio como los años activos son “insumos” del proceso — se generan mientras el pozo opera — y la producción acumulada es la consecuencia que resulta de ir sumando esa producción año tras año.
• Direccionalidad: no tendría sentido decir que la producción acumulada “causa” el promedio de producción o los años activos; es al revés — el ritmo de producción y el tiempo de vida del pozo son los que determinan cuánto termina acumulando en total.

Por lo tanto: AVG_PRODUCTION (X1) y YEARS_ACTIVE (X2) son las variables independientes (causas), y CUMULATIVE_PRODUCTION (Y) es la variable dependiente (efecto).

4. Tabla de pares de valores

4.1 Tabla condensada por YEARS_ACTIVE (todos los valores originales, celdas contraíbles)

condensado <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    n_valores = n(),
    valores_avg = paste(AVG_PRODUCTION, collapse = ", "),
    valores_cum = paste(CUMULATIVE_PRODUCTION, collapse = ", "),
    .groups   = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

condensado_html <- condensado %>%
  mutate(
    AVG_PRODUCTION = paste0(
      "<details><summary>", n_valores, " valor(es)</summary>", valores_avg, "</details>"
    ),
    CUMULATIVE_PRODUCTION = paste0(
      "<details><summary>", n_valores, " valor(es)</summary>", valores_cum, "</details>"
    )
  ) %>%
  select(YEARS_ACTIVE, AVG_PRODUCTION, CUMULATIVE_PRODUCTION)

datatable(condensado_html, escape = FALSE, rownames = FALSE,
          colnames = c("Años Activos (X2)", "Producción Promedio (valores, X1)", "Producción Acumulada (valores, Y)"),
          class = "stripe hover compact",
          options = list(pageLength = 10))

4.2 Tabla de tríadas (X1, X2, Y) usando la mediana por YEARS_ACTIVE — todos los valores

tabla_mediana <- datos %>%
  group_by(YEARS_ACTIVE) %>%
  summarise(
    AVG_PRODUCTION        = median(AVG_PRODUCTION),
    CUMULATIVE_PRODUCTION = median(CUMULATIVE_PRODUCTION),
    .groups = "drop"
  ) %>%
  arrange(YEARS_ACTIVE)

cat("Tamaño muestral original:", nrow(datos), "registros\n")

## Tamaño muestral original: 47757 registros

cat("Tamaño muestral de la tabla de tríadas (una por cada YEARS_ACTIVE):", nrow(tabla_mediana), "pares\n")

## Tamaño muestral de la tabla de tríadas (una por cada YEARS_ACTIVE): 89 pares

datatable(tabla_mediana, rownames = FALSE,
          colnames = c("Años Activos (X2)", "Mediana Producción Promedio (X1)", "Mediana Producción Acumulada (Y)"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = c("AVG_PRODUCTION", "CUMULATIVE_PRODUCTION"), digits = 2)

write.csv(condensado[, c("YEARS_ACTIVE", "valores_avg", "valores_cum")], "tabla_condensada_por_anios_oilgas.csv", row.names = FALSE)
write.csv(tabla_mediana, "tabla_mediana_por_anios_oilgas.csv", row.names = FALSE)

5. Gráfica de dispersión (3D)

plot_ly(
  data = tabla_mediana,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers",
  marker = list(size = 5, color = "#1b2a4a", opacity = 0.9)
) %>%
  layout(
    title = "Nube de puntos 3D — tabla de tríadas (mediana)",
    scene = list(
      xaxis = list(title = "Producción Promedio (X1)"),
      yaxis = list(title = "Años Activos (X2)"),
      zaxis = list(title = "Producción Acumulada (Y)")
    )
  )

6. Conjetura

La nube de puntos 3D de la sección 5 ya muestra una tendencia relativamente clara y creciente en ambos ejes, aunque persisten algunos puntos alejados (outliers) que conviene tratar antes de ajustar el modelo, según lo indicado en 6.1.

6.1 Tratamiento de los datos

Se aplican, en orden, los siguientes pasos:

1. Único x, único y: ya se resolvió en el paso 4 — cada valor de YEARS_ACTIVE (X2) aparece una sola vez en tabla_mediana, emparejado con la mediana de AVG_PRODUCTION (X1) y de CUMULATIVE_PRODUCTION (Y) de ese grupo.
2. Omisión de outliers: se usa el criterio de rango intercuartílico (IQR) sobre CUMULATIVE_PRODUCTION y sobre AVG_PRODUCTION: se descartan los pares cuyo valor esté fuera de \([Q1 - 1.5 \cdot IQR,\ Q3 + 1.5 \cdot IQR]\) en cualquiera de las dos variables.
3. Segmentar por partes: no fue necesario — tras quitar los outliers la nube de puntos restante mantiene una tendencia uniforme en todo el rango de X1 y X2.

Q1_y <- quantile(tabla_mediana$CUMULATIVE_PRODUCTION, 0.25)
Q3_y <- quantile(tabla_mediana$CUMULATIVE_PRODUCTION, 0.75)
IQR_y <- Q3_y - Q1_y

Q1_x1 <- quantile(tabla_mediana$AVG_PRODUCTION, 0.25)
Q3_x1 <- quantile(tabla_mediana$AVG_PRODUCTION, 0.75)
IQR_x1 <- Q3_x1 - Q1_x1

datos_limpios <- tabla_mediana %>%
  filter(
    CUMULATIVE_PRODUCTION >= (Q1_y - 1.5 * IQR_y),
    CUMULATIVE_PRODUCTION <= (Q3_y + 1.5 * IQR_y),
    AVG_PRODUCTION >= (Q1_x1 - 1.5 * IQR_x1),
    AVG_PRODUCTION <= (Q3_x1 + 1.5 * IQR_x1)
  )

cat("Pares antes del tratamiento de outliers:", nrow(tabla_mediana), "\n")

## Pares antes del tratamiento de outliers: 89

cat("Pares después del tratamiento de outliers:", nrow(datos_limpios), "\n")

## Pares después del tratamiento de outliers: 85

cat("Pares eliminados:", nrow(tabla_mediana) - nrow(datos_limpios), "\n")

## Pares eliminados: 4

datatable(datos_limpios, rownames = FALSE,
          colnames = c("Años Activos (X2)", "Producción Promedio (X1)", "Producción Acumulada (Y)"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = c("AVG_PRODUCTION", "CUMULATIVE_PRODUCTION"), digits = 2)

6.2 Nueva gráfica de dispersión (3D)

plot_ly(
  data = datos_limpios,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers",
  marker = list(size = 5, color = "#1b2a4a", opacity = 0.9)
) %>%
  layout(
    title = "Nube de puntos 3D — datos sin outliers",
    scene = list(
      xaxis = list(title = "Producción Promedio (X1)"),
      yaxis = list(title = "Años Activos (X2)"),
      zaxis = list(title = "Producción Acumulada (Y)")
    )
  )

6.3 Nueva conjetura

Tras quitar los valores extremos, la nube de puntos 3D muestra una tendencia claramente plana y creciente: la producción acumulada crece de forma proporcional tanto con la producción promedio como con los años activos. Se conjetura un modelo de regresión lineal múltiple, que geométricamente corresponde a un plano en el espacio (X1, X2, Y):

\[y = b_0 + b_1 \cdot x_1 + b_2 \cdot x_2\]

7. Cálculo de parámetros

modelo_lineal <- lm(CUMULATIVE_PRODUCTION ~ AVG_PRODUCTION + YEARS_ACTIVE,
                     data = datos_limpios)
summary(modelo_lineal)

## 
## Call:
## lm(formula = CUMULATIVE_PRODUCTION ~ AVG_PRODUCTION + YEARS_ACTIVE, 
##     data = datos_limpios)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -88907 -21943  -9819  20320 115119 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.479e+05  1.035e+04  -14.30   <2e-16 ***
## AVG_PRODUCTION  4.931e+01  3.651e+00   13.51   <2e-16 ***
## YEARS_ACTIVE    3.433e+03  1.921e+02   17.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35990 on 82 degrees of freedom
## Multiple R-squared:  0.9378, Adjusted R-squared:  0.9363 
## F-statistic: 618.2 on 2 and 82 DF,  p-value: < 2.2e-16

b0 <- unname(coef(modelo_lineal)[1])
b1 <- unname(coef(modelo_lineal)[2])
b2 <- unname(coef(modelo_lineal)[3])

kable(data.frame(Parámetro = c("b0 (intercepto)", "b1 (pendiente X1)", "b2 (pendiente X2)"),
                  Valor = round(c(b0, b1, b2), 4)),
      align = "lr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Parámetro	Valor
b0 (intercepto)	-147949.4176
b1 (pendiente X1)	49.3065
b2 (pendiente X2)	3432.5675

8. Realidad y modelo (plano de regresión en 3D)

x1_seq <- seq(min(datos_limpios$AVG_PRODUCTION),
              max(datos_limpios$AVG_PRODUCTION), length.out = 30)
x2_seq <- seq(min(datos_limpios$YEARS_ACTIVE),
              max(datos_limpios$YEARS_ACTIVE), length.out = 30)

# Matriz z: filas = x2_seq (eje Y del gráfico), columnas = x1_seq (eje X del gráfico)
z_plano <- outer(x2_seq, x1_seq, function(x2, x1) b0 + b1 * x1 + b2 * x2)

plot_ly() %>%
  add_markers(
    data = datos_limpios,
    x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
    marker = list(size = 5, color = "#5b6b8c", opacity = 0.9),
    name = "Datos reales (sin outliers)"
  ) %>%
  add_surface(
    x = x1_seq, y = x2_seq, z = z_plano,
    opacity = 0.6,
    colorscale = list(c(0, 1), c("#c9ced6", "#1b2a4a")),
    showscale = FALSE,
    name = "Plano del modelo lineal"
  ) %>%
  layout(
    title = "Realidad (puntos) vs. modelo lineal múltiple (plano)",
    scene = list(
      xaxis = list(title = "Producción Promedio (X1)"),
      yaxis = list(title = "Años Activos (X2)"),
      zaxis = list(title = "Producción Acumulada (Y)")
    )
  )

En este gráfico 3D se observan los puntos reales (datos sin outliers) junto con el plano estimado por el modelo lineal (\(y = b_0 + b_1 x_1 + b_2 x_2\)). Los puntos quedan muy cerca del plano, lo que anticipa un buen ajuste (confirmado en la sección 9).

9. Test

r_x1 <- cor(datos_limpios$AVG_PRODUCTION, datos_limpios$CUMULATIVE_PRODUCTION)
r_x2 <- cor(datos_limpios$YEARS_ACTIVE, datos_limpios$CUMULATIVE_PRODUCTION)
r2         <- summary(modelo_lineal)$r.squared
r2_ajust   <- summary(modelo_lineal)$adj.r.squared
p_modelo   <- pf(summary(modelo_lineal)$fstatistic[1],
                  summary(modelo_lineal)$fstatistic[2],
                  summary(modelo_lineal)$fstatistic[3], lower.tail = FALSE)

resultado_x1 <- ifelse(abs(r_x1) > 0.7, "Aceptado", "Rechazado")
resultado_x2 <- ifelse(abs(r_x2) > 0.7, "Aceptado", "Rechazado")
resultado_modelo <- ifelse(p_modelo < 0.05, "Significativo", "No significativo")

tabla_test <- data.frame(
  Test = c("Correlación de Pearson: X1 (AVG_PRODUCTION) vs Y",
           "Correlación de Pearson: X2 (YEARS_ACTIVE) vs Y",
           "Significancia global del modelo (F-test)"),
  Valor = c(round(r_x1, 4), round(r_x2, 4), format(p_modelo, scientific = TRUE, digits = 4)),
  Criterio = c("|r| > 0.7", "|r| > 0.7", "p-valor < 0.05"),
  Resultado = c(resultado_x1, resultado_x2, resultado_modelo)
)

kable(tabla_test, align = "lrrc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  column_spec(4, bold = TRUE,
              color = ifelse(tabla_test$Resultado %in% c("Aceptado", "Significativo"), "#1b2a4a", "#6b6b6b"))

Test	Valor	Criterio	Resultado
Correlación de Pearson: X1 (AVG_PRODUCTION) vs Y	0.8341	|r| > 0.7	Aceptado
Correlación de Pearson: X2 (YEARS_ACTIVE) vs Y	0.8941	|r| > 0.7	Aceptado
Significancia global del modelo (F-test)	3.512e-50	p-valor < 0.05	Significativo

cat("R² del modelo:", round(r2, 4), "\n")

## R² del modelo: 0.9378

cat("R² ajustado del modelo:", round(r2_ajust, 4), "\n")

## R² ajustado del modelo: 0.9363

Tanto AVG_PRODUCTION como YEARS_ACTIVE superan el umbral de |r| > 0.7 frente a CUMULATIVE_PRODUCTION, y el modelo en conjunto explica más del 90% de la variabilidad (R² > 0.9): estas 3 variables sí pasan el test de Pearson.

10. Restricciones

• El modelo es válido dentro del rango de datos con el que fue ajustado (después de la limpieza de outliers):
  • AVG_PRODUCTION (X1): entre 608.56 y 6605.7.
  • YEARS_ACTIVE (X2): entre 1 y 89.
• Al ser un modelo lineal, extrapolar fuera de estos rangos puede producir estimaciones poco confiables (por ejemplo, para pozos con tasas de producción extremadamente altas o muy pocos años activos).
• El modelo se ajustó sobre datos agregados por mediana (una tríada por cada valor de YEARS_ACTIVE) y sin outliers; no debe usarse para predecir pozos individuales atípicos, sino comportamientos típicos por año de actividad.
• Aunque la relación tiene una base física casi determinística (CUMULATIVE_PRODUCTION ≈ AVG_PRODUCTION × YEARS_ACTIVE), el modelo ajustado es lineal aditivo, no multiplicativo, por lo que es una aproximación válida dentro del rango de datos observado.

11. Estimación

nuevos_x1 <- c(500, 1000, 5000, 10000, 50000, 100000)
nuevos_x2 <- c(5, 10, 20, 30, 40, 50)

estimaciones <- expand.grid(AVG_PRODUCTION = nuevos_x1, YEARS_ACTIVE = nuevos_x2)
estimaciones$Produccion_Acumulada_Estimada <- predict(modelo_lineal, newdata = estimaciones)
estimaciones <- estimaciones %>% arrange(AVG_PRODUCTION, YEARS_ACTIVE)

datatable(estimaciones, rownames = FALSE,
          colnames = c("Producción Promedio (X1)", "Años Activos (X2)", "Producción Acumulada Estimada"),
          class = "stripe hover compact",
          options = list(pageLength = 10)) %>%
  formatRound(columns = "Produccion_Acumulada_Estimada", digits = 2)

12. Conclusión

De las variables disponibles en el dataset de leases de petróleo y gas, la combinación AVG_PRODUCTION, YEARS_ACTIVE y CUMULATIVE_PRODUCTION fue la que mejor pasó el test de Pearson. Tras agrupar por mediana y quitar outliers, la relación se ajusta a un modelo de regresión lineal múltiple, representado geométricamente por un plano en el espacio 3D:

\[y = -147949.4 + 49.3065 \cdot x_1 + 3432.568 \cdot x_2\]

siendo y la producción acumulada, x1 la producción promedio y x2 los años activos. El modelo es válido dentro de los rangos de X1 y X2 indicados en la sección 10.

Por cada unidad adicional de producción promedio, la producción acumulada cambia en 49.3065 unidades (manteniendo constantes los años activos); por cada año activo adicional, cambia en 3432.568 unidades (manteniendo constante la producción promedio).

Cuando la producción promedio es de 2,930.19 y los años activos son 45, se espera una producción acumulada de 150,993.5.

La producción acumulada está explicada en un 93.8% conjuntamente por la producción promedio y los años activos (R² del modelo), y en un 6.2% por otros factores no incluidos en el modelo.

De forma individual, el test de Pearson dio r = 0.834 (Aceptado) para la producción promedio, y r = 0.894 (Aceptado) para los años activos — ambas variables pasan el test (|r| > 0.7). La prueba F global del modelo resultó Significativo (p = 3.51e-50).