El objetivo es aplicar modelos de regresión lineal para estimar la estatura a partir de longitud (es) de la mano

Primero defino mi directorio de trabajo, se instalan librerías para poder leer el archivo con la variable a analizar

## Defino mi directorio de trabajo
setwd("~/Rstudio")
##Abriendo paquete pacman
library(pacman)
#Abriendo las librerias
## El archivo esta en formato SPSS, lo abrimos mediante la libreria haven
p_load(haven,dplyr,ggplot2,tinytex,tidyr,GGally,purrr,labelled)

Después de leer el archivo SPSS, selecciono las variables de interés

Cholula <- read_sav("Cholula.sav")
View(Cholula)
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                                28
# 1. Seleccionar variables
vars <- Cholula %>% dplyr::select(X11, X14:X21)
View(vars)

Se extraen etiquetas de la base original para poder renombrar las variables en la base delimitada “vars”

# 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)

Se da la instrucción de generar gráficas para observar la correlación de Estatura con diferentes variables, y para poder elegir la variable a analizar y generar un modelo de regresión lineal

# 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_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)
}

# 5. Generar todas las gráficas
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()`).

Interpretación: Se observó en el estadístico de prueba r (0.561) que no hay una correlación significativa de la Estatura total con la Longitud de la mano, lo que nos indica que no existe relación entre la estatura y la longitud de la mano. Sin embargo, se realizará el modelo de regresión lineal ajustado con el fin de ejemplificar el ajuste de un modelo en el programa R.

Se selecciona de la base de datos Cholula las variables Estatura total y Longitud de la mano, creando un objeto que contenga las etiquetas de interés para poder comenzar a ajustar el modelo de Estimación de estatura a partir de la longitud de la mano

# Ajustaré el modelo en función de la longitud de la mano

## Estatura total (X11)∼Longitud de la mano (X17)

datos <- Cholula %>% dplyr::select(X11, X17)
# Extraer etiquetas SPSS
lab_X11 <- attr(datos$X11, "label")
lab_X17 <- attr(datos$X17, "label")

Se creará el modelo con los datos de las variables Estatura total y Longitud de la mano, se obtienen los residuales y los coeficientes para poder construir la ecuación del modelo.

## Aquí estoy ajustando el modelo de mínimos cuadrados o clásico
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-16
anova(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 ' ' 1
coef <- 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        0
sqrt(1247) #raíz cuadrada de la varianza de los errores (desviación estándar)
## [1] 35.31289
op <- par()   # guarda la configuración actual del panel de gráfica
par(mar = c(3, 3, 2, 1))  # reduce márgenes
par(mfrow = c(1, 1))
plot(modelo)

par(op) # regresa a la configuración original
## Warning in par(op): el parámetro del gráfico "cin" no puede ser especificado
## Warning in par(op): el parámetro del gráfico "cra" no puede ser especificado
## Warning in par(op): el parámetro del gráfico "csi" no puede ser especificado
## Warning in par(op): el parámetro del gráfico "cxy" no puede ser especificado
## Warning in par(op): el parámetro del gráfico "din" no puede ser especificado
## Warning in par(op): el parámetro del gráfico "page" no puede ser especificado

Se graficó el modelo de regresión de acuerdo a la ecuación

## Graficando modelo de regresión 

ggplot(datos, aes(x = X17, 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_X17),
    x = lab_X17,
    y = lab_X11
  ) +
  theme_minimal(base_size = 14)
## `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()`).

Interpretación Se observó en la gráfica que no hay una correlación significativa entre ambas variables, porque los datos se encuentran muy dispersos enytre sí.

El diagnóstico calcula los valores dentro del modelo creado y nos permite ver las secciones en donde se puede mejorar el modelo

### 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,
  mano = datos_modelo$X17,
  residuales = residuals(modelo),
  residuales_est = rstandard(modelo),
  leverage = hatvalues(modelo),
  cooks = cooks.distance(modelo)
)

Se identificarán los datos atípicos e influyebtes, y a qué distancia se encuentra cada valor de la recta.

#  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 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.026745873
cat("\n---- PUNTOS CON ALTO LEVERAGE ----\n")
## 
## ---- PUNTOS 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-03
cat("\n---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----\n")
## 
## ---- PUNTOS 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.02674587
# Ordenar 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.02082712

Se procede a ajustar el modelo para hacerlo más robusto y que sea más significativo, a su vez, realizaremos una comparación gráfica de los dos modelos para poder observar qué es lo que cambia.

### Ajustar el modelo robusto
library(MASS)
## 
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
modelo_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.284907
coef(modelo_robusto)
## (Intercept)         X17 
##  990.679408    3.274143
library(ggplot2)

# Datos del modelo robusto
coef_rob <- coef(modelo_robusto)

# Datos del modelo OLS
coef_ols <- coef(modelo)

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.4904967
## 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
##   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.605
View(comparacion)

anova(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 ' ' 1

Se ajustará el modelo original, eliminando valores atípicos, y se procederá a la comparación entre el modelo original y el modelo robusto

## 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
## [1] 44
datos_limpios <- datos[-indices_fuera, ]

modelo_limpio <- lm(X11 ~ X17, data = datos_limpios)
summary(modelo_limpio)
## 
## Call:
## lm(formula = X11 ~ X17, data = datos_limpios)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -143.765  -29.564   -0.954   25.200  157.869 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 978.0698    52.5027   18.63   <2e-16 ***
## X17           3.3353     0.2779   12.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 46.01 on 290 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.3318, Adjusted R-squared:  0.3295 
## F-statistic:   144 on 1 and 290 DF,  p-value: < 2.2e-16
## Comparando los modelos
summary(modelo)$r.squared          # R² original
## [1] 0.3142753
summary(modelo_limpio)$r.squared   # R² sin atípicos
## [1] 0.3318124
anova(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 ' ' 1

Interpretación: En el ajuste de modelos No se observaron cambios significativos, aún eliminando los datos atípicos y los extremos, lo que significa que la variable Longitud de la mano no es recomendable para la estimación de Estatura, ya que no se encontró correlación entre las mismas.