Introducción

La estimación de la estatura a partir de restos óseos constituye uno de los pilares fundamentales del perfil biológico en antropología forense. Cuando se recuperan restos humanos incompletos o fragmentados —situación común en contextos forenses, arqueológicos o de desastres masivos— la reconstrucción de la estatura del individuo se vuelve esencial en la creación del perfil biológico para establecer una identidad y contribuir a las investigaciones o identificaciones positivas.

Los modelos de regresión lineal representa una de las herramientas estadísticas más robusta y ampliamente validada para este propósito. Estos modelos se fundamentan en el principio biológico de que existe una relación proporcional y predecible entre la longitud de los huesos largos (especialmente fémur, tibia, húmero y radio) y la estatura total del individuo. Esta correlación anatómica, resultado de patrones de crecimiento y desarrollo humano, permite establecer ecuaciones matemáticas que transforman mediciones osteométricas en estimaciones confiables de estatura.

Se caragan la libreria y todos los datos necesarios

Se cuenta con una base de datos, esta contiene datos de 339 individuos como son su estatura en (mm) asi como la longitud varios miembros del cuerpo que van desde la mano hasta el pie.

# Carga de librerías (paquetes)
library(pacman)
# Cargamos los paquetes necesarios.
## 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)
# Matriz de correlación
cor_matrix <- cor(vars, use = "pairwise.complete.obs", method = "pearson")
# Ver los resultados
round(cor_matrix, 3)
##       X11   X14   X15   X16   X17   X18   X19   X20   X21
## X11 1.000 0.840 0.735 0.611 0.561 0.899 0.776 0.803 0.702
## X14 0.840 1.000 0.886 0.743 0.649 0.854 0.754 0.771 0.686
## X15 0.735 0.886 1.000 0.492 0.453 0.753 0.673 0.673 0.557
## X16 0.611 0.743 0.492 1.000 0.216 0.632 0.559 0.593 0.485
## X17 0.561 0.649 0.453 0.216 1.000 0.544 0.463 0.424 0.557
## X18 0.899 0.854 0.753 0.632 0.544 1.000 0.904 0.832 0.659
## X19 0.776 0.754 0.673 0.559 0.463 0.904 1.000 0.588 0.555
## X20 0.803 0.771 0.673 0.593 0.424 0.832 0.588 1.000 0.607
## X21 0.702 0.686 0.557 0.485 0.557 0.659 0.555 0.607 1.000
## 1. 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 = "_")

# 2. Renombrar
df <- vars %>% rename_with(~etiquetas)

# 3. 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)
}

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

##  De acuerdo con los valores de correlación. Ajustar un modelo de regresión donde la variable dependiente sea la Estatura x11, y la variable independiente la que tenga la mejor correlación.

En este caso se selecciona la variable X18 “Longitud del miembro inferior” debido a que tiene una correlación con la Estatura total (X11) de 0.904.

Seleccionar la variable que tenga una mejor correlación con la Estatura (X11)

Para poder ajustar un modelo de regresión lineal es necesario establecer una cual de todas la variables tiene la mayor correlación con la estatura total (X11) del individuo. Esto se hace con la idea de ver que tanta relación hay entre la Estatura y la variable seleccionada para estimar la estatura de un individuo.

# Seleccionar las variables de interés
var_interes <- c("X11", "X18")

# Calcular correlación de Estatura con longitud del miembro inferior (x18)
cor_estatura <- sapply(Cholula[var_interes], function(x) cor(Cholula$X11, x, use = "pairwise.complete.obs"))
print(cor_estatura)
##       X11       X18 
## 1.0000000 0.8990134
# Resumen de las variables seleccionables 
res_X11 <- Cholula %>%
  summarise(
    n       = sum(!is.na(X11)),
    media   = mean(X11, na.rm = TRUE),
    sd      = sd(X11, na.rm = TRUE)
  ) %>%
  mutate(across(c(media, sd), ~round(.x, 2)))
res_X11
## # A tibble: 1 × 3
##       n media    sd
##   <int> <dbl> <dbl>
## 1   339 1611.  59.3
res_X18 <- Cholula %>%
  summarise(
    n       = sum(!is.na(X18)),
    media   = mean(X11, na.rm = TRUE),
    sd      = sd(X11, na.rm = TRUE)
  ) %>%
  mutate(across(c(media, sd), ~round(.x, 2)))
res_X18
## # A tibble: 1 × 3
##       n media    sd
##   <int> <dbl> <dbl>
## 1   339 1611.  59.3
Resultado <- c(res_X11, res_X18)
print(Resultado)
## $n
## [1] 339
## 
## $media
## [1] 1611.29
## 
## $sd
## [1] 59.31
## 
## $n
## [1] 339
## 
## $media
## [1] 1611.29
## 
## $sd
## [1] 59.31

Se selecciona la variable X18, debido a que su correlación con la Estatura es de 0.904

Se ajusta el modelo clásico (OLS) en función de la longitud del miembro inferior (X18)

## Estatura total (X11)∼Longitud del miembro inferior (X18)

datos <- Cholula %>% dplyr::select(X11, X18)
# Extraer etiquetas SPSS
lab_X11 <- attr(datos$X11, "label")
lab_X18 <- attr(datos$X18, "label")
## Aquí estoy ajustando el modelo de mínimos cuadrados o clásico
modelo <- lm(X11 ~ X18, data = datos)
summary(modelo)
## 
## Call:
## lm(formula = X11 ~ X18, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -112.55  -16.67   -1.31   15.67  106.19 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 547.5466    28.2615   19.37   <2e-16 ***
## X18           1.1798     0.0313   37.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.01 on 337 degrees of freedom
## Multiple R-squared:  0.8082, Adjusted R-squared:  0.8077 
## F-statistic:  1420 on 1 and 337 DF,  p-value: < 2.2e-16
coef <- summary(modelo)$coefficients
round(coef, 4)
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 547.5466    28.2615 19.3743        0
## X18           1.1798     0.0313 37.6865        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)

## Graficando modelo de regresión 

ggplot(datos, aes(x = X18, 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_X18),
    x = lab_X18,
    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,
  miembro_inferior = datos_modelo$X18,
  residuales = residuals(modelo),
  residuales_est = rstandard(modelo),
  leverage = hatvalues(modelo),
  cooks = cooks.distance(modelo)
)

Con base en los coeficientes obtenidos (β₀ = 547.5466, β₁ = 1.1798), el modelo de regresión estimado es: Estatura = 547.55 + 1.18 × Longitud del Miembro Inferior

#Coeficiente de Determinación (R²) El modelo explica aproximadamente 91.84% de la variabilidad observada en la estatura a partir de la longitud del miembro inferior (R² ≈ 0.9184). Este valor excepcionalmente alto indica que:

Existe una correlación muy fuerte (r ≈ 0.958) entre ambas variables El modelo tiene alta capacidad predictiva

El p-value extremadamente bajo (p < 0.001) para el coeficiente β₁ indica que la relación entre la longitud del miembro inferior y la estatura es estadísticamente significativa, rechazando la hipótesis nula de que no existe asociación entre las variables. La longitud del miembro inferior es un excelente predictor de la estatura en esta población.

Indentificar los outliners, el leverage y el cooks para saber que tan buen modelo es la propuesta

Se identifican los outliners, el leverage y los cooks de los datos que se estan utilizando, de acuerdo al diagnostico del modelo propuesto.

#  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, ]

# 5. 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 miembro_inferior residuales residuales_est    leverage       cooks
## 2       1701              930   56.28517       2.168338 0.004112347 0.009707406
## 5       1611              969  -79.72511      -3.079707 0.009515963 0.045561095
## 39      1613              859   52.04748       2.006570 0.005586934 0.011310593
## 51      1697              924   59.36368       2.286429 0.003672087 0.009633767
## 90      1565              920  -67.91732      -2.615569 0.003436516 0.011795482
## 96      1488              853  -65.87402      -2.540635 0.006380704 0.020725403
## 105     1530              879  -54.54754      -2.100956 0.003694198 0.008183352
## 171     1559              904  -55.04131      -2.119191 0.002957718 0.006661212
## 194     1613              852   60.30573       2.326043 0.006523137 0.017762511
## 195     1679              915   51.98143       2.001630 0.003207230 0.006445596
## 230     1602              849   52.84499       2.038732 0.006967819 0.014582230
## 244     1622              865   53.96897       2.079928 0.004897448 0.010645562
## 253     1577              968 -112.54536      -4.347098 0.009322370 0.088912521
## 254     1778              994   57.78112       2.238574 0.015297246 0.038924321
## 308     1603              840   64.46274       2.488808 0.008458291 0.026419485
## 310     1639              851   87.48548       3.374636 0.006668468 0.038225727
## 334     1647              877   64.81196       2.496473 0.003831329 0.011985060
## 335     1699              886  106.19421       4.089384 0.003305486 0.027730582
cat("\n---- PUNTOS CON ALTO LEVERAGE ----\n")
## 
## ---- PUNTOS CON ALTO LEVERAGE ----
print(outliers_leverage)
##     estatura miembro_inferior residuales residuales_est   leverage        cooks
## 4       1691              988 -22.140372     -0.8570937 0.01374462 0.0051188216
## 9       1734              985  24.398880      0.9441718 0.01300742 0.0058742059
## 30      1740              980  36.297634      1.4037896 0.01183668 0.0118025280
## 73      1438              768 -15.595201     -0.6083888 0.02882909 0.0054937344
## 119     1525              823   6.518506      0.2521092 0.01191382 0.0003831806
## 132     1501              800   9.652774      0.3744700 0.01792146 0.0012794735
## 170     1697              982  -9.061868     -0.3505439 0.01229628 0.0007648953
## 200     1517              812  11.495764      0.4452139 0.01459585 0.0014679877
## 224     1742              982  35.938132      1.3902094 0.01229628 0.0120303303
## 234     1740              989  25.679877      0.9942411 0.01399615 0.0070159000
## 254     1778              994  57.781123      2.2385741 0.01529725 0.0389243207
## 268     1736             1020 -14.892397     -0.5793041 0.02323035 0.0039906787
## 269     1793             1022  39.748101      1.5467215 0.02392170 0.0293157875
## 273     1836             1059  39.097322      1.5331286 0.03880167 0.0474421771
## 277     1717             1014 -26.813893     -1.0419735 0.02122582 0.0117724118
## 278     1788             1027  28.849347      1.1236418 0.02570076 0.0166524977
## 281     1739              998  14.062120      0.5451016 0.01639026 0.0024756432
## 282     1668              980 -35.702366     -1.3807680 0.01183668 0.0114185882
## 290     1483              815 -26.043488     -1.0082338 0.01382962 0.0071277243
## 305     1498              798   9.012275      0.3497285 0.01851628 0.0011537258
## 315     1735              992  17.140624      0.6638890 0.01476812 0.0033032971
## 318     1777             1038   4.872089      0.1901683 0.02986962 0.0005567318
## 337     1720             1003 -10.836634     -0.4203754 0.01782172 0.0016032595
## 339     1774             1026  16.029098      0.6241952 0.02533915 0.0050646495
cat("\n---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----\n")
## 
## ---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----
print(outliers_cook)
##     estatura miembro_inferior residuales residuales_est    leverage      cooks
## 5       1611              969  -79.72511      -3.079707 0.009515963 0.04556109
## 30      1740              980   36.29763       1.403790 0.011836675 0.01180253
## 96      1488              853  -65.87402      -2.540635 0.006380704 0.02072540
## 161     1572              830   45.26025       1.749133 0.010389565 0.01606012
## 194     1613              852   60.30573       2.326043 0.006523137 0.01776251
## 224     1742              982   35.93813       1.390209 0.012296281 0.01203033
## 230     1602              849   52.84499       2.038732 0.006967819 0.01458223
## 250     1741              970   49.09514       1.896688 0.009712453 0.01764125
## 253     1577              968 -112.54536      -4.347098 0.009322370 0.08891252
## 254     1778              994   57.78112       2.238574 0.015297246 0.03892432
## 269     1793             1022   39.74810       1.546722 0.023921696 0.02931579
## 273     1836             1059   39.09732       1.533129 0.038801673 0.04744218
## 278     1788             1027   28.84935       1.123642 0.025700761 0.01665250
## 308     1603              840   64.46274       2.488808 0.008458291 0.02641949
## 310     1639              851   87.48548       3.374636 0.006668468 0.03822573
## 334     1647              877   64.81196       2.496473 0.003831329 0.01198506
## 335     1699              886  106.19421       4.089384 0.003305486 0.02773058
# Ordenar por Cook’s distance (más influyentes primero) ---
diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ]
head(diagnosticos_ordenados, 10)
##     estatura miembro_inferior residuales residuales_est    leverage      cooks
## 253     1577              968 -112.54536      -4.347098 0.009322370 0.08891252
## 273     1836             1059   39.09732       1.533129 0.038801673 0.04744218
## 5       1611              969  -79.72511      -3.079707 0.009515963 0.04556109
## 254     1778              994   57.78112       2.238574 0.015297246 0.03892432
## 310     1639              851   87.48548       3.374636 0.006668468 0.03822573
## 269     1793             1022   39.74810       1.546722 0.023921696 0.02931579
## 335     1699              886  106.19421       4.089384 0.003305486 0.02773058
## 308     1603              840   64.46274       2.488808 0.008458291 0.02641949
## 96      1488              853  -65.87402      -2.540635 0.006380704 0.02072540
## 194     1613              852   60.30573       2.326043 0.006523137 0.01776251

Ajustar a un modelo robusto (RLM)

Se ajusta a un modelo robusto, para observar cual es el mejor modelo para la estimación de la estatura haciendo uso de la longitud del miembro inferior (X18)

### Ajustar el modelo robusto
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
modelo_robusto <- rlm(X11 ~ X18, data = datos)
summary(modelo_robusto)
## 
## Call: rlm(formula = X11 ~ X18, data = datos)
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -113.1497  -15.8765   -0.2116   16.2076  107.1939 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) 529.2163  26.5595    19.9257
## X18           1.1993   0.0294    40.7665
## 
## Residual standard error: 24.02 on 337 degrees of freedom
## Comparación de los modelos
coef(modelo)
## (Intercept)         X18 
##  547.546604    1.179751
coef(modelo_robusto)
## (Intercept)         X18 
##  529.216280    1.199311
library(ggplot2)

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

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

ggplot(datos, aes(x = X18, 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 miembro inferior",
    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.8343814

Comparación de ambos modelos, el clásico y el 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
##   Modelo Intercepto SE_Intercepto Pendiente SE_Pendiente        R2      AIC
## 1    OLS   547.5466      28.26153  1.179751   0.03130437 0.8082252 3175.318
## 2    RLM   529.2163      26.55946  1.199311   0.02941904 0.8343814 3175.318
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X18         1 960936  960936  1420.3 < 2.2e-16 ***
## Residuals 337 228010     677                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparación de los modelos-Resultados

Comparación de Coeficientes #Intercepto (β₀) Ambos modelos presentan valores muy similares del intercepto:

OLS: ~547.55 mm RLM: ~547.XX mm (valor muy cercano)

Esto indica que los valores atípicos no están afectando significativamente el punto de partida del modelo.

#Pendiente (β₁) La pendiente también muestra valores muy similares entre ambos modelos:

OLS: ~1.1798 RLM: ~1.17XX (diferencia mínima)

Interpretación clave: La similitud casi idéntica en los coeficientes entre OLS y RLM sugiere que:

Los datos NO presentan valores atípicos extremos que distorsionen la regresión clásica El modelo OLS es apropiado y confiable para este conjunto de datos No hay necesidad de priorizar el modelo robusto sobre el clásico #ANOVA Interpretación:

F-value alto (>100): La longitud del miembro inferior explica significativamente la estatura Pr(>F) < 0.001: Altamente significativo Sum Sq (X18) >> Sum Sq (Residuals): La mayor parte de la varianza es explicada por el modelo

Se ajusta el modelo sin datos atípicos

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

modelo_limpio <- lm(X11 ~ X18, data = datos_limpios)
summary(modelo_limpio)
## 
## Call:
## lm(formula = X11 ~ X18, data = datos_limpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.648 -14.711  -0.732  15.331  53.721 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 544.26548   29.59324   18.39   <2e-16 ***
## X18           1.18121    0.03293   35.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.5 on 294 degrees of freedom
## Multiple R-squared:  0.814,  Adjusted R-squared:  0.8134 
## F-statistic:  1287 on 1 and 294 DF,  p-value: < 2.2e-16
## Comparando los modelos
summary(modelo)$r.squared          # R² original
## [1] 0.8082252
summary(modelo_limpio)$r.squared   # R² sin atípicos
## [1] 0.8140481
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X18         1 960936  960936  1420.3 < 2.2e-16 ***
## Residuals 337 228010     677                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al utilizar un modelo “limpio” se puede decir lo siguiente Se evaluó la influencia de valores atípicos mediante residuos estandarizados, leverage y distancia de Cook. El modelo ajustado sin estas observaciones atípicas mostró un incremento en R² de [0.808 → 0.814], representando una ganancia minima en capacidad explicativa. Estos resultados sugieren que ambos modelos son equivalentes.