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
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
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 4.867, df = 1, p-value = 0.02737
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)
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
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:
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.
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.
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)
β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.
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
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.