Regresion Lineal Multiple - Oil & Gas Leases

1. Librerias

suppressMessages(suppressWarnings({
  library(dplyr)
  library(knitr)
  library(kableExtra)
  library(DT)
  library(plotly)
}))

2. Carga de datos

# usar Run All, no Knit
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)

3. Seleccion de variables

Y = CUMULATIVE_PRODUCTION (efecto) X1 = AVG_PRODUCTION, X2 = YEARS_ACTIVE (causas)

La produccion acumulada de un pozo depende de cuanto produce en promedio por año y de cuantos años ha estado activo, por eso X1 y X2 son las causas y Y es el efecto.

4. Tabla de pares de valores

4.1 Tripleta de valores originales

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)

nrow(datos)

## [1] 47757

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))

plot_ly(
  data = datos,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers",
  marker = list(size = 2, opacity = 0.4)
) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

Estrategia

1. Ordenar valores de X1, X2 e Y (un X para cada Y)
2. Llenar espacios donde no haya valores
3. Donde se encuentren valores diversos usar indicadores de posición (mediana)
4. Presentar nueva tabla con la mediana de los valores

4.2 Tabla de tripletas (X1, X2, Y) usando la mediana por YEARS_ACTIVE

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

nrow(tabla_mediana)

## [1] 89

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. Grafica de dispersion

plot_ly(
  data = tabla_mediana,
  x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION,
  type = "scatter3d", mode = "markers"
) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

6. Conjetura

Los puntos se ven mas o menos planos y creciendo, no es una nube caotica, entonces se prueba con regresion lineal multiple (un plano):

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

7. Calculo de parametros

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

## 
## Call:
## lm(formula = CUMULATIVE_PRODUCTION ~ AVG_PRODUCTION + YEARS_ACTIVE, 
##     data = tabla_mediana)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -111120  -28177  -13323   29092  169373 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.459e+05  1.145e+04  -12.73   <2e-16 ***
## AVG_PRODUCTION  4.229e+01  2.978e+00   14.20   <2e-16 ***
## YEARS_ACTIVE    3.921e+03  2.093e+02   18.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45170 on 86 degrees of freedom
## Multiple R-squared:  0.9209, Adjusted R-squared:  0.9191 
## F-statistic: 500.9 on 2 and 86 DF,  p-value: < 2.2e-16

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

c(b0 = b0, b1 = b1, b2 = b2)

##            b0            b1            b2 
## -145857.37971      42.29425    3920.59028

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)	-145857.3797
b1 (pendiente X1)	42.2943
b2 (pendiente X2)	3920.5903

8. Realidad y modelo

x1_seq <- seq(min(tabla_mediana$AVG_PRODUCTION), max(tabla_mediana$AVG_PRODUCTION), length.out = 30)
x2_seq <- seq(min(tabla_mediana$YEARS_ACTIVE), max(tabla_mediana$YEARS_ACTIVE), length.out = 30)
z_plano <- outer(x2_seq, x1_seq, function(x2, x1) b0 + b1 * x1 + b2 * x2)

plot_ly() %>%
  add_markers(data = tabla_mediana, x = ~AVG_PRODUCTION, y = ~YEARS_ACTIVE, z = ~CUMULATIVE_PRODUCTION) %>%
  add_surface(x = x1_seq, y = x2_seq, z = z_plano, opacity = 0.6, showscale = FALSE) %>%
  layout(scene = list(
    xaxis = list(title = "X1 AVG_PRODUCTION"),
    yaxis = list(title = "X2 YEARS_ACTIVE"),
    zaxis = list(title = "Y CUMULATIVE_PRODUCTION")
  ))

Los puntos quedan cerca del plano, buen ajuste.

9. Test

r_x1 <- cor(tabla_mediana$AVG_PRODUCTION, tabla_mediana$CUMULATIVE_PRODUCTION)
r_x2 <- cor(tabla_mediana$YEARS_ACTIVE, tabla_mediana$CUMULATIVE_PRODUCTION)
r2 <- summary(modelo_lineal)$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("Pearson X1", "Pearson X2", "F-test modelo"),
  Valor = c(round(r_x1, 4), round(r_x2, 4), format(p_modelo, scientific = TRUE, digits = 4)),
  Resultado = c(resultado_x1, resultado_x2, resultado_modelo)
)

kable(tabla_test, align = "lrc") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Test	Valor	Resultado
Pearson X1	0.7735	Aceptado
Pearson X2	0.8576	Aceptado
F-test modelo	4.105e-48	Significativo

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

## R² del modelo: 0.9209

Las dos variables pasan el test de Pearson (|r| > 0.7) y el modelo explica mas del 90% (R2 > 0.9).

10. Restricciones

• Valido para:
  • X1 (AVG_PRODUCTION) desde 3499 hasta infinito positivo
  • X2 (YEARS_ACTIVE) desde 37.30 hasta infinito positivo
• La inferencia de la muestra a la poblacion (reales positivos) se respalda en el R2 > 90% dentro del rango observado, pero entre mas se aleje la prediccion de ese rango, menor es la confiabilidad de la estimacion, ya que la relacion podria dejar de comportarse de forma lineal a medida que cambia la variabilidad de los datos en funcion del tiempo (por ejemplo, un pozo puede agotarse o cambiar de fase de produccion).
• Se uso la mediana por año (sin quitar outliers), no sirve para predecir pozos individuales atipicos.

11. Estimacion

nuevos_x1 <- c(3499, 5000, 10000, 25000, 50000, 100000)
nuevos_x2 <- c(38, 40, 45, 50, 60, 70)

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

# solo estimaciones positivas
estimaciones <- estimaciones %>% filter(Produccion_Acumulada_Estimada > 0)
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. Conclusion

Con AVG_PRODUCTION y YEARS_ACTIVE se obtiene el mejor ajuste con regresion lineal multiple (un plano en 3D):

\[y = -145857.4 + 42.29425 x_1 + 3920.59 x_2\]

Por cada unidad de produccion promedio, y cambia en 42.29425 (con años activos constante); por cada año activo, y cambia en 3920.59 (con produccion promedio constante).

Con produccion promedio = 2,992.51 y años activos = 45, se estima y = 157,135.2.

El modelo explica 92.1% de la variabilidad de Y (R2), el resto (7.9%) es por otros factores.

Pearson: r = 0.773 (Aceptado) para X1, r = 0.858 (Aceptado) para X2. F-test: Significativo (p = 4.11e-48).