Practica 3
Ruth Sua Cruz Mendez
2025-12-14
Introducción
Seleccionar una de las variables de X14-X21, para ajustar un modelo respecto de la estatura.
Ajustaré el modelo en función de la longitud de mano.
Estatura total (X11)∼Longitud de la pierna (X20)
##Abriendo paquete pacman
library(pacman)## Warning: package 'pacman' was built under R version 4.5.2## Abriendo las librerias o paquetes que se usan
p_load(haven,dplyr,ggplot2,tinytex,tidyr,GGally,purrr,labelled,tidyverse,FactoMineR,factoextra)
Cholula <- read_sav("Cholula.sav")Procedencia de los datos
table(Cholula$ORIGEN)##
## Adulta diversa procedencia Adulta Valle de Cholula
## 22 37
## Juvenil Cholula Juvenil diversa procedencia
## 165 58
## Juvenil San Nicolás de los Rancho Juvenil Santa Isabel Cholula
## 29 28Seleccionar variables
Seleccionamos una variable correspondiente a una parte del brazo Longitud de la mano
datos <- Cholula %>% dplyr::select(X11, X17)Extraer etiquetas SPSS
lab_X11 <- attr(datos$X11, "label")
lab_X17 <- attr(datos$X17, "label")Ajuste del modelo de minimos cuadrados
modelo <- lm(X11~X17, data = datos)
summary(modelo)##
## Call:
## lm(formula = X11 ~ X17, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -147.594 -31.861 -3.309 23.847 181.423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 991.0258 50.3650 19.68 <2e-16 ***
## X17 3.2849 0.2663 12.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.3 on 332 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3143, Adjusted R-squared: 0.3122
## F-statistic: 152.2 on 1 and 332 DF, p-value: < 2.2e-16anova(modelo)## Analysis of Variance Table
##
## Response: X11
## Df Sum Sq Mean Sq F value Pr(>F)
## X17 1 369763 369763 152.16 < 2.2e-16 ***
## Residuals 332 806795 2430
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coef <- summary(modelo)$coefficients
round(coef,4)## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 991.0258 50.3650 19.6769 0
## X17 3.2849 0.2663 12.3353 0sqrt(2430) #Raíz cuadrada de la varianza de los errores## [1] 49.29503op <- par() # guarda la configuración actual del panel de gráfica
par(mar = c(4, 4, 2, 1)) # reduce márgenes
par(mfrow = c(2, 2))
plot(modelo)
Residual vs Fitted
Denota que las observaciones 269 y 273, asi como el 73 en la inferior de la linea podrian sesgar la linealidad de la grafica
Se observa heterocedasticidad con valores muy dispersos en los extremos (colas)
Q-Q Plot
Se observan desviaciones notables con colas pesadas(las observaciones que se salen de la linea)
Además de que volvemos a encontrarnos con que las observaciones: 269, 273, 73 continuan siendo los que presentan una discordancia importante con el resto de las observaciones
Los extremos tambien muestran valores atipicos y los residuos centrales son aproximadamente normales por lo que tambien presenta una ligera heterocedasticidad
Scale-Location
La linea tiene una ligera pendiente descendente así como persistencia de heterocedasticidad moderada ya que la variabilidad de los residuos aumenta conforme aumenta el valor ajustado
Leverage
No hay muchas observaciones que afecten el modelo, pero si datos con influencia moderada, que me entregan residuos estandarizados grandes, afectando la pendiente del modelo.
Por lo que podemos concluir que este modelo es razonablemente util aunque presenta ciertas consideraciones:
- Ajuste lineal aceptable
- Normalidad aceptable salvo colas
- Heterocedasticidad moderada
- Valores atipicos y puntos influyentes
Nota: Hay que evitar la eliminacion de los valores atipicos sin una revisión previa, ya que pueden sesgar el modelo hacia valores promedio
En el caso de la observación 73, esta en la gráfica Residual vs Fitted y la Scale-Location donde en uno aparece por debajo de la linea y en la otra por arriba igual que la observación 269y 273, esto se debe a que en la primera se coloca con respecto a sus valores mientras que en la segunda se ubica considerando su raíz indicando que este valor es particularmente influyente
Consecuente en este caso debemos evaluar que opciones tenemos para perfeccionar nuestros supuestos.
Modelo lineal simple vs Modelo robusto
El modelo \(Estatura= \beta_0 + \beta_1(Longitud de mano)\) es util pero es deficiente debido a las consideraciones antes mencionadas. provocando que sus supuestos no se cumplen a la perfección, lo cual un Modelo robusto (MASS::rlm) compensa reduciendo la influencia de los valores atipicos sin eliminarlos, tolera la heterocedasticidad, sus estimaciones son más estables y mejora la precisión para casos reales.
En ciencias forenses el modelo robusto es particularmente útil porque las observaciones sulen contener muestras con mediciones imperfectas o erroneas.*
Extraemos los datos usados por el modelo
# Extraer solo datos usados por el modelo
datos_modelo <- model.frame(modelo)Diagnóstico del modelo
Para determinar cual seria el modelo más optimo realizare una serie de diagnosticos para evaluar los modelos.
Diagnósticos
diagnosticos <- data.frame(
estatura = datos_modelo$X11,
mano = datos_modelo$X17,
residuales = residuals(modelo),
residuales_est = rstandard(modelo),
leverage = hatvalues(modelo),
cooks = cooks.distance(modelo)
)Identificar posibles outliers por residuales estandarizados
outliers_residuales <- diagnosticos[abs(diagnosticos$residuales_est) > 2, ]Identificar puntos con alto leverage
Regla: leverage > 2p/n (p = nº de parámetros del modelo)
p <- length(coefficients(modelo))
n <- nrow(datos)
umbral_leverage <- 2 * p / n
outliers_leverage <- diagnosticos[diagnosticos$leverage > umbral_leverage, ]Identificar puntos influyentes por distancia de Cook
Regla común: Cook > 4/n
umbral_cook <- 4 / n
outliers_cook <- diagnosticos[diagnosticos$cooks > umbral_cook, ]Mostrar resultados
OUTLIERS POR RESIDUALES ESTANDARIZADOS
print(outliers_residuales)## estatura mano residuales residuales_est leverage cooks
## 73 1438 181 -147.5940 -3.001235 0.004795193 0.021700190
## 80 1487 195 -144.5827 -2.938969 0.004095512 0.017760306
## 96 1488 186 -114.0186 -2.316679 0.003232094 0.008701454
## 132 1501 193 -124.0129 -2.520082 0.003495086 0.011137241
## 224 1742 187 136.6965 2.777269 0.003094569 0.011971598
## 234 1740 189 128.1267 2.603025 0.002994615 0.010175839
## 238 1661 161 141.1041 2.899794 0.025638844 0.110632518
## 250 1741 193 115.9871 2.356988 0.003495086 0.009742337
## 254 1778 204 116.8531 2.381998 0.009686496 0.027748959
## 269 1793 197 154.8475 3.148941 0.004929398 0.024560606
## 273 1836 202 181.4229 3.695145 0.008035500 0.055303143
## 277 1717 187 111.6965 2.269343 0.003094569 0.007993123
## 278 1788 201 136.7078 2.783372 0.007297549 0.028475443
## 280 1642 164 112.2494 2.301365 0.021023988 0.056870095
## 318 1777 206 109.2833 2.229813 0.011570952 0.029102513
## 335 1699 180 116.6909 2.373420 0.005282908 0.014958657
## 337 1720 180 137.6909 2.800547 0.005282908 0.020827115
## 339 1774 205 109.5682 2.234529 0.010599542 0.026745873PUNTOS CON ALTO LEVERAGE
print(outliers_leverage)## estatura mano residuales residuales_est leverage cooks
## 4 1691 207 19.9983860 0.408259190 0.01260073 1.063518e-03
## 7 1701 209 23.4285712 0.478827210 0.01483537 1.726304e-03
## 17 1710 214 16.0040340 0.328188915 0.02144338 1.180117e-03
## 21 1673 214 -20.9959660 -0.430556652 0.02144338 2.031131e-03
## 30 1740 207 68.9983860 1.408574928 0.01260073 1.265997e-02
## 58 1711 216 10.4342191 0.214305174 0.02449513 5.766147e-04
## 117 1649 210 -31.8563363 -0.651470241 0.01604024 3.459337e-03
## 167 1536 171 -16.7449463 -0.341789238 0.01229877 7.273154e-04
## 200 1517 170 -32.4600388 -0.662917264 0.01337013 2.977626e-03
## 202 1536 166 -0.3204091 -0.006559782 0.01823924 3.997146e-07
## 204 1553 170 3.5399612 0.072295088 0.01337013 3.541352e-05
## 223 1604 168 61.1097760 1.249486925 0.01568796 1.244134e-02
## 228 1602 165 68.9644983 1.412901519 0.01960243 1.995729e-02
## 230 1602 170 52.5399612 1.073000789 0.01337013 7.801023e-03
## 238 1661 161 141.1041281 2.899794097 0.02563884 1.106325e-01
## 242 1702 216 1.4342191 0.029456980 0.02449513 1.089424e-05
## 243 1649 215 -48.2808734 -0.990836372 0.02294007 1.152517e-02
## 280 1642 164 112.2494058 2.301364707 0.02102399 5.687009e-02
## 281 1739 212 51.5738489 1.056086495 0.01862508 1.058357e-02
## 282 1668 210 -12.8563363 -0.262915371 0.01604024 5.634243e-04
## 284 1704 209 26.4285712 0.540140451 0.01483537 2.196712e-03
## 290 1483 171 -69.7449463 -1.423598000 0.01229877 1.261772e-02
## 293 1637 208 -37.2865214 -0.761609450 0.01368887 4.025207e-03
## 306 1546 164 16.2494058 0.333149282 0.02102399 1.191766e-03
## 309 1632 171 79.2550537 1.617713426 0.01229877 1.629330e-02
## 314 1518 171 -34.7449463 -0.709195987 0.01229877 3.131400e-03PUNTOS INFLUYENTES (COOK’S DISTANCE)
print(outliers_cook)## estatura mano residuales residuales_est leverage cooks
## 30 1740 207 68.99839 1.408575 0.012600728 0.01265997
## 50 1555 200 -93.00726 -1.892980 0.006617964 0.01193632
## 73 1438 181 -147.59402 -3.001235 0.004795193 0.02170019
## 80 1487 195 -144.58272 -2.938969 0.004095512 0.01776031
## 223 1604 168 61.10978 1.249487 0.015687958 0.01244134
## 224 1742 187 136.69653 2.777269 0.003094569 0.01197160
## 228 1602 165 68.96450 1.412902 0.019602433 0.01995729
## 238 1661 161 141.10413 2.899794 0.025638844 0.11063252
## 254 1778 204 116.85311 2.381998 0.009686496 0.02774896
## 269 1793 197 154.84746 3.148941 0.004929398 0.02456061
## 273 1836 202 181.42292 3.695145 0.008035500 0.05530314
## 278 1788 201 136.70783 2.783372 0.007297549 0.02847544
## 280 1642 164 112.24941 2.301365 0.021023988 0.05687009
## 290 1483 171 -69.74495 -1.423598 0.012298767 0.01261772
## 309 1632 171 79.25505 1.617713 0.012298767 0.01629330
## 318 1777 206 109.28329 2.229813 0.011570952 0.02910251
## 329 1667 177 94.54561 1.924754 0.007096243 0.01323859
## 335 1699 180 116.69089 2.373420 0.005282908 0.01495866
## 337 1720 180 137.69089 2.800547 0.005282908 0.02082712
## 339 1774 205 109.56820 2.234529 0.010599542 0.02674587Ordenar por Cook’s distance (más influyentes primero)
diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ]
head(diagnosticos_ordenados, 10)## estatura mano residuales residuales_est leverage cooks
## 238 1661 161 141.1041 2.899794 0.025638844 0.11063252
## 280 1642 164 112.2494 2.301365 0.021023988 0.05687009
## 273 1836 202 181.4229 3.695145 0.008035500 0.05530314
## 318 1777 206 109.2833 2.229813 0.011570952 0.02910251
## 278 1788 201 136.7078 2.783372 0.007297549 0.02847544
## 254 1778 204 116.8531 2.381998 0.009686496 0.02774896
## 339 1774 205 109.5682 2.234529 0.010599542 0.02674587
## 269 1793 197 154.8475 3.148941 0.004929398 0.02456061
## 73 1438 181 -147.5940 -3.001235 0.004795193 0.02170019
## 337 1720 180 137.6909 2.800547 0.005282908 0.02082712Ajustar el modelo robusto
library(MASS)##
## Adjuntando el paquete: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectmodelo_robusto <- rlm(X11 ~ X17, data = datos)
summary(modelo_robusto)##
## Call: rlm(formula = X11 ~ X17, data = datos)
## Residuals:
## Min 1Q Median 3Q Max
## -145.299 -29.310 -1.025 26.314 183.944
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 990.6794 48.7312 20.3295
## X17 3.2741 0.2577 12.7071
##
## Residual standard error: 40.96 on 332 degrees of freedom
## (5 observations deleted due to missingness)Comparación de los modelos
coef(modelo)## (Intercept) X17
## 991.025775 3.284907coef(modelo_robusto)## (Intercept) X17
## 990.679408 3.274143Datos del modelo robusto
coef_rob <- coef(modelo_robusto)Datos del modelo OLS
coef_ols <- coef(modelo)Gráficamos
ggplot(datos, aes(x = X17, y = X11)) +
geom_point(alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, color = "gray40", linetype = "dashed") +
geom_smooth(method = MASS::rlm, se = FALSE, color = "red") +
labs(
title = "Modelo robusto (RLM) vs modelo lineal (OLS)",
x = "Longitud de la pierna",
y = "Estatura"
) +
theme_minimal(base_size = 14)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
Función de pseudo-R² robusto
pseudoR2 <- function(modelo) {
y <- modelo$fitted + modelo$residuals
mad_y <- mad(y) # variabilidad total
mad_res <- mad(modelo$residuals) # variabilidad residual
1 - (mad_res^2 / mad_y^2)
}
pseudoR2(modelo_robusto)## [1] 0.4904967Comparativa OLS vs Robusto
Coeficientes y SE de OLS
ols_sum <- summary(modelo)
ols_tab <- data.frame(
Modelo = "OLS",
Intercepto = ols_sum$coefficients[1,1],
SE_Intercepto = ols_sum$coefficients[1,2],
Pendiente = ols_sum$coefficients[2,1],
SE_Pendiente = ols_sum$coefficients[2,2]
)Coeficientes y SE de RLM
rlm_sum <- summary(modelo_robusto)
rlm_tab <- data.frame(
Modelo = "RLM",
Intercepto = rlm_sum$coefficients[1,1],
SE_Intercepto = rlm_sum$coefficients[1,2],
Pendiente = rlm_sum$coefficients[2,1],
SE_Pendiente = rlm_sum$coefficients[2,2]
)
R2_ols <- summary(modelo)$r.squared
R2_rlm <- pseudoR2(modelo_robusto)
AIC_rlm <- AIC(lm(modelo_robusto$wresid ~ modelo_robusto$fitted))
AIC_ols <- AIC(modelo)
comparacion <- data.frame(
Modelo = c("OLS", "RLM"),
Intercepto = c(ols_tab$Intercepto, rlm_tab$Intercepto),
SE_Intercepto = c(ols_tab$SE_Intercepto, rlm_tab$SE_Intercepto),
Pendiente = c(ols_tab$Pendiente, rlm_tab$Pendiente),
SE_Pendiente = c(ols_tab$SE_Pendiente, rlm_tab$SE_Pendiente),
R2 = c(R2_ols, R2_rlm),
AIC = c(AIC_ols, AIC_rlm)
)
comparacion## Modelo Intercepto SE_Intercepto Pendiente SE_Pendiente R2 AIC
## 1 OLS 991.0258 50.36503 3.284907 0.2663017 0.3142753 3555.605
## 2 RLM 990.6794 48.73121 3.274143 0.2576630 0.4904967 3555.605Conclusión
Los valores de ambos modelos no se diferencian significativamente OLS vs Robusto por lo que el uso de uno u otro no genera diferencia significativa a pesar de la situaciones que se presenta en el OLS