Diagnóstico del modelo lineal: trees

Author

Jorge Che, Cristina Cruz

1. Introducción

Se presenta un análisis de regresión simple para modelar la relación entre dos variables cuantitativas. Se utilizó la base de datos trees disponible en R, el cual contiene mediciones del diámetro, altura y volumen de madera de 31 cerezos negros talados.

Se analizará la relación entre la altura ( \(X=\text{Height}\) , en pies) sobre el volumen del árbol ( \(Y=\text{Volume}\) , en pies cúbicos).

2. Análisis Exploratorio de Datos (AED)

Previo al ajuste del modelo, se obtienen datos descriptivos de las variables de interés.

Variable explicativa \(X\): Altura (Height)

favstats(~ Height, data = trees)

 min Q1 median Q3 max mean       sd  n missing
  63 72     76 80  87   76 6.371813 31       0

Esta variable presenta un comportamiento altamente simétrico. La media y la mediana coinciden exactamente, lo que indica que la distribución no presenta sesgos marcados hacia ninguno de los extremos.

La dispersión es relativamente baja, lo cual indica que las alturas están concentradas alrededor de su promedio. No posee datos faltantes.

Variable respuesta \(Y\): Volumen (Volume)

favstats(~ Volume, data = trees)

  min   Q1 median   Q3 max     mean       sd  n missing
 10.2 19.4   24.2 37.3  77 30.17097 16.43785 31       0

Presenta una distribución sesgada a la derecha. La media se encuentra notablemente por encima de la mediana, probablemente influencia por la presencia de valores atípicos superiores.

Muestra una variablidad alta, extendiéndose entre 10.2 y 77 pies cúbicos y una desviación estándar considerable. No posee datos faltantes.

Gráfica de dispersión

ggplot(trees, aes(x = Height, y = Volume)) +
  geom_point(color = "#16a085", size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se = TRUE, color = "#2c3e50", fill = "#16a085", alpha = 0.2) +
  labs(
    title = "Relación entre la altura y el volumen de los árboles",
    subtitle = "Dataset trees",
    x = "Altura (pies)",
    y = "Volumen (pies cúbicos)"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

El gráfico de dispersión sugiere una correlación positiva, pero evidencia un aumento de la variabilidad a medida que la altura crece. Para alturas bajas los puntos están muy agrupados cerca de la recta de regresión, mientras que a mayores alturas el volumen se dispersa sensiblemente (alcanzando el punto atípico de \(77\) en la esquina superior derecha).

3. Ajuste del modelo de Regresión Lineal Simple

El modelo matemático poblacional se define como: \[y_i=\beta_0+\beta_1x_i+\epsilon_i\]

Donde:

\(Y_i\) es el volumen (Volume).
\(X_i\) es la altura (Height).
\(B_0\) es el intercepto.
\(B_1\) es la pendiente.
\(\epsilon_i\) es el término de error aleatorio.

Estimación de los coeficientes

Ajustamos el modelo en R utilizando la función lm():

modelo <- lm(Volume ~ Height, data = trees) 

summary(modelo)


Call:
lm(formula = Volume ~ Height, data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.274  -9.894  -2.894  12.068  29.852 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.1236    29.2731  -2.976 0.005835 ** 
Height        1.5433     0.3839   4.021 0.000378 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.4 on 29 degrees of freedom
Multiple R-squared:  0.3579,    Adjusted R-squared:  0.3358 
F-statistic: 16.16 on 1 and 29 DF,  p-value: 0.0003784

Interpretación de los parámetros

Intercepto \(\beta_0\): -87.1236 Representa el valor estimado del volumen \(Y\) cuando la altura del árbol (\(X\)) es igual a \(0\) pies.

En la práctica, no tiene un sentido real, ya que un árbol no puede tener un volumen negativo ni una altura de cero pies.

Dado que \(\text{p-value}=0.005835< \alpha =0.05\), el intercepto es estadísticamente distinto de cero.
Pendiente \(\beta_1\): 1.5433 Por cada pie adicional que aumenta la altura de un árbol (\(1\) \(\text{ft}\)), se estima que su volumen aumenta en promedio \(1.5433\) pies cúbicos (\(\text{ft}^3\)).

Dado que \(\text{p-value}=0.000378< \alpha =0.05\), se infiere que existe una relación lineal positiva y estadísticamente significativa entre la altura y el volumen de los árboles.

4. Diagnóstico del modelo

Gráficos de diagnóstico estándar

par(mfrow = c(2, 2))
plot(modelo)

Análisis detallados de los supuestos

Normalidad de los errores (Q-Q normal): En el gráfico, la mayoría de los residuales se alinean a la recta teórica. Únicamente se observan ligeras desviaciones en los extremos. Esto indica que posiblemente se cumpla el supuesto de normalidad de los residuales.
Puntos influyentes (Residuales vs Leverage): La mayoría de las observaciones presentan valores bajos de Leverage (apalancamiento) y residuos acotados, sin embargo, la observación 31 se unica cercano al límite de la Distancia de Cook (0.5). Esto indica que la observación 31 es un punto influyente ya que ejerce un peso considerable sobre la estimación de la pendiente de la recta de regresión.
Homocedasticidad (Scale-Location): La línea roja presenta una marcada tendencia ascendente a lo largo de los valores ajustados. Es decir, que a medida que el volumen predicho aumenta, la magnitud de los residuos estandarizados también incrementa significativamente. Por lo tanto, esto podría indicar que el supuesto de homocedasticidad no se cumple.
Independencia (Residuales vs valores ajustados): Los residuos no se distribuyen de manera completamente aleatoria alrededor de la línea cero. Se aprecia una curvatura leve y un patrón de dispersión en forma de “embudo”, donde los errores se vuelven más variables en valores predichos altos (alrededor de \(35\) a \(45\)). Esto sugiere que el supuesto de independencia no se cumple.

5. Análisis de la variable respuesta

En este apartado se pretende realizar un análisis exploratorio de las variables volumen, en el cual se realizan diversos gráficos e indicadores para observar el comportamiento y la distribución de cada variable.

grafico_1<-ggplot(data=trees)+
  geom_histogram(aes(x=Volume),binwidth = 2,fill = "#4D6291", color = "black")+
  theme(axis.text.x = element_text(size = 10))+
  labs(title = "Histograma variable volumen", x = "Volumen en pies cúbicos", y = "Frecuencia")+
  theme(plot.title = element_text(hjust = 0.5))

grafico_1

En relación con los cuartiles, encontramos que el 25% de los árboles tienen un volumen menor o igual a 19.4 pies cúbicos, el 50% de los árboles tienen un volumen menor o igual a 24.2 pies cúbicos y por último el 75% de los árboles tienen un volumen menor o igual a 37.3 pies cúbicos.

# 1. Creamos la estructura Vertical (Indicador vs Valor)

tabla_vertical <- data.frame(
  Valor = quantile(trees$Volume,prob=c(0.25,0.5,0.75)))

# 2. Presentación con Kable
kable(
  tabla_vertical,
  digits = 2, # Redondea automáticamente a 2 decimales
  caption = "Indicadores de posición",
  align = c("l", "c"),
  col.names = c("Cuartil", "Valor")
)

Indicadores de posición
Cuartil	Valor
25%	19.4
50%	24.2
75%	37.3

Esta representación de lo cuartiles se evidencia de una mejor forma en el siguiente gráfico, donde podemos afirmar lo anteriormente dicho y anexándole que hay una pequeña cantidad de árboles con volumenes grandes (puntos rojos):

grafico_2<-ggplot(data=trees)+
  geom_boxplot(aes(x=Volume), col= "blue", outlier.color = "red")+
  labs(title = "Diagrama de cajas - Variable volumen", x = "Pies cúbicos", y="Porcentaje")+
  theme(plot.title = element_text(hjust = 0.5))
grafico_2

La desviación de los volúmenes varía por encima o por debajo de la media en 16.44 pies cúbicos. Por otro lado, el coeficiente de variación (CV = 54.48%) indica que el volumen de los árboles presenta una variabilidad bastante alta alrededor de su media (30.17 pies cúbicos). Esto refleja que los datos muestran heterogeneidad en su distribución.

En cuanto a la forma de la distribución, el coeficiente de curtosis es de 0.25. Al presentar un valor positivo, la distribución se clasifica como levemente leptocúrtica, lo que indica que la concentración de los datos alrededor de la media es ligeramente mayor a la de una distribución normal, presentando una forma relativamente más larga y alta.

# 1. Obtenemos los valores 
promedio <- mean(trees$Volume)
desviacion <- sd(trees$Volume)

# 2. Creamos la estructura Vertical (Indicador vs Valor)

tabla_vertical <- data.frame(
  Indicador = c("Moda", "Varianza", "Desviación Estándar", "Coeficiente de Variación (%)", "Curtosis"),
  Valor = c(
    mfv(trees$Volume)[1],
    var(trees$Volume),
    desviacion,
    (desviacion / promedio) * 100,
    kurtosis(trees$Volume)
  )
)

# 3. Presentación con Kable
kable(
  tabla_vertical,
  digits = 2, # Redondea automáticamente a 2 decimales
  caption = "Indicadores Descriptivos de Volumen",
  align = c("l", "c"),
  col.names = c("Estadístico", "Valor")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE, position = "center")

Indicadores Descriptivos de Volumen
Estadístico	Valor
Moda	10.30
Varianza	270.20
Desviación Estándar	16.44
Coeficiente de Variación (%)	54.48
Curtosis	0.25

Se evalua la presencia de valores atípicos mediante el método del Rango Intercuartílico (1.5×IQR1.5×IQR). Los datos oscilan entre un mínimo de -7.45 pies cúbicos y un máximo de 64.15 pies cúbicos.

# Cálculo de cuartiles y Rango Intercuartílico (IQR)
Q1 <- quantile(trees$Volume, 0.25)
Q3 <- quantile(trees$Volume, 0.75)
IQR_val <- IQR(trees$Volume)

# Límites para detectar atípicos
limite_inferior <- Q1 - 1.5 * IQR_val
print(limite_inferior)

  25% 
-7.45

limite_superior <- Q3 + 1.5 * IQR_val
print(limite_superior)

  75% 
64.15

Se detectó la presencia de un valor atípico (outlier) en la variable volumen, siendo este de 77 pies cúbicos.

# Identificar valores fuera de los límites
atipicos <- trees$Volume[trees$Volume < limite_inferior | trees$Volume> limite_superior]
atipicos

[1] 77

6. Supuesto de Normalidad de los residuales

Planteamiento de las Hipótesis

Para cualquier prueba de normalidad sobre los residuales del modelo:

\(H_0\) (Hipótesis Nula): Los residuales del modelo siguen una distribución normal
\(H_1\) (Hipótesis Alternativa): Los residuales del modelo no siguen una distribución normal

Decisión (⍺=0.05):
Si p-valor > 0.05: No se rechaza \(H_0\) . Se concluye que existe evidencia suficiente de que se cumple el supuesto de normalidad en los residuales.
Si p-valor < 0.05: Se rechaza \(H_0\) . Existe evidencia de que los residuales no se distribuyen de manera normal.

Debido a que el número de observaciones es 31, se utilizó la prueba Shapiro-Wilk para probar la normalidad.

# 3. Prueba Principal: Shapiro-Wilk (Recomendada para n < 50)
shapiro_res <- shapiro.test(modelo$residuals)
print(shapiro_res)


    Shapiro-Wilk normality test

data:  modelo$residuals
W = 0.95078, p-value = 0.164

Conclusion de la prueba de normalidad

Como p-valor= 0.164, es mayor al nivel de significancia previamente establecido ⍺=0.05, por lo que no se rechaza la hipótesis nula \(H_0\). Por lo tanto, con una confianza del 95%, existe evidencia estadística suficiente para afirmar que los residuales del modelo satisfacen el supuesto de normalidad.

7. Supuesto de Homocedasticidad

Planteamiento de las Hipótesis

Para la prueba de varianza constante sobre los residuales del modelo ( \(\sigma_i^2\) \(\neq\) \(\sigma^2\)) se tiene:

\(H_0\) (Hipótesis Nula): Los residuales son homocedásticos (la varianza de los errores es constante)
\(H_1\) (Hipótesis Alternativa): Los residuales son heterocedásticos (la varianza de los errores no es constante)

Decisión (⍺=0.05):
Si p-valor > 0.05: No se rechaza \(H_0\) . Se concluye que existe evidencia suficiente de que se cumple el supuesto de homocedasticidad.
Si p-valor < 0.05: Se rechaza \(H_0\) . Existe evidencia de que los residuales no tienen varianza constante.

prueba_bp <- bptest(modelo)
print(prueba_bp)


    studentized Breusch-Pagan test

data:  modelo
BP = 12.207, df = 1, p-value = 0.0004762

Conclusiones sobre varianza constante

Mediante la prueba Breusch-Pagan se obtuvo un estadístico 𝒳²=12.207, junto con un p-valor= 0.0004. Al ser el p-valor menor que la significancia ⍺=0.05, existe evidencia suficiente para rechazar la hipótesis nula, siendo así que los residuales son heterocedásticos.

8.Supuesto de independencia

Para la prueba Durbin-Watson:

\(H_0\) (Hipótesis Nula): ⍴=0, hay independencia (No autocorrelación entre residuales)
\(H_1\) (Hipótesis Alternativa): ⍴≠0, no hay independencia (Autocorrelación entre residuales)
Si p-valor > 0.05: No se rechaza \(H_0\) . Se concluye que existe evidencia suficiente de que se cumple el supuesto independencia (No autocorrelación)
Si p-valor < 0.05: Se rechaza \(H_0\) . Existe evidencia de que existe autocorrelación.

prueba_dw <- dwtest(modelo)
print(prueba_dw)


    Durbin-Watson test

data:  modelo
DW = 0.5009, p-value = 8.957e-08
alternative hypothesis: true autocorrelation is greater than 0

# 3. Gráfico de autocorrelación (ACF) y residuales en función del orden
datos_indep <- data.frame(
  Orden = 1:nrow(trees),
  Residuales = residuals(modelo)
)
# Gráfico de residuales seguidos por orden de observación
ggplot(datos_indep, aes(x = Orden, y = Residuales)) +
  geom_point(color = "#2b5c8f", size = 2.5) +
  geom_line(color = "#2b5c8f", linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#d9534f", linewidth = 1) +
  labs(
    title = "Residuales por Orden de Observación",
    subtitle = "Evaluación de Independencia de Errores",
    x = "Orden de la Observación",
    y = "Residuales"
  ) +
  theme_minimal()

Resultado prueba DW

Al ejecutar la prueba al modelo ajustado se obtuvo el estadístico DW=.5009 con un p-valor menor al nivel de significancia de 0.05. Rechazando la independencia entre residuales. Esto pudiese ocurrir al no considerar otras variables como la medida del diámetro de la circunferencia (pulgadas).

CONCLUSIONES

El modelo satisface el supuesto de normalidad de los residuales, sin embargo, la prueba Breusch-Pagan indica que hay evidencia significativa de heterocedasticidad entre los residuales. Asímismo el supuesto de independencia de los errores, tampoco se cumple debido a que la prueba Durbin-Watson presenta un p-valor menos al nivel de significancia (omisión de covariables relevantes como el diámetro de los árboles).