library(ggplot2)
#cargar las librerías necesarias
library(ggplot2)
library(MASS)
library(mosaic)
library(knitr)
library(kableExtra)
library(modeest)
library(e1071)
library(nortest)
library(lmtest)
library(palmerpenguins)Diagnóstico de Regresión Lineal y Efecto de la Winsorización en pingüinos machos Adelie
Análisis de Regresión Lineal: Diagnóstico de Valores Atípicos y Winsorización en pingüinos Adelie
En este estudio se analizará la relación entre la masa corporal y la longitud de aleta de los pingüinos Adelie, diferenciando entre machos y hembras. Se busca determinar si una transformación de la variable respuesta permite mejorar el cumplimiento de los supuestos del modelo de regresión lineal.
La transformación de Box-Cox permite estimar, mediante máxima verosimilitud, un parámetro \(\lambda\) que determina la potencia aplicada a la variable respuesta. Su propósito es encontrar una transformación que contribuya a mejorar las propiedades de los residuos, particularmente la normalidad y la homocedasticidad.
La transformación se define como:
\(Y^{(\lambda)}= \begin{cases} \dfrac{Y^\lambda-1}{\lambda}, & \lambda\neq 0\\ \ln(Y), & \lambda=0 \end{cases}\)
En nuestro caso, la variable longitud de aleta contiene únicamente valores positivos, por lo que es posible aplicar este procedimiento.
data("penguins")
#Eliminamos NA en penguis
data = subset(
penguins,
species == "Adelie" &
!is.na(sex) &
!is.na(flipper_length_mm ) &
!is.na(body_mass_g)
)
#Separar subconjubtos
adelie_machos <- subset(
data,
sex == "male")
adelie_hembras <- subset(
data,
sex == "female")Modelo de regresión linea para pingüinos machos
#MODELO DE REGRESIÓN LINEAL MACHOS
modelo_machos <- lm(
flipper_length_mm ~ body_mass_g+0,
data = adelie_machos
)
summary(modelo_machos)
Call:
lm(formula = flipper_length_mm ~ body_mass_g + 0, data = adelie_machos)
Residuals:
Min 1Q Median 3Q Max
-35.908 -8.355 1.645 12.014 32.478
Coefficients:
Estimate Std. Error t value Pr(>|t|)
body_mass_g 0.047292 0.000443 106.7 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.36 on 72 degrees of freedom
Multiple R-squared: 0.9937, Adjusted R-squared: 0.9936
F-statistic: 1.14e+04 on 1 and 72 DF, p-value: < 2.2e-16Interpretación de los parámetros del modelo machos
Pendiente: 0.047292 indica que, por cada incremento de 1,000 gramos en la masa corporal, se estima que la longitud de la aleta aumente en promedio 47.292 milímetros. Esto representa la relación estimada entre ambas variables en el modelo sin intercepto.
Coeficiente de determinación: 0.9937 indica que el modelo explica aproximadamente el 99.37 % de la variabilidad observada en la longitud de la aleta a partir de la masa corporal.
Intercepto: El modelo no incluye un intercepto, ya que se ajustó mediante la fórmula
flipper_length_mm ~ body_mass_g + 0. Por lo tanto, se supone que la longitud de la aleta es igual a 0 mm cuando la masa corporal es de 0 gramos. Esta es una restricción matemática del modelo y no una interpretación biológica realista.
Nota: Aunque el \(R^2\) es muy alto, el error estándar residual es de 15.36 mm, por lo que conviene considerar también el comportamiento de los residuos al evaluar el ajuste del modelo.
Residuals:
Verificación de los supuestos
Supuesto de normalidad
\(H_0\) Hipótesis Nula: Los residuales del modelo siguen una distribución normal .
\(H_1\) : Los residuales del modelo no siguen una distribución normal.
lillie.test(residuals(modelo_machos))
Lilliefors (Kolmogorov-Smirnov) normality test
data: residuals(modelo_machos)
D = 0.054428, p-value = 0.8557Como \(p-valor=0.8557 >0.05\) no se rechaza \(H_0\) . No existen evidencias de que los residuales no se distribuyen de forma normal.
Supuesto de homocedasticidad
\(H_0\) : Los residuales son homocedásticos (la varianza de los errores es constante, ).
\(H_1\): Los residuales son heterocedásticos (la varianza de los errores no es constante).
bptest(modelo_machos,
~ body_mass_g + I(body_mass_g^2),
data = adelie_machos)
studentized Breusch-Pagan test
data: modelo_machos
BP = 48.547, df = 2, p-value = 2.871e-11Como \(p-valor=2.871e-11 < 0.05\) se rechaza \(H_0\) . Los residuales son heterocedásticos (la varianza de los errores no es constante).
Supuesto de Independencia
- \(H_0\) :: Los errores del modelo no presentan autocorrelación de primer orden positiva; se consideran independientes respecto a esa estructura de autocorrelación.
- \(H_1\) : Los errores del modelo presentan autocorrelación de primer orden positiva.
# 1. Modelo de machos
prueba_dw_machos <- dwtest(modelo_machos)
print(prueba_dw_machos)
Durbin-Watson test
data: modelo_machos
DW = 1.7398, p-value = 0.1314
alternative hypothesis: true autocorrelation is greater than 0Como \(p-valor=0.1314 > 0.05\) no se rechaza \(H_0\) . No se rechaza la hipotesis nula, es decir que no hay evidencia estadistica suficiente de autocorrelación positiva. Se cumple el supuesto de independencia
Transformación de BoxCox
# Transformación de Box-Cox
boxcox_result <- boxcox(modelo_machos,lambda = seq(-2, 2, by = 0.1))
Valor de lambda
lambda_optimo <- boxcox_result$x[
which.max(boxcox_result$y)
]
lambda_optimo[1] 0.9090909No se le aplica ninguna transformación sugerida dado a que 0.9090 es cercano a 1
Detección de Outliers
Gráfica de regresión lineal: machos Adelie
# Gráfica de regresión lineal: machos Adelie
ggplot(adelie_machos,
aes(x = body_mass_g, y = flipper_length_mm)) +
geom_point(color = "navy", size = 2) +
geom_smooth(
method = "lm",
formula = y ~ x + 0,
color = "red",
se = FALSE
) +
labs(
title = "Regresión Lineal: Peso vs Longitud de Aleta",
subtitle = "Modelo OLS sin intercepto para pingüinos Adelie machos",
x = "Peso corporal (g)",
y = "Longitud de aleta (mm)"
) +
theme_minimal()
Observaciones:
Tal como se observa en la gráfica, al anclar la línea roja en el origen —un punto muy por debajo de la nube real de observaciones—, la pendiente se distorsiona de forma artificial.
La recta se ve obligada a partir desde abajo, provocandol que subestime o sobrestime los valores reales dentro del rango observado, alejándose por completo de la tendencia natural que siguen los puntos.
Detección Numérica de Valores Atípicos para modelo machos
# Número de observaciones y parámetros
n <- nobs(modelo_machos)
p <- 1 # Solo pendiente, sin intercepto
# 1. Leverage
hat_vals_machos <- hatvalues(modelo_machos)
umbral_hat_machos <- 3 * (p / n)
# 2. Distancia de Cook
cook_vals_machos <- cooks.distance(modelo_machos)
umbral_cook_machos <- 4 / (n - p)
# 3. Residuos estudentizados
r_student_machos <- rstudent(modelo_machos)
# Identificar observaciones problemáticas
atipicos_machos <- which(abs(r_student_machos) > 3)
leverage_machos <- which(hat_vals_machos > umbral_hat_machos)
cook_machos <- which(cook_vals_machos > umbral_cook_machos)
cat("Observaciones con RStudent > 3:", atipicos_machos, "\n")Observaciones con RStudent > 3: cat("Observaciones con Leverage alto:", leverage_machos, "\n")Observaciones con Leverage alto: cat("Observaciones con Distancia de Cook alta:", cook_machos, "\n")Observaciones con Distancia de Cook alta: 17 20 38 52 plot(modelo_machos, which = 4,
col = "blue", pch = 16,
main = "Distancia de Cook - Adelie machos")
abline(h = umbral_cook_machos,
col = "red", lty = 2)
cook_machos <- cooks.distance(modelo_machos)
indices_machos <- which( cook_machos > 4/(nrow(adelie_machos) - 2))
data.frame(
Observacion = indices_machos,
Peso = adelie_machos$body_mass_g[indices_machos],
Longitud_aleta = adelie_machos$flipper_length_mm[indices_machos],
Cook = cook_machos[indices_machos]) Observacion Peso Longitud_aleta Cook
17 17 4650 184 0.10191356
20 20 4600 190 0.05863423
52 52 4775 197 0.06936300No es posible determinar si estas observaciones corresponden a errores de medición o a variaciones biológicas extraordinarias ajenas al comportamiento general de los datos. Por tanto, no se justifica su eliminación y se propone evaluar la winsorización para reajustar el modelo y analizar su efecto sobre las estimaciones.
Winsorización
# 1. Definir límites de percentiles
q_inf_machos <- quantile( adelie_machos$flipper_length_mm, 0.05, na.rm = TRUE)
q_sup_machos <- quantile( adelie_machos$flipper_length_mm, 0.95, na.rm = TRUE)
# 2. Winsorizar longitud de aleta
adelie_machos$flipper_winsor <- pmin( pmax(adelie_machos$flipper_length_mm,q_inf_machos), q_sup_machos)
# 3. Ajustar el modelo winsorizado
modelo_machos_winsor <- lm( flipper_winsor ~ body_mass_g+0, data = adelie_machos)
# 4. Comparar modelos
summary(modelo_machos)
Call:
lm(formula = flipper_length_mm ~ body_mass_g + 0, data = adelie_machos)
Residuals:
Min 1Q Median 3Q Max
-35.908 -8.355 1.645 12.014 32.478
Coefficients:
Estimate Std. Error t value Pr(>|t|)
body_mass_g 0.047292 0.000443 106.7 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.36 on 72 degrees of freedom
Multiple R-squared: 0.9937, Adjusted R-squared: 0.9936
F-statistic: 1.14e+04 on 1 and 72 DF, p-value: < 2.2e-16summary(modelo_machos_winsor)
Call:
lm(formula = flipper_winsor ~ body_mass_g + 0, data = adelie_machos)
Residuals:
Min 1Q Median 3Q Max
-35.708 -8.171 0.916 12.185 32.628
Coefficients:
Estimate Std. Error t value Pr(>|t|)
body_mass_g 0.0472491 0.0004394 107.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.23 on 72 degrees of freedom
Multiple R-squared: 0.9938, Adjusted R-squared: 0.9937
F-statistic: 1.157e+04 on 1 and 72 DF, p-value: < 2.2e-16Comparando los resultados de los parametros del modelo, se tiene una ligera mejoría con el modelo de machos winsorizado, sin embargo permanece alto el error estandar residual, por lo que se duda la pertinencia del modelo.
valores_winsorizados <- which(
adelie_machos$flipper_length_mm !=
adelie_machos$flipper_winsor
)
cat("Observaciones winsorizadas:",
valores_winsorizados, "\n")Observaciones winsorizadas: 8 10 12 13 43 45 48 62 plot(
adelie_machos$body_mass_g,
adelie_machos$flipper_length_mm,
pch = 16,
col = "gray70",
main = "Efecto de la Winsorización en la Regresión Lineal",
xlab = "Masa corporal (g)",
ylab = "Longitud de aleta (mm)"
)
# Puntos originales que fueron modificados
points(
adelie_machos$body_mass_g[valores_winsorizados],
adelie_machos$flipper_length_mm[valores_winsorizados],
col = "purple",
pch = 17,
cex = 1.3
)
# Recta del modelo original (sin intercepto)
abline(
a = 0,
b = coef(modelo_machos)[1],
col = "red",
lwd = 2
)
# Recta del modelo winsorizado (sin intercepto)
abline(
a = 0,
b = coef(modelo_machos_winsor)[1],
col = "forestgreen",
lwd = 2,
lty = 2
)
# Leyenda
legend(
"topleft",
legend = c(
"OLS Original",
"OLS Winsorizado",
"Valores modificados"
),
col = c("red", "forestgreen", "purple"),
lty = c(1, 2, NA),
pch = c(NA, NA, 17),
lwd = 2
)
Supuesto de homocedasticidad modelo original
\(H_0\) : Los residuales son homocedásticos (la varianza de los errores es constante, ).
\(H_1\): Los residuales son heterocedásticos (la varianza de los errores no es constante).
bptest(modelo_machos, ~ body_mass_g + I(body_mass_g^2), data = adelie_machos)
studentized Breusch-Pagan test
data: modelo_machos
BP = 48.547, df = 2, p-value = 2.871e-11Como \(p-valor=2.871e-11 < 0.05\) se rechaza \(H_0\) . Los residuales son heterocedásticos (la varianza de los errores no es constante).
Supuesto de homocedasticidad modelo winsorizado
\(H_0\) : Los residuales son homocedásticos (la varianza de los errores es constante, ).
\(H_1\): Los residuales son heterocedásticos (la varianza de los errores no es constante).
bptest(modelo_machos_winsor, ~ body_mass_g + I(body_mass_g^2), data = adelie_machos)
studentized Breusch-Pagan test
data: modelo_machos_winsor
BP = 50.234, df = 2, p-value = 1.235e-11Como \(p-valor=2.871e-11 < 0.05\) se rechaza \(H_0\) . Los residuales son heterocedásticos (la varianza de los errores no es constante).
Conclusión:
Una vez ajustado el modelo con intercepto, los residuos continuaron mostrando heterocedasticidad,por lo tanto se debe volver a evaluar Box-Cox o considerar un modelo que estime un intercepto (\(\beta_0 \neq 0\)).