Introducción

Un modelo de regresión lineal estadístico es una técnica que describe la relación entre una variable dependiente y una o más variables independientes mediante una línea recta. Su objetivo es modelar esta relación para hacer predicciones o comprender la influencia de las variables independientes sobre la dependiente.

Componentes principales Variable dependiente (y): La variable que se intenta predecir o explicar. Variable independiente (x): La variable que se utiliza para predecir o explicar la variable dependiente.

Idea básica

La regresión lineal asume que la relación entre las variables se puede aproximar mediante una línea recta (si hay una sola variable independiente) o un plano/hiperplano (si hay varias).

Forma general

Regresión lineal simple (una sola variable predictora):

Y = β0 + β1 (X)

Donde:

Y = variable dependiente Xi = variables independientes β0 = intercepto βi = coeficientes que indican cuánto cambia Y cuando Xi cambia en una unidad

Desarrollo

De acuerdo a lo anterior, crearé un modelo de regresión lineal para la relación entre la estatura y la longitud del miembro superior (mano, antebrazo y brazo), donde la variable dependiente es la estatura y la variable independiente es la longitud del miembro superior.

El objetivo es aplicar modelos de regresión lineal para estimar la altura a partir de la longitud de la extremidad superior.

Los datos utilizados son de una base de datos de mediciones a hombres que realizaron su servicio militar en Cholula, Puebla. Las medidas se encuentran en milmetros

La base está en formato .sav y para trabajar en R se le cambió el nombre a Cholula1, después se eliminaron las entradas no númericas para poder manejar los datos de mejor forma. La variable estatura está codificada con el nombre de X11 y la variable longitud del miembro superior tiene como nombre X14.

setwd("C:/Users/house/OneDrive/Escritorio/ENAH/A. Forense/Estadistica/Carpeta trabajos R")
Cholula1 <- read_sav("Cholula.sav")
Cholula_num <- Cholula1 %>%
  dplyr::select("X11", "X14","X15","X16","X17","X18","X19","X20","X21")
cor_matrix <- cor(Cholula_num, use = "pairwise.complete.obs", method = "pearson")
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

Se observa que las variables tienen un alto valor de correlación con 0.840

##De acuerdo a los valores de correlacion observados ajustaré un modelo de regresion con la variable dependiente V.D= X11 y la variable independiente V.1= X14 (LONGITUD DEL MIEMBRO SUPERIOR)
vars_interes <- c("X11", "X14")
# CORRELACION DE ESTATURA CON LONGITUD DEL MIEMBRO SUPERIOR
cor_longitud <- sapply(Cholula1[vars_interes], function(x) cor(Cholula1$X11, x, use = "pairwise.complete.obs"))
print(cor_longitud)
##       X11       X14 
## 1.0000000 0.8401463
# Crear un data frame para los resultados
res_X11 <- Cholula1 %>%
  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_X14 <- Cholula1 %>%
  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)))

resultados <- c(res_X11, res_X14)
print(resultados)   
## $n
## [1] 339
## 
## $media
## [1] 1611.29
## 
## $sd
## [1] 59.31
## 
## $n
## [1] 339
## 
## $media
## [1] 1611.29
## 
## $sd
## [1] 59.31

Se va a utilizar la prueba Shapiro-Wilk de normalidad, es decir, evalúa si los valores de una variable se comportan como lo harían si provinieran de una distribución normal.

# Normalidad (Shapiro-Wilk)
p_norm_X11 <- shapiro.test(Cholula1$X11)
p_norm_X14 <- shapiro.test(Cholula1$X14)
table(p_norm_X11)
## , , method = Shapiro-Wilk normality test, data.name = Cholula1$X11
## 
##                    p.value
## statistic           0.00173911189385936
##   0.985484575672881                   1
table(p_norm_X14)
## , , method = Shapiro-Wilk normality test, data.name = Cholula1$X14
## 
##                    p.value
## statistic           0.234950359545882
##   0.994250738833757                 1

La variable estatura (X11) tiene P.VALUE= 0.0017 < 0.005, por lo tanto se rechaza la hipótesis de normalidad

La variable longitud del miembro superior (X14) tiene P.VALUE= 0.2349 > 0.005,por lo tanto se confirma la hipótesis de normalidad.

Gráfica de dispersión con línea de tendencia

ggplot(Cholula1, aes(x = X14, y = X11)) +
  geom_point(alpha = 0.6, color = "steelblue") +   # puntos
  geom_smooth(method = "lm", se = FALSE, color = "red") +  # línea de regresión
  labs(title = "Relación entre Estatura y Longitud de miembro superior",
       x = "Longitud de miembro superior (mm.)",
       y = "Estatura (mm.)") +
  theme_minimal(base_size = 13)
## `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()`).

Ajuste del modelo lineal

modelo <- lm(X11 ~ X14, data = Cholula1)

summary(modelo)
## 
## Call:
## lm(formula = X11 ~ X14, data = Cholula1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -81.700 -20.537  -2.229  19.974 102.219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 489.59116   39.60080   12.36   <2e-16 ***
## X14           1.55286    0.05477   28.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.24 on 335 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7058, Adjusted R-squared:  0.705 
## F-statistic: 803.9 on 1 and 335 DF,  p-value: < 2.2e-16

Los resultados son: β0 = 489.5911 βi = 1.5528 R-cuadrada = 0.7058

Lo cual nos deja con el modelo siguiente:

Estatura = 489.5911 + 1.5528 (longitud del miembro superior)

Gráfico 1 Residuos vs ajustados

plot(modelo, which = 1)

Gráfico 2 Q-Q plot

plot(modelo, which = 2)

Gráfico 3 Distancia de Cook

plot(modelo, which = 4)

Gráfico 4 Escala-Localización

plot(modelo, which = 3)

Pruebas de los supuestos 1. Independemcia de los errores Dubin-Watson

dwtest(modelo)
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.7237, p-value = 0.00525
## alternative hypothesis: true autocorrelation is greater than 0
  1. Homocedasticidad Prueba de Breusch-Pagan
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 4.867, df = 1, p-value = 0.02737
  1. Normalidad de los residuos
shapiro.test(residuals(modelo))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo)
## W = 0.9932, p-value = 0.1318

Ajuste del modelo

# Seleccionar variables

vars <- Cholula1 %>% dplyr::select(X11, X14:X20)

# 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 = "_")

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

#  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)
}
datos <- Cholula1 %>% dplyr::select(X11, X14)
# Extraer etiquetas SPSS
lab_X11 <- attr(datos$X11, "label")
lab_X14 <- attr(datos$X14, "label")
modelo <- lm(X11 ~ X14, data = datos)
summary(modelo)
## 
## Call:
## lm(formula = X11 ~ X14, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -81.700 -20.537  -2.229  19.974 102.219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 489.59116   39.60080   12.36   <2e-16 ***
## X14           1.55286    0.05477   28.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.24 on 335 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7058, Adjusted R-squared:  0.705 
## F-statistic: 803.9 on 1 and 335 DF,  p-value: < 2.2e-16

Volviendo a graficar el modelo

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

Gráfico 5 Residuales estandarizados con el modelo ajustado

op <- par() 
par(mar = c(4, 4, 2, 1)) 
par(mfrow = c(2, 2))
plot(modelo)

Diagnóstico del modelo

Extracción de los datos usados por el modelo

datos_modelo <- model.frame(modelo)

Diagnósticos

diagnosticos <- data.frame(
  estatura = datos_modelo$X11,
  lmsp = datos_modelo$X14,
  residuales = residuals(modelo),
  residuales_est = rstandard(modelo),
  leverage = hatvalues(modelo),
  cooks = cooks.distance(modelo)
)

Posibles outliners 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, ]

Identificar puntos influyentes por distancia de Cook — Regla común: Cook > 4/n

umbral_cook <- 4 / n
outliers_cook <- diagnosticos[diagnosticos$cooks > umbral_cook, ]

Resultados

cat("---- OUTLIERS POR RESIDUALES ESTANDARIZADOS ----\n")
## ---- OUTLIERS POR RESIDUALES ESTANDARIZADOS ----
print(outliers_residuales)
##     estatura lmsp residuales residuales_est    leverage       cooks
## 9       1734  756   70.44548       2.191929 0.006239960 0.015084236
## 34      1649  697   77.06431       2.396163 0.004818594 0.013900190
## 50      1555  732  -71.28584      -2.214734 0.003237437 0.007965679
## 72      1620  681   72.91010       2.270510 0.007896481 0.020516029
## 96      1488  686  -66.85421      -2.080747 0.006775904 0.014768235
## 143     1541  728  -79.07440      -2.456493 0.003060262 0.009261704
## 157     1631  779  -68.27034      -2.130688 0.012237272 0.028121712
## 224     1742  761   70.68117       2.200419 0.007283973 0.017763316
## 238     1661  689  101.48721       3.157694 0.006172826 0.030965863
## 249     1578  650   79.04881       2.474390 0.018064962 0.056319724
## 269     1793  783   87.51822       2.733281 0.013591983 0.051471270
## 273     1836  826   63.74517       2.011732 0.033987979 0.071195597
## 278     1788  792   68.54246       2.144333 0.016977760 0.039707399
## 305     1498  702  -81.69999      -2.539460 0.004159796 0.013468993
## 320     1567  737  -67.05015      -2.083505 0.003588782 0.007817495
## 323     1622  684   70.25151       2.186958 0.007206818 0.017359440
## 335     1699  713  102.21853       3.175731 0.003218399 0.016281614
## 339     1774  784   66.96536       2.091769 0.013945092 0.030939825
cat("\n---- PUNTOS CON ALTO LEVERAGE ----\n")
## 
## ---- PUNTOS CON ALTO LEVERAGE ----
print(outliers_leverage)
##     estatura lmsp residuales residuales_est   leverage        cooks
## 16      1672  778 -25.717475     -0.8024997 0.01191302 3.882278e-03
## 19      1702  787  -9.693230     -0.3029515 0.01503905 7.006768e-04
## 21      1673  783 -32.481783     -1.0144385 0.01359198 7.090024e-03
## 52      1535  666  11.203023      0.3496215 0.01212405 7.500868e-04
## 73      1438  629 -28.341098     -0.8917049 0.02810493 1.149676e-02
## 75      1523  663   3.861608      0.1205734 0.01312541 9.667708e-05
## 103     1530  658  18.625916      0.5820949 0.01490979 2.564208e-03
## 152     1493  663 -26.138392     -0.8161351 0.01312541 4.429401e-03
## 157     1631  779 -68.270337     -2.1306881 0.01223727 2.812171e-02
## 198     1687  780 -13.823199     -0.4314882 0.01256729 1.184792e-03
## 202     1536  663  16.861608      0.5264803 0.01312541 1.843254e-03
## 245     1652  787 -59.693230     -1.8656481 0.01503905 2.657241e-02
## 249     1578  650  79.048809      2.4743898 0.01806496 5.631972e-02
## 250     1741  778  43.282525      1.3506075 0.01191302 1.099652e-02
## 254     1778  791  60.095324      1.8796846 0.01657847 2.978138e-02
## 264     1522  632  51.000317      1.6033279 0.02651480 3.500852e-02
## 269     1793  783  87.518217      2.7332812 0.01359198 5.147127e-02
## 273     1836  826  63.745168      2.0117322 0.03398798 7.119560e-02
## 278     1788  792  68.542462      2.1443328 0.01697776 3.970740e-02
## 281     1739  790  22.648185      0.7082569 0.01618496 4.126196e-03
## 282     1668  797 -59.221846     -1.8547059 0.01906078 3.342093e-02
## 290     1483  653 -20.609776     -0.6447257 0.01683848 3.559574e-03
## 306     1546  662  28.414469      0.8873577 0.01347074 5.375872e-03
## 315     1735  797   7.778154      0.2435957 0.01906078 5.765113e-04
## 318     1777  799  46.672431      1.4623354 0.01993439 2.174762e-02
## 334     1647  780 -53.823199     -1.6800798 0.01256729 1.796238e-02
## 339     1774  784  66.965355      2.0917693 0.01394509 3.093982e-02
cat("\n---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----\n")
## 
## ---- PUNTOS INFLUYENTES (COOK'S DISTANCE) ----
print(outliers_cook)
##     estatura lmsp residuales residuales_est    leverage      cooks
## 9       1734  756   70.44548       2.191929 0.006239960 0.01508424
## 14      1622  768  -60.18886      -1.875387 0.008988025 0.01594913
## 30      1740  773   50.04683       1.560473 0.010378372 0.01276857
## 34      1649  697   77.06431       2.396163 0.004818594 0.01390019
## 72      1620  681   72.91010       2.270510 0.007896481 0.02051603
## 96      1488  686  -66.85421      -2.080747 0.006775904 0.01476824
## 157     1631  779  -68.27034      -2.130688 0.012237272 0.02812171
## 165     1615  763  -59.42455      -1.850409 0.007741983 0.01335775
## 224     1742  761   70.68117       2.200419 0.007283973 0.01776332
## 234     1740  764   64.02259       1.993825 0.007979647 0.01598848
## 238     1661  689  101.48721       3.157694 0.006172826 0.03096586
## 245     1652  787  -59.69323      -1.865648 0.015039051 0.02657241
## 249     1578  650   79.04881       2.474390 0.018064962 0.05631972
## 254     1778  791   60.09532       1.879685 0.016578474 0.02978138
## 264     1522  632   51.00032       1.603328 0.026514798 0.03500852
## 268     1736  770   50.70542       1.580328 0.009526847 0.01201077
## 269     1793  783   87.51822       2.733281 0.013591983 0.05147127
## 273     1836  826   63.74517       2.011732 0.033987979 0.07119560
## 278     1788  792   68.54246       2.144333 0.016977760 0.03970740
## 282     1668  797  -59.22185      -1.854706 0.019060776 0.03342093
## 305     1498  702  -81.69999      -2.539460 0.004159796 0.01346899
## 318     1777  799   46.67243       1.462335 0.019934388 0.02174762
## 323     1622  684   70.25151       2.186958 0.007206818 0.01735944
## 334     1647  780  -53.82320      -1.680080 0.012567291 0.01796238
## 335     1699  713  102.21853       3.175731 0.003218399 0.01628161
## 339     1774  784   66.96536       2.091769 0.013945092 0.03093982

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

diagnosticos_ordenados <- diagnosticos[order(-diagnosticos$cooks), ]
head(diagnosticos_ordenados, 10)
##     estatura lmsp residuales residuales_est    leverage      cooks
## 273     1836  826   63.74517       2.011732 0.033987979 0.07119560
## 249     1578  650   79.04881       2.474390 0.018064962 0.05631972
## 269     1793  783   87.51822       2.733281 0.013591983 0.05147127
## 278     1788  792   68.54246       2.144333 0.016977760 0.03970740
## 264     1522  632   51.00032       1.603328 0.026514798 0.03500852
## 282     1668  797  -59.22185      -1.854706 0.019060776 0.03342093
## 238     1661  689  101.48721       3.157694 0.006172826 0.03096586
## 339     1774  784   66.96536       2.091769 0.013945092 0.03093982
## 254     1778  791   60.09532       1.879685 0.016578474 0.02978138
## 157     1631  779  -68.27034      -2.130688 0.012237272 0.02812171

Creación de modelo robusto

La regresión robusta es un conjunto de métodos diseñados para abordar los problemas que surgen cuando se violan los supuestos tradicionales del modelo clásico, particularmente cuando los datos contienen valores atípicos u observaciones influyentes.

La característica principal es su capacidad para manejar valores extremos sin que estos sesguen fuertemente los coeficientes del modelo, a diferencia del modelo clásico, que intenta minimizar el cuadrado de todos los residuos y es fuertemente afectado por valores grandes y distantes. Los métodos robustos a menudo funcionan asignando diferentes ponderaciones a las observaciones. A las observaciones atípicas se les da menos peso o se ignoran en el cálculo del ajuste, reduciendo su influencia.

El método robusto mejora la confiabilidad de las predicciones y las inferencias cuando los datos reales no se ajustan perfectamente a una distribución normal o presentan variabilidad desigual (heterocedasticidad).

modelo_robusto <- rlm(X11 ~ X14, data = datos)
summary(modelo_robusto)
## 
## Call: rlm(formula = X11 ~ X14, data = datos)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.632 -19.468  -1.752  20.667 103.130 
## 
## Coefficients:
##             Value    Std. Error t value 
## (Intercept) 478.5071  39.7489    12.0383
## X14           1.5671   0.0550    28.5063
## 
## Residual standard error: 29.7 on 335 degrees of freedom
##   (2 observations deleted due to missingness)

Los resultados son: β0 = 478.5071 βi = 1.5671

Comparación de los modelos

coef(modelo)
## (Intercept)         X14 
##  489.591156    1.552862
coef(modelo_robusto)
## (Intercept)         X14 
##  478.507145    1.567129

Gráfico 6 modelo robusto vs modelo clásico

coef_rob <- coef(modelo_robusto)

coef_ols <- coef(modelo)

ggplot(datos, aes(x = X14, 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 miembro superior",
    y = "Estatura"
  ) +
  
  theme_minimal(base_size = 14)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## `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()`).

Pseudo R2

El pseudo R2 es una medida utilizada en modelos de regresión, especialmente en la regresión logística, para evaluar la bondad de ajuste del modelo. Se calcula comparando la verosimilitud del modelo estimado con la de un modelo nulo (sin variables predictoras). Al igual que el R2 convencional, un valor más alto de pseudo R2 indica que el modelo explica una mayor proporción de la varianza en la variable dependiente.

Evalúa la capacidad predictiva de un modelo. Un valor más alto sugiere un mejor ajuste.Se basa en la razón de verosimilitud, que compara el modelo con una versión simplificada, como un modelo nulo sin variables predictoras.Indica la proporción de variación explicada por el modelo.En la regresión logística, no necesariamente alcanza un valor máximo de 1, pero los valores más altos son preferibles para comparar modelos.

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

Comparación de 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   489.5912      39.60080  1.552862   0.05477002 0.7058458 3301.286
## 2    RLM   478.5071      39.74886  1.567129   0.05497480 0.7467848 3301.286
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X14         1 835509  835509  803.86 < 2.2e-16 ***
## Residuals 335 348190    1039                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

modelo_limpio <- lm(X11 ~ X14, data = datos_limpios)
summary(modelo_limpio)
## 
## Call:
## lm(formula = X11 ~ X14, data = datos_limpios)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.646 -19.726  -2.458  19.487 102.739 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 484.6642    42.8473   11.31   <2e-16 ***
## X14           1.5590     0.0593   26.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.23 on 291 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.7037, Adjusted R-squared:  0.7027 
## F-statistic: 691.2 on 1 and 291 DF,  p-value: < 2.2e-16

Comparación de modelos con y sin atípicos

summary(modelo)$r.squared          # R² original
## [1] 0.7058458
summary(modelo_limpio)$r.squared   # R² sin atípicos
## [1] 0.7037178
anova(modelo)
## Analysis of Variance Table
## 
## Response: X11
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## X14         1 835509  835509  803.86 < 2.2e-16 ***
## Residuals 335 348190    1039                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación

Los resultados del primer modelo son:

β0 = 489.5911 βi = 1.5528 R2 = 0.7058

Lo cual nos deja con el modelo siguiente:

Estatura = 489.5911 + 1.5528 (longitud del miembro superior)

De las pruebas realizadas a los supuestos y a los residuos se obtuvo:

  1. Durbin-Watson test

DW = 1.7237, p-value = 0.00525 La prueba de Durbin-Watson es un test estadístico para detectar la autocorrelación en los residuos de un modelo de regresión lineal. Si el valor es menor que 2, existe autocorrelación positiva. Los residuos grandes tienden a ir seguidos de residuos grandes, y los pequeños por pequeños. Lo cual sucede en este caso.

  1. Breusch-Pagan

BP = 4.867, df = 1, p-value = 0.02737 Esta prueba mide la homocesdasticidad. La homocedasticidad es un supuesto en estadística, especialmente en modelos de regresión, que establece que la varianza de los errores (residuos) es constante a lo largo de todas las observaciones. Hipótesis nula (H0): Existe homocedasticidad. Si el valor p es bajo (menor a 0.05), se rechaza la hipótesis nula y se concluye que existe heterocedasticidad. En este modelo se rechaza la H0 ya que el p-value es de 0.02737, lo que dice que la varianza de los errores no es constante.

  1. Shapiro-Wilks

W = 0.9932, p-value = 0.1318 El p-value es mayor a 0.05, por lo tanto existe normalidad de los residuos.

Debido a estas pruebas se hizo un nuevo modelo ajustando los datos, pero los resultados continuaban igual. Se realizó un modelo robusto para conseguir datos más confiables y se obtuvó:

β0 = 478.5071 βi = 1.5671 R2 = 0.7467

Lo que nos indica que el modelo robusto tiene un mayor porcentaje de precicción con el 74%, siendo una mejoría en el porcentaje de predicción.Sin embargo, al momento de quitar valores atípicos, el valor de R2 cae a 0.7037 en comparación del modelo inicial; lo cual nos indica que el modelo inicial es un mejor predictor para la altura.

Entonces, a partir del modelo:

altura = 489.5911 + 1.5528 × (longitud del miembro superior)

  1. Interpretación de los coeficientes

β0 = 489.5911 Es la altura estimada cuando la longitud del miembro superior es 0. No tiene interpretación biológica, pero es necesario matemáticamente para calcular predicciones.

β1 = 1.5528 Significa que por cada 1 mm adicional de longitud del miembro superior, la altura aumenta en promedio 1.5528 mm.

Es decir, al aumentar el miembro superior en 10 mm, la altura aumentaría alrededor de 15.5 mm.

  1. Predicción con una longitud hipotética en mm

Longitud del miembro superior = 720 mm (es decir, 72 cm, una medida posible en adultos)

Sustituimos en el modelo:

altura estimada = 489.5911 + 1.5528 (720)

altura estimada = 1607.6071mm = 160.7 cm = 1.67 m

Conclusión

El modelo de regresresión lineal altura = 489.5911 + 1.5528 (longitud del miembro superior), con una R2 = 0.7058, constituye un predictor sólido y estadísticamente adecuado de la estatura individual. El valor de la R2 indica que aproximadamente el 70.58 % de la variación total observada en la altura puede explicarse a partir de las diferencias en la longitud del miembro superior.

Asimismo, el coeficiente de β1 = 1.5528 revela una relación positiva y consistente entre ambas variables: a medida que aumenta la longitud del miembro superior, la altura tiende a incrementarse proporcionalmente. Esta asociación es congruente con los principios de proporcionalidad corporal antropométrica.

Estos resultados permiten afirmar que el modelo posee una buena capacidad predictiva y que la longitud del miembro superior constituye un indicador confiable para la estimación de la estatura en el contexto poblacional analizado.