Practica 3
Nohemi Angel
2025-12-09
INTRODUCCION
En esta práctica se presenta el análisis realizado en la sesión del 28 de noviembre de la materia Métodos Cuantitativos, impartida por el profesor José Luis Castrejón. El ejercicio consistió en explorar la relación entre la estatura y distintas variables corporales, así como en ajustar e interpretar modelos de regresión lineal y robusta.
setwd("~/Blanca")
##Abriendo paquete pacman
library(pacman)
## El archivo esta en formato SPSS, lo abrimos mediante la libreria haven
p_load(haven,dplyr,ggplot2,tinytex,tidyr,GGally,purrr,labelled)
### abriendo los paquetes o las librerias que se usan
Cholula <- read_sav("Cholula.sav")Se selecciono varias variables que se comprenden desde la x11, x14 a la x21
vars <- Cholula %>% dplyr::select(X11, X14:X21)Se extraen etiquetas spss
# 2. Extraer etiquetas SPSS
etiquetas_raw <- map_chr(names(vars), ~{
lab <- attr(vars[[.x]], "label")
if (is.null(lab)) .x else lab
})etiquetas <- make.unique(etiquetas_raw, sep = "_")Se renombran etiquetas
df <- vars %>% rename_with(~etiquetas)Se generan graficas por variable de pares con respecto a la estatura
graficar_cor <- function(y_var){
datos <- df %>% dplyr::select(`Estatura total`, all_of(y_var))
prueba <- cor.test(datos[[1]], datos[[2]], use = "pairwise.complete.obs")
r <- round(prueba$estimate, 3)
p <- formatC(prueba$p.value, digits = 3, format = "e")
ggplot(datos, aes(x = `Estatura total`, y = .data[[y_var]])) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE) +
labs(
title = paste0("Estatura vs ", y_var),
subtitle = paste0("r = ", r, " | p = ", p),
x = "Estatura total",
y = y_var
) +
theme_minimal(base_size = 14)
}Ahora se generan todas las graficas de la variable x11 a la x19
graficas <- map(names(df)[-1], graficar_cor)
graficas ## [[1]]## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## [[2]]## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## [[3]]## `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()`).
##
## [[4]]## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## [[5]]## `geom_smooth()` using formula = 'y ~ x'
##
## [[6]]## `geom_smooth()` using formula = 'y ~ x'
##
## [[7]]## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
##
## [[8]]## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
En estas graficas podemos observar que entre mas alejada la medicion de las extremidades hay una menor relacion a la estatura, esto sugiere que las variables más próximas anatómicamente a la estatura presentan relaciones más fuertes, mientras que las más distales muestran mayor dispersión, lo que es consistente con la variabilidad morfológica del cuerpo humano.
Se selecciono la variable longitud del muslo para ajustar el modelo en funcion con respecto a la estatura.
## Estatura total (X11)∼Longitud del muslo (X19)
datos <- Cholula %>% dplyr::select(X11, X19)
# Extraer etiquetas SPSS
lab_X11 <- attr(datos$X11, "label")
lab_X19 <- attr(datos$X19, "label")
## Aquí estoy ajustando el modelo de mínimos cuadrados o clásico
modelo <- lm(X11 ~ X19, data = datos)
summary(modelo)##
## Call:
## lm(formula = X11 ~ X19, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -100.053 -25.388 -1.753 20.212 136.853
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 771.80913 37.22876 20.73 <2e-16 ***
## X19 1.76676 0.07823 22.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 37.47 on 337 degrees of freedom
## Multiple R-squared: 0.6021, Adjusted R-squared: 0.6009
## F-statistic: 510 on 1 and 337 DF, p-value: < 2.2e-16El modelo dice que los residuales varian entre -100.053 y 136.853, mientras que el Beta o es 771.80913 y el x19 que es la longitud del muslo es de 1.76676, esto indica que la pendiente en promedio, por cada milímetro adicional en la longitud del muslo, la estatura total aumenta aproximadamente 1.77 mm, lo que confirma una relación positiva y fuerte entre ambas variables.
anova(modelo)## Analysis of Variance Table
##
## Response: X11
## Df Sum Sq Mean Sq F value Pr(>F)
## X19 1 715892 715892 510 < 2.2e-16 ***
## Residuals 337 473054 1404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El presente modelo nos da los grados de residuales, donde la raiz del promedio de los residuales y los observados es 1404, es decir, el análisis muestra que la longitud del muslo tiene una proporción significativa de la variabilidad en la estatura, lo que confirma la relevancia del predictor dentro del modelo.
Raiz cuadrada de la varianza de los errores
sqrt(1404) ## [1] 37.46999op <- 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)
Ahora sacaremos los puntos que son atipicos por residuos estandarizados,apalancamiento, distancia de Cook, exclusion.
# Extraer solo datos usados por el modelo
datos_modelo <- model.frame(modelo)
# Diagnósticos
diagnosticos <- data.frame(
estatura = datos_modelo$X11,
muslo = datos_modelo$X19,
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, ]
# 6. 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 ---
cat("---- OUTLIERS POR RESIDUALES ESTANDARIZADOS ----\n")## ---- OUTLIERS POR RESIDUALES ESTANDARIZADOS ----print(outliers_residuales)## estatura muslo residuales residuales_est leverage cooks
## 2 1701 460 116.48040 3.115099 0.003951070 0.019246326
## 56 1588 508 -81.32417 -2.178951 0.007654096 0.018310313
## 80 1487 448 -76.31846 -2.043298 0.006164675 0.012948790
## 96 1488 462 -100.05312 -2.675444 0.003704222 0.013306704
## 105 1530 479 -88.08808 -2.354681 0.003014368 0.008381880
## 254 1778 523 82.17440 2.207610 0.012931704 0.031924431
## 273 1836 556 81.87126 2.220392 0.031449034 0.080041279
## 308 1603 393 136.85345 3.713318 0.032377787 0.230693693
## 310 1639 425 116.31707 3.126412 0.013917397 0.068977435
## 315 1735 495 88.64373 2.371501 0.004667292 0.013186004
## 316 1696 458 115.01392 3.076316 0.004232800 0.020114149
## 320 1567 493 -75.82274 -2.028164 0.004338590 0.008962165
## 335 1699 461 112.71364 3.014169 0.003823286 0.017434347
## 339 1774 514 94.07526 2.522981 0.009529687 0.030622123cat("\n---- PUNTOS CON ALTO LEVERAGE ----\n")##
## ---- PUNTOS CON ALTO LEVERAGE ----print(outliers_leverage)## estatura muslo residuales residuales_est leverage cooks
## 9 1734 523 38.174402 1.02555274 0.01293170 6.889608e-03
## 73 1438 338 69.025351 1.92597383 0.08497025 1.722275e-01
## 119 1525 424 4.083829 0.10979120 0.01435912 8.780396e-05
## 132 1501 423 -18.149409 -0.48804713 0.01480956 1.790257e-03
## 152 1493 426 -31.449695 -0.84513059 0.01348440 4.881409e-03
## 184 1670 524 -27.592360 -0.74142527 0.01335331 3.719906e-03
## 200 1517 419 4.917638 0.13236482 0.01669853 1.487670e-04
## 201 1497 427 -29.216457 -0.78494922 0.01306012 4.076707e-03
## 203 1552 422 34.617353 0.93109575 0.01526872 6.721150e-03
## 214 1507 426 -17.449695 -0.46891619 0.01348440 1.502755e-03
## 232 1695 531 -14.959694 -0.40262954 0.01654871 1.363931e-03
## 234 1740 536 21.206497 0.57149749 0.01909275 3.178625e-03
## 254 1778 523 82.174402 2.20760979 0.01293170 3.192443e-02
## 268 1736 530 27.807068 0.74822393 0.01606606 4.570636e-03
## 269 1793 539 68.906211 1.85851101 0.02072382 3.654811e-02
## 273 1836 556 81.871259 2.22039169 0.03144903 8.004128e-02
## 277 1717 530 8.807068 0.23697785 0.01606606 4.584891e-04
## 278 1788 537 67.439735 1.81794067 0.01962772 3.308326e-02
## 305 1498 410 1.818495 0.04906615 0.02145885 2.639741e-05
## 306 1546 425 23.317067 0.62672444 0.01391740 2.771839e-03
## 308 1603 393 136.853447 3.71331839 0.03237779 2.306937e-01
## 310 1639 425 116.317067 3.12641162 0.01391740 6.897743e-02
## 313 1638 531 -71.959694 -1.93674409 0.01654871 3.155918e-02
## 318 1777 545 42.305640 1.14309412 0.02422142 1.621744e-02
## 324 1647 522 -47.058837 -1.26396816 0.01251882 1.012691e-02
## 337 1720 535 2.973259 0.08010537 0.01856650 6.069634e-05cat("\n---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----\n")##
## ---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----print(outliers_cook)## estatura muslo residuales residuales_est leverage cooks
## 2 1701 460 116.48040 3.115099 0.003951070 0.01924633
## 5 1611 513 -67.15798 -1.800790 0.009195288 0.01504781
## 30 1740 507 72.44259 1.940708 0.007372020 0.01398589
## 56 1588 508 -81.32417 -2.178951 0.007654096 0.01831031
## 73 1438 338 69.02535 1.925974 0.084970254 0.17222749
## 80 1487 448 -76.31846 -2.043298 0.006164675 0.01294879
## 96 1488 462 -100.05312 -2.675444 0.003704222 0.01330670
## 224 1742 510 69.14231 1.853109 0.008244412 0.01427338
## 250 1741 516 57.54173 1.543739 0.010224648 0.01230920
## 254 1778 523 82.17440 2.207610 0.012931704 0.03192443
## 269 1793 539 68.90621 1.858511 0.020723822 0.03654811
## 273 1836 556 81.87126 2.220392 0.031449034 0.08004128
## 278 1788 537 67.43974 1.817941 0.019627721 0.03308326
## 308 1603 393 136.85345 3.713318 0.032377787 0.23069369
## 310 1639 425 116.31707 3.126412 0.013917397 0.06897743
## 313 1638 531 -71.95969 -1.936744 0.016548709 0.03155918
## 315 1735 495 88.64373 2.371501 0.004667292 0.01318600
## 316 1696 458 115.01392 3.076316 0.004232800 0.02011415
## 318 1777 545 42.30564 1.143094 0.024221416 0.01621744
## 335 1699 461 112.71364 3.014169 0.003823286 0.01743435
## 339 1774 514 94.07526 2.522981 0.009529687 0.03062212# Ordenar por Cook’s distance (más influyentes primero) ---
diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ]
head(diagnosticos_ordenados, 10)## estatura muslo residuales residuales_est leverage cooks
## 308 1603 393 136.85345 3.713318 0.032377787 0.23069369
## 73 1438 338 69.02535 1.925974 0.084970254 0.17222749
## 273 1836 556 81.87126 2.220392 0.031449034 0.08004128
## 310 1639 425 116.31707 3.126412 0.013917397 0.06897743
## 269 1793 539 68.90621 1.858511 0.020723822 0.03654811
## 278 1788 537 67.43974 1.817941 0.019627721 0.03308326
## 254 1778 523 82.17440 2.207610 0.012931704 0.03192443
## 313 1638 531 -71.95969 -1.936744 0.016548709 0.03155918
## 339 1774 514 94.07526 2.522981 0.009529687 0.03062212
## 316 1696 458 115.01392 3.076316 0.004232800 0.02011415En este caso se puede observar que en el caso de los residuales estandarizados tinen un valor en residuales de 116.48040 y ya estandarizados de 3.115099 esto en el caso 2 pero tambien en en caso 56, 80, 96, 105,254, 273, 308, 310, 315,316, 320, 335, 339.
En el caso de Leverage, se identifican varias observaciones con valores relativamente altos, como los casos 9, 73 y 119. No obstante, estos puntos no presentan simultáneamente residuales extremos ni valores críticos de distancia de Cook, por lo que no comprometen la estabilidad del modelo.
Respecto a la distancia de Cook, se identifican algunas observaciones con valores relativamente elevados; sin embargo, ninguno supera umbrales críticos que indiquen una influencia desproporcionada sobre el ajuste del modelo.
Esto nos permite identificar observaciones con residuales estandarizados elevados, distancia de leverage moderados y algunos puntos con mayor distancia de Cook.
Estos resultados indican que el modelo cumple adecuadamente los supuestos de la regresión lineal y se mantiene estable y robusto, por lo que las estimaciones obtenidas pueden considerarse confiables.
Ajustar el modelo robusto
library(MASS)##
## Adjuntando el paquete: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectmodelo_robusto <- rlm(X11 ~ X19, data = datos)
summary(modelo_robusto)##
## Call: rlm(formula = X11 ~ X19, data = datos)
## Residuals:
## Min 1Q Median 3Q Max
## -97.3417 -22.8522 -0.3585 20.7392 144.2066
##
## Coefficients:
## Value Std. Error t value
## (Intercept) 738.0185 36.4919 20.2242
## X19 1.8340 0.0767 23.9163
##
## Residual standard error: 33.12 on 337 degrees of freedom## Comparación de los modelos
coef(modelo)## (Intercept) X19
## 771.809130 1.766762coef(modelo_robusto)## (Intercept) X19
## 738.018470 1.834033El modelo original es de 771.809130 y bajo a 738.018470, mientras que la variable de 1.766762 subio a 1.834033, lo que indica que la estimación del modelo es estable y poco sensible a la presencia de valores atípicos.
# Datos del modelo robusto
coef_rob <- coef(modelo_robusto)
# Datos del modelo OLS
coef_ols <- coef(modelo)
ggplot(datos, aes(x = X19, 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 del muslo",
y = "Estatura"
) +
theme_minimal(base_size = 14)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Los modelos siguieren que los valores atípicos detectados no distorsionan de manera importante el modelo y que la relación estimada es estable.
# 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.6919875Vamos a comparar los dos modelos el clásico y el 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 771.8091 37.22876 1.766762 0.07823391 0.6021228 3422.728
## 2 RLM 738.0185 36.49192 1.834033 0.07668548 0.6919875 3422.728resid_est <- rstandard(modelo)
lev <- hatvalues(modelo)
cook <- cooks.distance(modelo)
umbral_cook <- 4 / n
umbral_lev <- 2 * p / n
indices_fuera <- which(abs(resid_est) > 2 |
lev > umbral_lev |
cook > umbral_cook)
length(indices_fuera) # cuántos eliminarás## [1] 40datos_limpios <- datos[-indices_fuera, ]dando como resultado solo 299 variables, el modelo quita 40 datos.
modelo_limpio <- lm(X11 ~ X19, data = datos_limpios)
summary(modelo_limpio)##
## Call:
## lm(formula = X11 ~ X19, data = datos_limpios)
##
## Residuals:
## Min 1Q Median 3Q Max
## -71.167 -22.289 -0.597 16.837 74.384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 760.41631 42.80511 17.77 <2e-16 ***
## X19 1.78505 0.09016 19.80 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.91 on 297 degrees of freedom
## Multiple R-squared: 0.5689, Adjusted R-squared: 0.5675
## F-statistic: 392 on 1 and 297 DF, p-value: < 2.2e-16La eliminación de observaciones atípicas no mejora el poder explicativo del modelo y reduce la variabilidad de la muestra, por lo que no resulta metodológicamente conveniente excluir dichos casos.
Comparacion de los modelos
## Comparando los modelos
summary(modelo)$r.squared # R² original## [1] 0.6021228summary(modelo_limpio)$r.squared # R² sin atípicos## [1] 0.5689381Este análisis confirma una relación positiva, significativa y estable entre la estatura total y la longitud del muslo. Los diagnósticos del modelo indican que, aunque existen valores atípicos y observaciones influyentes, estos no distorsionan de manera sustancial el ajuste. La comparación entre el modelo clásico, el robusto y el modelo limpio demuestra que la estimación es consistente.