El objetivo es aplicar modelos de regresión lineal para

estimar la estatura a partir de longitud (es) de las extremidades

Defino mi directorio de trabajo

setwd(“C:/Users/linda/OneDrive/Escritorio/ESTADISTICA”) ##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) Cholula <- read_sav(“Cholula.sav”)

1. Seleccionar variables

vars <- Cholula %>% dplyr::select(X11, X19)

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 = “_“)

3. Renombrar

df <- vars %>% rename_with(~etiquetas)

4. Función para generar gráfica por variable

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_X19]])) + 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_X19 ) + theme_minimal(base_size = 14) }

5. Generar todas las gráficas

graficas <- map(names(df)[-1], graficar_cor)

graficas

Ajustaré el modelo en función de la longitud de pierna

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) datos <- datos[-73, ]

coef <- summary(modelo)$coefficients round(coef, 4)

op <- 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) par(op) # regresa a la configuración original

Graficando modelo de regresión

ggplot(datos, aes(x = X19, y = X11)) + geom_point(alpha = 0.6) + geom_smooth(method = “lm”, se = TRUE, color = “blue”) + labs( title = paste0(lab_X11, ” en función de “, lab_X19), x = lab_X19, y = lab_X11 ) + theme_minimal(base_size = 14)

Se va a realizar el diagnóstico del modelo

Calcular valores diagnósticos

Extraer solo datos usados por el modelo

datos_modelo <- model.frame(modelo)

Diagnósticos

diagnosticos <- data.frame( estatura = datos_modelo\(X11, pierna = 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 —-”) print(outliers_residuales)

cat(“—- PUNTOS CON ALTO LEVERAGE —-”) print(outliers_leverage)

cat(“—- PUNTOS INFLUYENTES (COOK’S DISTANCE) —-”) print(outliers_cook)

Ordenar por Cook’s distance (más influyentes primero) —

diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ] head(diagnosticos_ordenados, 10)

Ajustar el modelo robusto

library(MASS) modelo_robusto <- rlm(X11 ~ X19, data = datos) summary(modelo_robusto)

Comparación de los modelos

coef(modelo) coef(modelo_robusto)

library(ggplot2)

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)

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)

Vamos 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

anova(modelo)

Vamos a ajustar el modelo sin valores atípicos o influyentes

resid_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

datos_limpios <- datos[-indices_fuera, ]

modelo_limpio <- lm(X11 ~ X19, data = datos_limpios) summary(modelo_limpio)

Comparando los modelos

summary(modelo)\(r.squared # R² original summary(modelo_limpio)\)r.squared # R² sin atípicos

anova(modelo) sqrt(1257)