Taller 2: Regresión Lineal

Especialización en Estadística Aplicada - Fundación Universitaria Los Libertadores

24 de febrero de 2026

Nombre: Danuil Elias Dueñas Criado

Código: 202540002444

  • Encontrar (idealmente) mínimo cinco (5) variables continuas (columnas). En general mínimo veinte (20) observaciones (filas)
  • Relacionadas con su actividad profesional o con su experiencia profesional
  • Cada una debe tener información de mínimo dos (2) años si es Mensual o cuatro (4) años trimestral o equivalente. Seleccionar dos (2) de esas variables

Calcular:

  1. Ecuación de Regresión
  2. ANOVA
  3. Pruebas de normalidad: Histograma con curva normal
  4. Gráfica Normal QQ Plot
  5. Asimetría
  6. Curtosis
  7. Test de Shapiro
  8. Test de Kolmogorov
  9. Test Jarque-Bera
  10. Homocedasticidad (2)
  11. Independencia (2)
  12. Pronosticar los valores de Y
  13. Analizar los resultados

Carga de datos

En esta sección se realiza la carga de los datos disponibles, los cuales corresponden a registros eléctricos tomados a lo largo de un pozo petrolero. Para esto se almacena la base de datos del pozo en la variable logs:

# Importación de los registros de pozo
logs <- read.csv("data/WELL1.csv")
head(logs)
##        MD BITSIZE   CALI    DRHO      DTC DTS     GRC NPHI     PE RHOB   RT10
## 1 1520.00    8.75 9.1026 -0.5681 143.4967  NA 37.9366   NA 3.9224   NA 0.3666
## 2 1520.25    8.75 9.1354 -0.3431 143.6970  NA 42.2258   NA 3.9757   NA 0.3397
## 3 1520.50    8.75 9.1643 -0.1532 143.6766  NA 46.3579   NA 4.1935   NA 0.3662
## 4 1520.75    8.75 9.1934 -0.0245 146.1810  NA 49.6346   NA 4.4278   NA 0.4508
## 5 1521.00    8.75 9.2193  0.0344 148.2777  NA 51.2478   NA 4.6215   NA 0.5736
## 6 1521.25    8.75 9.2385  0.0511 150.2207  NA 50.8796   NA 4.6870   NA 0.6643
##     RT20   RT30   RT60   RT90      SP
## 1 0.3927 0.4371 0.6567 1.0936 27.5609
## 2 0.3511 0.3688 0.5507 0.8884 23.2524
## 3 0.3647 0.3626 0.5535 0.9355 18.6427
## 4 0.4333 0.4104 0.6566 1.2250 14.5209
## 5 0.5358 0.4897 0.8075 1.6205 11.2540
## 6 0.6334 0.5939 0.9826 1.9795  8.9685

Se calcula el tamaño de la base de datos ingresada:

dim(logs)
## [1] 27865    16

El set de datos contiene 27865 registros para 16 variables. Se observa que los registros fueron tomados cada 0.25 pies a lo largo de la trayectoria del pozo.

library(dplyr)
glimpse(logs)
## Rows: 27,865
## Columns: 16
## $ MD      <dbl> 1520.00, 1520.25, 1520.50, 1520.75, 1521.00, 1521.25, 1521.50,…
## $ BITSIZE <dbl> 8.75, 8.75, 8.75, 8.75, 8.75, 8.75, 8.75, 8.75, 8.75, 8.75, 8.…
## $ CALI    <dbl> 9.1026, 9.1354, 9.1643, 9.1934, 9.2193, 9.2385, 9.2614, 9.2955…
## $ DRHO    <dbl> -0.5681, -0.3431, -0.1532, -0.0245, 0.0344, 0.0511, 0.0531, 0.…
## $ DTC     <dbl> 143.4967, 143.6970, 143.6766, 146.1810, 148.2777, 150.2207, 15…
## $ DTS     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ GRC     <dbl> 37.9366, 42.2258, 46.3579, 49.6346, 51.2478, 50.8796, 49.9535,…
## $ NPHI    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ PE      <dbl> 3.9224, 3.9757, 4.1935, 4.4278, 4.6215, 4.6870, 4.6609, 4.5997…
## $ RHOB    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ RT10    <dbl> 0.3666, 0.3397, 0.3662, 0.4508, 0.5736, 0.6643, 0.7672, 0.9177…
## $ RT20    <dbl> 0.3927, 0.3511, 0.3647, 0.4333, 0.5358, 0.6334, 0.7508, 0.9205…
## $ RT30    <dbl> 0.4371, 0.3688, 0.3626, 0.4104, 0.4897, 0.5939, 0.7287, 0.9246…
## $ RT60    <dbl> 0.6567, 0.5507, 0.5535, 0.6566, 0.8075, 0.9826, 1.1800, 1.4620…
## $ RT90    <dbl> 1.0936, 0.8884, 0.9355, 1.2250, 1.6205, 1.9795, 2.2461, 2.5892…
## $ SP      <dbl> 27.5609, 23.2524, 18.6427, 14.5209, 11.2540, 8.9685, 7.5723, 6…

Es posible observar que todas la variables importadas corresponden a tipo numérico double.

Análisis Exploratorio

En el siguiente diccionario se presenta el nombre completo asociado a cada registro para robustecer los análisis posteriores:

Diccionario de Curvas

  • MD .ft : Profundidad

  • BITSIZE .in : Bit Size - Tamaño de broca

  • CALI .in : Caliper - diámetro de hueco

  • DRHO .g/cc : DensityCorr - Corrección de densidad

  • DTC .us/ft : sónico compresional

  • DTS .us/ft : sónico flexural

  • GRC .api : Gamma Ray corregido

  • NPHI .ft/ft : porosidad neutrón

  • PE . : Factor fotoeléctrico

  • RHOB .g/cc : Densidad

  • RT10 .ohmm : 10in resistividad

  • RT20 .ohmm : 20in resistividad

  • RT30 .ohmm : 30in resistividad

  • RT60 .ohmm : 60in resistividad

  • RT90 .ohmm : 90in resistividad

  • SP .mV : Potencial espontáneo

A continuación, se grafican los datos de registros importados. Para este tipo de aplicaciones los registros son convencionalmente visualizados en profundidad:

library(ggh4x) 
library(tidyr)

# Se obtiene un nuevo dataframe en formato long
df_long <- logs %>%
  pivot_longer(
    cols = -MD,
    names_to = "registro",
    values_to = "valor"
  )

# Gráficos de registros en formato vertical
ggplot(df_long, aes(x = valor, y = MD, color = registro)) +
  geom_path(na.rm = TRUE) +
  scale_y_reverse() +
  facet_grid(. ~ registro, scales = "free_x") + 
  # Sección para seleccionar diferentes escalas por registro
  facetted_pos_scales(
    x = list(
      registro %in% c("RT10","RT20","RT30","RT60","RT90") ~ scale_x_log10(
        position = "top", guide = "axis_logticks"),
      TRUE ~ scale_x_continuous(position = "top")
    )
  ) +
  theme_bw() +
  labs(
    title = "Registro de Pozos: Pozo 01",
    subtitle = "Intervalo 1468 - 8413 pies",
    x = NULL, 
    y = "Profundidad (MD)"
  ) +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none",
    strip.text = element_text(face = "bold"),
    panel.grid.minor = element_line(color = "gray95"),
    axis.text.x = element_text(size = 5),
    panel.spacing = unit(0.2, "lines")
  ) +
  labs(x = NULL, y = "Profundidad (MD)")

A partir de la figura se observa que el registro sónico de tensión DTS exhibe una sección no registrada, así como los registros de porosidad neutrón NPHI y densidad bulk RHOB. Así mismo, es posible identificar que el registro de tamaño de broca BIT SIZE no proporciona información relevante a lo largo del pozo, lo cual indica que todo el pozo fue perforado con el mismo diámetro. A continuación, se obtienen los estadísticos básicos de las variables.

library(skimr)

# Estadísticos básicos de los registros 
skim(logs)
Data summary
Name logs
Number of rows 27865
Number of columns 16
_______________________
Column type frequency:
numeric 16
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
MD 0 1.00 5003.00 2011.02 1520.00 3261.50 5003.00 6744.50 8486.00 ▇▇▇▇▇
BITSIZE 0 1.00 8.75 0.00 8.75 8.75 8.75 8.75 8.75 ▁▁▇▁▁
CALI 314 0.99 9.72 1.62 3.62 8.68 9.02 10.21 23.23 ▁▇▁▁▁
DRHO 422 0.98 0.10 0.07 -0.57 0.04 0.09 0.15 0.31 ▁▁▁▇▃
DTC 210 0.99 104.60 19.65 58.34 90.27 100.42 117.75 173.46 ▂▇▅▂▁
DTS 7608 0.73 194.26 33.91 124.94 167.29 187.74 219.39 350.90 ▆▇▅▁▁
GRC 532 0.98 67.28 17.44 28.48 55.44 65.03 77.86 132.82 ▂▇▅▂▁
NPHI 1381 0.95 0.35 0.09 0.17 0.28 0.33 0.41 0.81 ▆▇▃▁▁
PE 422 0.98 2.45 0.33 1.25 2.25 2.43 2.64 4.69 ▁▇▃▁▁
RHOB 1343 0.95 2.24 0.15 1.34 2.16 2.27 2.34 2.57 ▁▁▂▇▅
RT10 55 1.00 4.61 2.49 0.34 3.10 3.91 5.51 33.91 ▇▁▁▁▁
RT20 55 1.00 4.44 2.50 0.35 2.83 3.63 5.36 30.19 ▇▂▁▁▁
RT30 55 1.00 4.47 2.76 0.36 2.70 3.52 5.44 28.07 ▇▂▁▁▁
RT60 55 1.00 4.51 3.01 0.55 2.62 3.45 5.41 27.17 ▇▂▁▁▁
RT90 55 1.00 4.58 3.25 0.89 2.59 3.44 5.44 29.79 ▇▂▁▁▁
SP 15 1.00 -10.58 16.56 -153.66 -18.16 -6.89 -0.15 27.56 ▁▁▁▆▇

A partir de los estadísticos básicos para cada variable es posible observar que:

  • MD/ Profundidad: Esta variable exhibe un rango entre 1520 y 8486 pies, con una mediana cercana a la media de 5003 pies, lo cual indica una distribución uniforme de la profundidad y un aspecto positivo para el análisis estratigráfico y de modelado.
  • BITSIZE / Tamaño de broca: Como se mencionó previamente, esta variable exhibe un comportamiento constante con un valor de 8.75, indicando que el pozo fue perforado con un solo tamaño de herramienta.
  • CALI / Tamaño de orificio: Exhibe una media de 9.72 in y valor medio de 23.22 in. Esto indica que el huevo está en promedio ligeramente sobrecalibrado (por encima del valor del bitsize), reflejando zonas de lavado severo, lo cual implica posibles afectaciones sobre registros RHOB y NPHI.
  • DRHO / Corrección de densidad: Exhibe una media de 0.10 indicando una buena calidad en general. Se identifican 422 valores faltantes.
  • DTC/DTS: Se observan medias de 104 µs/ft y 194 µs/ft, respectivamente. En el caso de DTS se identifican 7608 valores faltantes, lo cual puede limitar el análisis geomecánico en algunos intervalos del pozo. La dispersión de la variable sugiere una heterogeneidad litológica clara.
  • GRC/Gamma Ray : Se observa una media de 67 API en un rango entre 28 - 133 API. Lo cual indica una variabilidad estratigráfica clara con una alternancia arenas y arcillas marcada.
  • NPHI/ RHOB: Se observan medias de 0.35 y 2.23 g/cc, respectivamente. A su vez, se identifican 1381 datos faltantes para porosidad neutrón y 1343 para la densidad.
  • RT10:RT90: Medias entre 4 a 4.6 ohm.m, con máximos entre 30 y 33 ohm.m. Las variables exhiben pocos datos faltantes. Los valores son relativamente similares entre las diferentes profundidades de investigación lo cual sugiere una invasión probablemente moderada.
  • SP: Se observa una media cercana a -10mV con un mínimo de -153mV.

Con base en la cantidad significativa de datos, se decide remover las filas con datos faltantes para el análisis posterior:

library(tidyverse)
# Se crea una base de datos sin datos faltantes
logs_cleaned <- logs %>% drop_na()
nrow(logs_cleaned)
## [1] 19935

Posterior a la remoción de datos faltantes se cuenta con un total de 19935 observaciones para los registros eléctricos en el pozo. Ahora, se analizan las relaciones entre los registros importados:

# Dada la cantidad de datos, se seleccionan aleatoriamente 200 observaciones para obtener los gráficos de dispersión:
# Seleccionamos 200 números de fila al azar
set.seed(44)
indices <- sample(1:nrow(logs_cleaned), 200)

# Gráfico de pares
plot(logs_cleaned[indices, ])

A partir del gráfico de pares es posible observar que:

  • Los registros de resistividad (RT) exhiben una alta correlación, indicando zonas con saturación de fluidos similares y bajo radio de invasión

  • Los registros sónicos también exhiben una alta correlación positiva, esto es esperado dada la naturaleza de ambas mediciones.

El propósito del ejercicio es encontrar variables que permitan explicar el comportamiento del registro sónico DTC, excluyendo como variable independiente al registro DTS, por tanto inicialmente se analizará la relación entre el registro DTC y los registros GRC y RHOB, para esto se obtiene la siguiente gráfica y mapa de calor incluyendo los coeficientes de correlación de Pearson para cada caso:

# Se obtiene el gráfico de dispersión para la muestra seleccionada de 200 observaciones:
library(PerformanceAnalytics)
chart.Correlation(logs_cleaned[indices, c(5, 7, 10) ], histogram = TRUE, method = "pearson")

library(corrplot)
# Heatmap 
cor_matrix <- cor(logs_cleaned[indices, c(5, 7, 10)])
corrplot(cor_matrix, method = "color", )

Con base en el gráfico de dispersión y coeficientes de correlación visualizados en el mapa de calor, se observa que el registro sónico exhibe una correlación negativa con los registros de gamma ray y densidad. Para efectos del ejercicio, se implementará la regresión lineal simple empleando inicialmente el registro GRC como variable predictora.

Regresión Lineal Simple

1. Ecuación de Regresión

Se realiza la definición del modelo de regresión lineal simple entre el registro DTC y GRC

#Definición del modelo
modelo_rls <- lm(DTC ~ GRC , data = logs_cleaned)

# Resumen de resultados 
summary(modelo_rls)
## 
## Call:
## lm(formula = DTC ~ GRC, data = logs_cleaned)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.187  -8.756  -1.292   8.181  41.171 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 114.259155   0.445268  256.61   <2e-16 ***
## GRC          -0.276657   0.006604  -41.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.24 on 19933 degrees of freedom
## Multiple R-squared:  0.08091,    Adjusted R-squared:  0.08086 
## F-statistic:  1755 on 1 and 19933 DF,  p-value: < 2.2e-16

El modelo de regresión obtenido se representa mediante la siguiente ecuación:

\[\hat{y} = \beta_{0} - \beta_{1} x\]

\[DTC = 114.25 - 0.277 \times GRC\]

A partir del coeficiente \(\beta_{1}\) es posible interpretar que por cada incremento en 1 API en el registro de gamma ray, la medida de DTC disminuye en promedio 0.277 \(\mu\)s/ft, esto significa que entre más arcillosa sea la formación, menor es el tiempo registrado. Dado el t-value y el valor-p, es posible inferir que existe evidencia estadística que indica que la relación lineal entre las variables es significativa. Sin embargo, dado que el coeficiente de determinación es solo del 8.1%, la variabilidad de DTC explicada por GRC es muy baja, y por tanto, el carácter predictivo del modelo es extremadamente limitado, y evidencia la necesidad de incluir otras variables o tipos de relaciones para mejorar la captura del comportamiento de la variable objetivo, indicando que el modelo univariado parece ser insuficiente.

# Intervalo de confianza
confint(modelo_rls, level = 0.95)
##                   2.5 %      97.5 %
## (Intercept) 113.3863922 115.1319174
## GRC          -0.2896021  -0.2637114

La significancia de la variable explicativa es también observada mediante el intervalo de confianza de \(\beta_{1}\) entre [-0.289, -0.2637], al no contener el cero soporta la observación que la relación lineal entre las variables es significativa.

En la siguiente figura se presenta la regresión obtenida y la ecuación de la recta asociada.

# install.packages("ggpubr")
library(ggpubr)

ggplot(logs_cleaned, aes(x = GRC, y = DTC)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", formula = y ~ x) +
  theme_bw() +
  stat_regline_equation(label.y = 60, aes(label = ..eq.label..)) + # Muestra la ecuación
  stat_regline_equation(label.y = 65, aes(label = ..rr.label..))   # Muestra R²

2. Análisis de Varianza: ANOVA

# Pruebas de hipótesis
# ANOVA
anova(modelo_rls)
## Analysis of Variance Table
## 
## Response: DTC
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## GRC           1  263024  263024  1754.7 < 2.2e-16 ***
## Residuals 19933 2987888     150                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir de los resultados del análisis de varianza es posible calcular nuevamente la proporción explicada, así:

\[R^{2} = \frac{SCR}{SCT} = \frac{263024}{2987888 + 263024} = 0.08091\]

Como se mencionó previamente, el modelo solo explica una pequeña fracción de la variabilidad total de DTC.

Por otra parte, se obtuvo un \(valor-p\) de 2.2e-16, por tanto, a un nivel de confianza del 5%, se rechaza la hipótesis nula (\(H_0: \beta_{GRC} = 0\)) y se infiere que GRC está asociado linealmente con DTC.

El análisis ANOVA confirma que existe relación estadística entre las variables, sin embargo, el comportamiento de DTC depende claramente de más variables.

A continuación, se presentan las pruebas para la validación de supuestos para que el modelo de regresión lineal simple \(y_i = \beta_0 + \beta_1 x_i + \epsilon_i\) sea válido. Estos supuestos son:

  • Linealidad: La relación entre la variable dependiente (\(y\)) y la independiente (\(x\)) debe ser lineal en sus parámetros. \[\mathbb{E}(y_i) = \beta_0 + \beta_1 x_i\]

  • Media cero de los errores: El valor esperado del término de error es cero para cualquier valor de la variable independiente. \[\mathbb{E}(\epsilon_i) = 0\]

  • Homocedasticidad (Varianza Constante): La varianza de los errores debe ser constante para todas las observaciones de \(x\). Si la varianza no es constante, existe heterocedasticidad. \[Var(\epsilon_i) = \sigma^2\]

  • Independencia (No Autocorrelación): Los errores de las observaciones son independientes entre sí; no existe correlación entre errores de diferentes niveles de \(x\). \[Cov(\epsilon_i, \epsilon_j) = 0, \quad \forall i \neq j\]

  • Normalidad: Para realizar inferencia estadística (pruebas \(t\) y \(F\)), se asume que los errores siguen una distribución normal. \[\epsilon_i \sim N(0, \sigma^2)\]

3. Pruebas de Normalidad: Histograma con curva normal

En esta sección se realiza la validación visual de la normalidad de los errores al comparar los residuales del modelo con la curva normal teórica.

# Gráfico conjunto de Histograma de datos y curva normal
hist(resid(modelo_rls), prob = TRUE, 
     main = "Histograma con curva normal", ylab = "Frecuencia Relativa")
x <- seq(min(resid(modelo_rls)), max(resid(modelo_rls)), length = 40)
f <- dnorm(x, mean = mean(resid(modelo_rls)), sd=sd(resid(modelo_rls)))
lines(x, f, col = "red", lwd=2)

Con base en la figura es posible observar que:

  • La distribución de los errores se encuentra aproximadamente centrada alrededor del cero, lo cual nos indica que el modelo no tiene un sesgo importante en sus predicciones.

  • El histograma exhibe una forma de campana con una simetría importante. No se observa un sesgo pronunciado hacia ningún lado, lo cual puede sugerir que los errores en ambas direcciones ocurren con frecuencias similares.

  • La curva roja muestra un ajuste aceptable a la mayoría de las barras del histograma. Sin embargo, se nota una ligera mayor concentración de datos alrededor de -10 en comparación con la curva teórica y colas se reducen gradualmente hacia +/-40.

  • Visualmente, los residuos parecen seguir una distribución normal, lo cual implica que los intervalos de confianza resultantes del modelo son válidos.

Estos hallazgos son revaluados en secciones posteriores mediante pruebas estadísticas.

4. Gráfica Normal QQ Plot

# Verificación de supuestos de normalidad
par(mfrow = c(2, 2))
plot(modelo_rls)

par(mfrow = c(1, 1))
# Gráfico Q-Q-norm
qqnorm(resid(modelo_rls))
qqline(resid(modelo_rls), col = "red")

Con base en los gráficos de diagnóstico de residuales, es posible realizar las siguiente observaciones:

  • Residuales vs. Valores Predichos: La línea roja se encuentra aproximadamente horizontal y está cerca del cero, lo cual sugiere que el supuesto de linealidad se cumple. No se observa algún patrón claro, aunque se nota una gran dispersión de puntos entre los valores ajustados de 90 y 100.

  • Normal Q-Q plot: Los puntos se ajustan de forma cercana a la línea roja en la mayor parte. Sin embargo, los extremos exhiben cierta desviación. La figura indica una distribución aproximadamente normal con una ligera desviación en los valores extremos.

  • Homocedasticidad (Gráfico Scale-Location): La línea roja es bastante plana, lo cual indica un comportamiento cercano a la homocedasticidad. La dispersión de los errores estandarizados no exhibe un aumento o disminución sistemática con los valores ajustados.

  • Residuales vs. Apalancamiento (puntos influyentes): La mayoría de los datos tienen un bajo nivel de apalancamiento. No se observan valores atípicos influyentes que deformen de forma importante los resultados del modelo.

En general, el modelo cumple razonablemente con los supuestos de la regresión lineal simple para inferencia.

5. Asimetría

# MÉTODOS ANALÍTICOS
# Coeficiente de Asimetría
library(moments)
skewness(resid(modelo_rls))
## [1] 0.3232495

El valor obtenido es ligeramente positivo, lo cual indica una asimetría leve a la derecha. Para fines prácticos es posible considerar como simétrico. Esta evidencia confirma lo observado en el histograma, la distribución de los errores exhibe un equilibrio alrededor de la media de cero.

6. Curtosis

# Coeficiente de curtosis
kurtosis((resid(modelo_rls)))
## [1] 2.746562

El valor obtenido es muy cercano a 3. Al ser ligeramente menor a 3, la distribución es técnicamente platicúrtica, indicando que tiene colas un poco más delgadas y un centro ligeramente más achatado que la normal teórica.

Teniendo en cuenta que la asimetría es casi nula y que la curtosis se aproxima significativamente a 3, es posible concluir que los residuales del modelo siguen una distribución aproximadamente normal.

Adicionalmente, dada la alineación de los puntos en el gráfico Q-Q, se valida que las inferencias estadísticas del modelo son confiables.

A continuación, se realizan las pruebas estadísticas de normalidad:

7. Test de Shapiro

Dado el tamaño de la muestra empleada excede el umbral para el uso de la prueba de Shapiro-Wilk, se define un modelo con una nueva muestra reducida:

# Test de Shapiro-Wilk tamaño de la muestra es menor de 50

# Seleccionamos 50 números de fila al azar
set.seed(44)
indices <- sample(1:nrow(logs_cleaned), 50)
logs_cleaned_red = logs_cleaned[indices, ]
modelo_rls_red = lm(formula = DTC ~ GRC, data = logs_cleaned_red)
shapiro.test(resid(modelo_rls_red))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(modelo_rls_red)
## W = 0.98452, p-value = 0.7506

La prueba de Shapiro-Wilk es un contraste de hipótesis diseñado para evaluar si una muestra de datos proviene de una población con distribución normal. Recomendada para muestras de tamaño moderado (\(n < 5000\)).

Definición de Hipótesis

Para el conjunto de residuales del modelo ajustado, se plantean las siguientes hipótesis estadísticas:

  • Hipótesis Nula (\(H_0\)): Los residuales siguen una distribución normal (\(\epsilon_i \sim N(0, \sigma^2)\)).
  • Hipótesis Alterna (\(H_1\)): Los residuales no siguen una distribución normal.

A partir de la ejecución de la prueba sobre los residuales del modelo, se observa que, considerando un nivel de significancia estándar de \(\alpha = 0.05\), el p-valor (0.7506) es significativamente mayor que \(\alpha\). Por tanto, no se rechaza la hipótesis nula (\(H_0\)) y existe evidencia estadística suficiente para afirmar que los residuales del modelo siguen una distribución normal.

El resultado del p-valor es coherente con los valores de asimetría (0.32) y curtosis (2.74), los cuales se encuentran dentro de los rangos esperados para una distribución normal (cercanos a 0 y 3, respectivamente).

8. Test de Kolmogorov

# Test de Kolmogorov-Smirnov-Lilliefors tamaño de la muestra es mayor de 50
library(nortest)
lillie.test(resid(modelo_rls))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  resid(modelo_rls)
## D = 0.044402, p-value < 2.2e-16

La prueba de Lilliefors (Kolmogorov-Smirnov) es una adaptación de la prueba de Kolmogorov-Smirnov diseñada para contrastar la normalidad de una muestra cuando los parámetros de la población (media \(\mu\) y varianza \(\sigma^2\)) son desconocidos y deben ser estimados a partir de los datos. Es especialmente útil en conjuntos de datos grandes (\(n > 50\)).

Hipótesis del Contraste

Se plantean las siguientes hipótesis:

  • \(H_0\): Los residuales siguen una distribución normal (\(\epsilon_i \sim N(\mu, \sigma^2)\)).
  • \(H_1\): Los residuales no siguen una distribución normal.

A partir de la ejecución de la prueba sobre los residuales del modelo, se obtuvieron los siguientes valores: El valor \(D = 0.0444\) representa la máxima diferencia absoluta entre la función de distribución acumulada de la muestra y la distribución normal teórica. Un valor tan bajo sugiere una forma visualmente cercana a la normal. Dado que el p-valor (< 2.2e-16) es extremadamente pequeño y mucho menor que el nivel de significancia, \(\alpha = 0.05\), Se rechaza la hipótesis nula (\(H_0\)).

A pesar de que el estadístico \(D\) es pequeño, el volumen masivo de datos (27,865 observaciones) otorga a la prueba un poder estadístico tan alto que cualquier desviación mínima de la normalidad perfecta resulta en un rechazo estadístico.

9. Test Jarque-Bera

# Test de Jarque-Bera
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
jarque.bera.test(resid(modelo_rls))
## 
##  Jarque Bera Test
## 
## data:  resid(modelo_rls)
## X-squared = 400.52, df = 2, p-value < 2.2e-16

La prueba de Jarque-Bera (JB) es un contraste de bondad de ajuste que determina si los datos de la muestra tienen la asimetría y la curtosis correspondientes a una distribución normal. A diferencia de otras pruebas, esta se basa específicamente en los momentos de la distribución.

Hipótesis del Contraste

  • \(H_0\): Los residuales tienen una asimetría de 0 y una curtosis de 3 (Distribución Normal).
  • \(H_1\): Los residuales no presentan una combinación de asimetría y curtosis normal.

A partir de los residuales del modelo, se observa que el valor del estadístico de prueba es significativamente alto. Este estadístico mide qué tan lejos están la asimetría (0.32) y la curtosis (2.74) de los valores teóricos (0 y 3). Por su parte, el p-valor (< 2.2e-16) es extremadamente pequeño y mucho menor al nivel de significancia \(\alpha = 0.05\). Por tanto, existe evidencia estadística para rechazar la hipótesis nula (\(H_0\)).

Tras aplicar las pruebas de Shapiro-Wilk (en muestra reducida), Lilliefors y Jarque-Bera, se observa una discrepancia técnica común en el análisis de grandes bases de datos (n = 27,865): Todas las pruebas aplicadas al dataset completo rechazan la normalidad estricta debido al gran poder estadístico (sensibilidad) ante muestras masivas.

A pesar del rechazo formal de las pruebas de hipótesis, la distribución de los residuales es lo suficientemente cuasi-normal para no invalidar el modelo.

10. Homocedasticidad (2)

##### Pruebas de Homocedasticidad
# Prueba de Heterocedasticidad, Ho = varianzas son constantes
#bptest
library(lmtest)
bptest(modelo_rls)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_rls
## BP = 78.037, df = 1, p-value < 2.2e-16

La prueba de Breusch-Pagan se emplea para validar el supuesto de varianza constante, homocedasticidad.

Definición de Hipótesis

  • \(H_0\): Existe homocedasticidad (Varianza constante de los errores).
  • \(H_1\): Existe heterocedasticidad (La varianza de los errores no es constante).

Al asumir un nivel de significancia (\(\alpha\)) de \(0.05\), y dado que el p-valor (0.0007119) es significativamente menor que \(\alpha\), se rechaza la hipótesis nula (\(H_0\)).

A partir de la prueba estadística, se concluye que el modelo presenta heterocedasticidad. Esto significa que la variabilidad de los errores cambia a medida que cambian los valores de la variable independiente. Esta heterocedasticidad puede afectar la precisión de los errores estándar y los intervalos de confianza.

#ncvtest
#install.packages("car")
library(car)
ncvTest(modelo_rls)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 68.1478, Df = 1, p = < 2.22e-16

La prueba de puntuación de varianza no constante (ncvTest) también es empleada para validar el supuesto de homocedasticidad. Este contraste evalúa si la varianza de los residuos es dependiente de los valores ajustados del modelo (\(\hat{y}\)).

  • \(H_0\): Homocedasticidad (Varianza constante de los errores: \(\sigma_i^2 = \sigma^2\)).
  • \(H_1\): Heterocedasticidad (Varianza no constante: \(\sigma_i^2 = \sigma^2 f(\hat{y}_i)\)).

Dado que el p-valor (< 2.22e-16) es pequeño y asumiendo un nivel de significancia (\(\alpha = 0.05\)), se rechaza la hipótesis nula (\(H_0\)). Por tanto, existe evidencia estadística suficiente para afirmar que los residuales presentan heterocedasticidad.

11. Independencia (2)

##### Pruebas de Independencia
dwtest(modelo_rls)
## 
##  Durbin-Watson test
## 
## data:  modelo_rls
## DW = 0.013267, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
bgtest(modelo_rls) ## Durbin - Watson test
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_rls
## LM test = 19672, df = 1, p-value < 2.2e-16

El supuesto de independencia establece que los errores de las observaciones deben ser aleatorios y no presentar una relación sistemática entre sí. Para su verificación, se aplican dos pruebas complementarias.

Prueba de Durbin-Watson (DW)

Evalúa la autocorrelación de primer orden (\(\rho\)) en los residuales del modelo.

  • \(H_0\): No existe autocorrelación (\(\rho = 0\)).
  • \(H_1\): Existe autocorrelación positiva (\(\rho > 0\)).

El estadístico \(DW\) se encuentra extremadamente cerca de 0 (en una escala de 0 a 4, donde 2 es independencia total). Dado que el p-valor (< 2.2e-16) es mucho menor a \(\alpha = 0.05\), se rechaza la hipótesis nula, por tanto se puede observar que existe una fuerte autocorrelación positiva.

Prueba de Breusch-Godfrey (LM Test)

A diferencia de DW, esta prueba es más general y robusta ante la falta de normalidad de los residuos.

  • \(H_0\): No existe correlación serial de orden 1.
  • \(H_1\): Existe correlación serial de orden 1.

El alto valor del estadístico de Lagrange (19672) y un p-valor (< 2.2e-16) confirman nuevamente el rechazo de la hipótesis nula. Por tanto, los errores del modelo están correlacionados serialmente.

El rechazo simultáneo de la independencia por ambas pruebas confirma que el modelo de regresión lineal simple presenta autocorrelación serial. Este fenómeno es intrínseco a los datos de pozos: los valores de los registros (como DTC o GRC) cambian gradualmente con la profundidad, lo que significa que el error en un punto está fuertemente ligado al error en el punto inmediatamente superior.

La autocorrelación causa que los errores estándar sean subestimados.

12. Pronosticar los valores de Y

En esta sección se pronostican los valores de DTC empleando los valores de GRC y el modelo de regresión definido.

# Pronósticos

# Se crean un dataframe con solo los datos de GRC y se emplean en la predicción. 
nuevo <- data.frame(GRC = logs_cleaned$GRC)
predictions = predict(modelo_rls, newdata = nuevo)

# Se presentan los diez primeros valores predichos.
predictions[0:10]
##        1        2        3        4        5        6        7        8 
## 94.35998 94.62126 95.07984 95.14555 94.96116 94.61735 93.81292 93.17570 
##        9       10 
## 92.72275 92.87956

A continuación, se presenta la gráfica de los recta estimada y el intervalo de predicción para el rango de datos de GRC en la base de datos empleada.

grid_datos <- data.frame(GRC = seq(min(logs_cleaned$GRC, na.rm=T), 
                                  max(logs_cleaned$GRC, na.rm=T), 
                                  length.out = 100))

predicciones <- predict(modelo_rls, newdata = grid_datos, 
                        interval = "prediction", level = 0.95)

df_grafico <- cbind(grid_datos, predicciones)

ggplot(df_grafico, aes(x = GRC, y = fit)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "lightblue", alpha = 0.4) +
  geom_point(data = logs_cleaned, aes(x = GRC, y = DTC), alpha = 0.5, size = 1) +
  geom_line(color = "red", linewidth = 1) +
  theme_bw() +
  labs(
    title = "Predicción de DTC vs GRC",
    subtitle = "Sombreado: Intervalo de Predicción al 95%",
    y = expression(DTC ~ (mu*s/ft)),
    x = expression(GRC ~ (API))
  )

13. Analizar resultados

El modelo de regresión lineal simple desarrollado permitió identificar una relación inversa y estadísticamente significativa (\(p < 2.2e-16\)) entre el registro sónico (DTC) y el gamma ray (GRC), donde el tiempo de tránsito disminuye en promedio \(0.277 \mu s/ft\) por cada unidad API de incremento en la arcillosidad. Sin embargo, la capacidad predictiva del modelo es limitada, explicando únicamente el 8.1% de la variabilidad total de los datos (\(R^2 = 0.081\)).

El análisis de supuestos reveló que, aunque los residuales presentan una normalidad aproximada respaldada por valores de asimetría (\(0.32\)) y curtosis (\(2.74\)) cercanos a los parámetros teóricos, las pruebas de ncvTest y Durbin-Watson confirmaron la presencia de heterocedasticidad y una fuerte autocorrelación serial (\(DW = 0.013\)). Estos resultados sugieren que, si bien el modelo es útil para describir la tendencia general del pozo, la naturaleza secuencial y compleja de los registros geofísicos requiere la inclusión de variables adicionales, como densidad o porosidad, para lograr un pronóstico técnicamente robusto.