### El objetivo es aplicar modelos de regresión lineal para 
### estimar la estatura a partir de longitud del muslo base Cholula

## 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, X14:X21)

# 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_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()`).

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

## 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")
# PRIMERO eliminar observaciones atípicas/influyentes
datos <- datos[-c(73, 273, 308, 310), ]

## 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 
## -97.427 -23.242  -1.481  19.462 119.314 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 721.34623   38.65175   18.66   <2e-16 ***
## X19           1.87030    0.08115   23.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.8 on 333 degrees of freedom
## Multiple R-squared:  0.6147, Adjusted R-squared:  0.6135 
## F-statistic: 531.2 on 1 and 333 DF,  p-value: < 2.2e-16
coef <- summary(modelo)$coefficients
round(coef, 4)
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 721.3462    38.6517 18.6627        0
## X19           1.8703     0.0811 23.0488        0
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 
## 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
## 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_X19, " en función de ", lab_X11),
    x = lab_X19,
    y = lab_X11
  ) +
  theme_minimal(base_size = 14)
## `geom_smooth()` using formula = 'y ~ x'

### 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 ----\n")
## ---- OUTLIERS POR RESIDUALES ESTANDARIZADOS ----
print(outliers_residuales)
##     estatura pierna residuales residuales_est    leverage       cooks
## 2       1701    460  119.31388       3.340025 0.004254192 0.023830770
## 50      1555    485  -73.44372      -2.055102 0.003427893 0.007263654
## 56      1588    508  -83.46071      -2.341174 0.008340056 0.023048552
## 79      1487    448  -72.24247      -2.025052 0.006932089 0.014312893
## 89      1565    490  -72.79524      -2.037575 0.004033334 0.008406526
## 95      1488    462  -97.42672      -2.726911 0.003951740 0.014750943
## 104     1530    479  -87.22189      -2.440168 0.003040472 0.009079732
## 253     1778    523   78.48473       2.208431 0.014472304 0.035810166
## 309     1638    531  -76.47771      -2.156574 0.018688229 0.044285329
## 311     1735    495   87.85324       2.460121 0.004895675 0.014887679
## 312     1696    458  118.05449       3.305341 0.004597748 0.025231843
## 316     1567    493  -76.40615      -2.139169 0.004519911 0.010388609
## 330     1647    455   74.66540       2.091137 0.005190152 0.011407089
## 331     1699    461  115.44358       3.231428 0.004097828 0.021483050
## 335     1774    514   91.31746       2.564380 0.010515503 0.034942644
cat("\n---- PUNTOS CON ALTO LEVERAGE ----\n")
## 
## ---- PUNTOS CON ALTO LEVERAGE ----
print(outliers_leverage)
##     estatura pierna residuales residuales_est   leverage        cooks
## 9       1734    523  34.484726     0.97034322 0.01447230 6.913367e-03
## 118     1525    424  10.644831     0.29987120 0.01672712 7.648683e-04
## 131     1501    423 -11.484865    -0.32362376 0.01726370 9.199149e-04
## 139     1645    519 -47.034058    -1.32221267 0.01261097 1.116433e-02
## 151     1493    426 -25.095777    -0.70658850 0.01568480 3.977845e-03
## 183     1670    524 -31.385578    -0.88335841 0.01496333 5.926792e-03
## 199     1517    419  11.996352     0.33842406 0.01951276 1.139644e-03
## 200     1497    427 -22.966081    -0.64645942 0.01517905 3.220622e-03
## 202     1552    422  41.385439     1.16649519 0.01781055 1.233724e-02
## 213     1507    426 -11.095777    -0.31240908 0.01568480 7.776103e-04
## 231     1695    531 -19.477707    -0.54924655 0.01868823 2.872538e-03
## 233     1740    536  16.170773     0.45668658 0.02165715 2.308430e-03
## 253     1778    523  78.484726     2.20843054 0.01447230 3.581017e-02
## 267     1736    530  23.392597     0.65945235 0.01812527 4.013888e-03
## 268     1793    539  63.559860     1.79677449 0.02356182 3.895123e-02
## 275     1717    530   4.392597     0.12383014 0.01812527 1.415309e-04
## 276     1788    537  62.300469     1.76001946 0.02228177 3.529725e-02
## 288     1483    432 -46.317602    -1.30219938 0.01280445 1.099721e-02
## 303     1498    410   9.829088     0.27808834 0.02517431 9.985415e-04
## 304     1546    425  29.774527     0.83854169 0.01620082 5.789618e-03
## 309     1638    531 -76.477707    -2.15657402 0.01868823 4.428533e-02
## 314     1777    545  36.338036     1.02939687 0.02764860 1.506557e-02
## 320     1647    522 -50.644970    -1.42471835 0.01399156 1.440169e-02
## 333     1720    535  -1.958923    -0.05530553 0.02104282 3.287361e-05
cat("\n---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----\n")
## 
## ---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----
print(outliers_cook)
##     estatura pierna residuales residuales_est    leverage      cooks
## 2       1701    460  119.31388       3.340025 0.004254192 0.02383077
## 5       1611    513  -69.81223      -1.960085 0.010127238 0.01965311
## 30      1740    507   70.40959       1.974749 0.008013448 0.01575098
## 36      1611    443   61.10905       1.714308 0.008484611 0.01257419
## 56      1588    508  -83.46071      -2.341174 0.008340056 0.02304855
## 79      1487    448  -72.24247      -2.025052 0.006932089 0.01431289
## 95      1488    462  -97.42672      -2.726911 0.003951740 0.01475094
## 193     1613    441   66.84966       1.876006 0.009177551 0.01629932
## 202     1552    422   41.38544       1.166495 0.017810550 0.01233724
## 223     1742    510   66.79868       1.874431 0.009024101 0.01599741
## 243     1622    444   70.23875       1.970097 0.008153554 0.01595319
## 249     1741    516   54.57685       1.533255 0.011322860 0.01346171
## 253     1778    523   78.48473       2.208431 0.014472304 0.03581017
## 268     1793    539   63.55986       1.796774 0.023561819 0.03895123
## 276     1788    537   62.30047       1.760019 0.022281765 0.03529725
## 279     1739    515   54.44716       1.529295 0.010914043 0.01290340
## 296     1592    435   57.07149       1.603484 0.011502997 0.01496011
## 306     1632    450   69.01692       1.934101 0.006383013 0.01201531
## 307     1598    436   61.20118       1.719152 0.011089733 0.01657154
## 309     1638    531  -76.47771      -2.156574 0.018688229 0.04428533
## 311     1735    495   87.85324       2.460121 0.004895675 0.01488768
## 312     1696    458  118.05449       3.305341 0.004597748 0.02523184
## 314     1777    545   36.33804       1.029397 0.027648603 0.01506557
## 320     1647    522  -50.64497      -1.424718 0.013991555 0.01440169
## 331     1699    461  115.44358       3.231428 0.004097828 0.02148305
## 335     1774    514   91.31746       2.564380 0.010515503 0.03494264
# Ordenar por Cook’s distance (más influyentes primero) ---
diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ]
head(diagnosticos_ordenados, 10)
##     estatura pierna residuales residuales_est    leverage      cooks
## 309     1638    531  -76.47771      -2.156574 0.018688229 0.04428533
## 268     1793    539   63.55986       1.796774 0.023561819 0.03895123
## 253     1778    523   78.48473       2.208431 0.014472304 0.03581017
## 276     1788    537   62.30047       1.760019 0.022281765 0.03529725
## 335     1774    514   91.31746       2.564380 0.010515503 0.03494264
## 312     1696    458  118.05449       3.305341 0.004597748 0.02523184
## 2       1701    460  119.31388       3.340025 0.004254192 0.02383077
## 56      1588    508  -83.46071      -2.341174 0.008340056 0.02304855
## 331     1699    461  115.44358       3.231428 0.004097828 0.02148305
## 5       1611    513  -69.81223      -1.960085 0.010127238 0.01965311
### Ajustar el modelo robusto
library(MASS)
## 
## Adjuntando el paquete: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
modelo_robusto <- rlm(X11 ~ X19, data = datos)
summary(modelo_robusto)
## 
## Call: rlm(formula = X11 ~ X19, data = datos)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -95.8927 -21.6479  -0.2948  20.9991 120.8908 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) 709.9023  38.9410    18.2302
## X19           1.8918   0.0818    23.1400
## 
## Residual standard error: 31.81 on 333 degrees of freedom
## Comparación de los modelos
coef(modelo)
## (Intercept)         X19 
##  721.346233    1.870304
coef(modelo_robusto)
## (Intercept)         X19 
##  709.902337    1.891754
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)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

# 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.7042647
## 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   721.3462      38.65175  1.870304   0.08114547 0.6146932 3351.883
## 2    RLM   709.9023      38.94098  1.891754   0.08175267 0.7042647 3351.883
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X19         1 680816  680816  531.25 < 2.2e-16 ***
## Residuals 333 426755    1282                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 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] 48
datos_limpios <- datos[-indices_fuera, ]

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 
## -62.619 -19.777  -0.846  15.549  72.190 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 700.03240   42.17430   16.60   <2e-16 ***
## X19           1.91017    0.08874   21.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.71 on 285 degrees of freedom
## Multiple R-squared:  0.6192, Adjusted R-squared:  0.6178 
## F-statistic: 463.3 on 1 and 285 DF,  p-value: < 2.2e-16
## Comparando los modelos
summary(modelo)$r.squared          # R² original
## [1] 0.6146932
summary(modelo_limpio)$r.squared   # R² sin atípicos
## [1] 0.619157
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X19         1 680816  680816  531.25 < 2.2e-16 ***
## Residuals 333 426755    1282                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sqrt(1257)
## [1] 35.4542
##Para el caso de la obtención de la estatura por medio de la longitud del muslo, se comprueba que después de quitar cuatro datos atípicos/influyentes el modelo robusto no modificó sustancialmente los coeficientes del modelo OLS, lo que significa que el presente modelo no cuenta con valores atípicos severos ni puntos influyentes que afecten la estimación, por lo tanto, para la obtención de la longitud de la estatura a través del miembro inferior proximal se selecciona el modelo original con la siguiente ecuación:

Estatura= 721.346+/-1.870 (Longitud de muslo)+/-35.45