Objetivo

Ajustar un modelo de regresión lineal para estimar la estatura a partir de la longitud del miembro inferior utilizando la base de datos Cholula.

1. Carga de datos y paquetes

# Definir directorio de trabajo
setwd("C:/Users/OSCAR/Desktop/Estadistica forense EMAC")

# Cargar paquetes
library(pacman)
p_load(haven, dplyr, ggplot2, tinytex, tidyr, GGally, lmtest, car)

# Leer base de datos
Cholula <- read_sav("Cholula.sav")
str(Cholula)
## tibble [339 × 13] (S3: tbl_df/tbl/data.frame)
##  $ Folio  : num [1:339] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "format.spss")= chr "F8.0"
##  $ CEDULA : num [1:339] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ ORIGEN : chr [1:339] "Juvenil Cholula" "Juvenil Cholula" "Juvenil Cholula" "Juvenil Cholula" ...
##   ..- attr(*, "format.spss")= chr "A34"
##   ..- attr(*, "display_width")= int 34
##  $ Origen2: dbl+lbl [1:339] 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
##    ..@ label      : chr "Lugar de origen"
##    ..@ format.spss: chr "F8.0"
##    ..@ labels     : Named num [1:6] 1 2 3 4 5 6
##    .. ..- attr(*, "names")= chr [1:6] "Juvenil Cholula" "Juvenil San Nicolás de los Rancho" "Juvenil Santa Isabel Cholula" "Juvenil diversa procedencia" ...
##  $ X11    : num [1:339] 1715 1701 1652 1691 1611 ...
##   ..- attr(*, "label")= chr "Estatura total"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "display_width")= int 12
##  $ X14    : num [1:339] 771 740 743 775 749 722 764 711 756 721 ...
##   ..- attr(*, "label")= chr "Longitud del miembro superior"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X15    : num [1:339] 318 324 322 342 320 335 329 293 321 316 ...
##   ..- attr(*, "label")= chr "Longitud del brazo"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X16    : num [1:339] 250 216 236 226 234 200 226 212 237 209 ...
##   ..- attr(*, "label")= chr "Longitud del antebrazo"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X17    : num [1:339] 203 200 185 207 195 187 209 206 198 196 ...
##   ..- attr(*, "label")= chr "Longitud de la mano"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X18    : num [1:339] 970 930 942 988 969 896 960 922 985 907 ...
##   ..- attr(*, "label")= chr "Longitud del miembro inferior"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "display_width")= int 12
##  $ X19    : num [1:339] 505 460 505 512 513 472 491 492 523 455 ...
##   ..- attr(*, "label")= chr "Longitud del muslo"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X20    : num [1:339] 378 382 355 400 363 356 383 347 381 380 ...
##   ..- attr(*, "label")= chr "Longitud de la pierna"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12
##  $ X21    : num [1:339] 260 267 255 256 242 262 259 247 271 242 ...
##   ..- attr(*, "label")= chr "Longitud del pie"
##   ..- attr(*, "format.spss")= chr "F4.0"
##   ..- attr(*, "display_width")= int 12

2. Análisis exploratorio

2.1 Matriz de correlación

# Seleccionar variables numéricas continuas
Cholula_num <- Cholula %>% select(X11, X14, X15, X16, X17, X18, X19, X20, X21)

# Matriz de correlación
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

Interpretación: La variable X18 (Longitud del miembro inferior) presenta la correlación más alta con X11 (Estatura) con un valor de 0.899.

2.2 Matriz de gráficos de dispersión

# Matriz de gráficos
ggpairs(Cholula_num)

Interpretación: Los gráficos de dispersión confirman visualmente la fuerte relación lineal entre X11 y X18.

2.3 Correlación específica X11 - X18

# Definir vector con la variable de interés
vars_longitud <- c("X18")

# Calcular correlación
cor_estatura <- sapply(Cholula[vars_longitud], 
                       function(x) cor(Cholula$X11, x, use = "pairwise.complete.obs"))
print(cor_estatura)
##       X18 
## 0.8990134
sort(round(cor_estatura, 3), decreasing = TRUE)
##   X18 
## 0.899

2.4 Estadísticas descriptivas

# Estadísticas de X11 (Estatura)
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)))

# Estadísticas de X18 (Longitud miembro inferior)
res_X18 <- Cholula %>%
  summarise(
    n = sum(!is.na(X18)),
    media = mean(X18, na.rm = TRUE),
    sd = sd(X18, na.rm = TRUE)
  ) %>%
  mutate(across(c(media, sd), ~round(.x, 2)))

cat("Estatura (X11):\n")
## Estatura (X11):
print(res_X11)
## # A tibble: 1 × 3
##       n media    sd
##   <int> <dbl> <dbl>
## 1   339 1611.  59.3
cat("\nLongitud del miembro inferior (X18):\n")
## 
## Longitud del miembro inferior (X18):
print(res_X18)
## # A tibble: 1 × 3
##       n media    sd
##   <int> <dbl> <dbl>
## 1   339  902.  45.2

Interpretación: La estatura promedio es de 1611 mm (SD=59.3) y la longitud promedio del miembro inferior es de 902 mm (SD=45.2). Ambas variables presentan variabilidad moderada.

2.5 Pruebas de normalidad

# Shapiro-Wilk
p_norm_X11 <- shapiro.test(Cholula$X11)
p_norm_X18 <- shapiro.test(Cholula$X18)

cat("Normalidad X11 (Estatura):\n")
## Normalidad X11 (Estatura):
print(p_norm_X11)
## 
##  Shapiro-Wilk normality test
## 
## data:  Cholula$X11
## W = 0.98548, p-value = 0.001739
cat("\nNormalidad X18 (Longitud miembro inferior):\n")
## 
## Normalidad X18 (Longitud miembro inferior):
print(p_norm_X18)
## 
##  Shapiro-Wilk normality test
## 
## data:  Cholula$X18
## W = 0.9888, p-value = 0.01044

Interpretación: Ambas variables rechazan la hipótesis nula de normalidad (p<0.05).

3. Modelo de regresión lineal simple

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

ggplot(Cholula, aes(x = X18, y = X11)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  labs(title = "Relación entre Longitud del Miembro Inferior y Estatura",
       x = "Longitud del Miembro Inferior (mm)",
       y = "Estatura (mm)") +
  theme_minimal(base_size = 13)

Interpretación: Se observa una relación lineal positiva clara entre la longitud del miembro inferior y la estatura. La banda de confianza (área sombreada) indica el rango de incertidumbre del modelo.

3.2 Ajuste del modelo

# Ajustar modelo lineal
modelo <- lm(X11 ~ X18, data = Cholula)

# Resumen del modelo
summary(modelo)
## 
## Call:
## lm(formula = X11 ~ X18, data = Cholula)
## 
## 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
# Coeficientes
coef(modelo)
## (Intercept)         X18 
##  547.546604    1.179751

Interpretación del modelo:

  • Intercepto (547.55): Cuando X18=0, la estatura estimada sería 547.55 mm (valor teórico sin interpretación práctica)
  • Pendiente (1.18): Por cada mm de aumento en la longitud del miembro inferior, la estatura aumenta 1.18 mm
  • R² = 0.8082: El modelo explica el 80.82% de la variabilidad de la estatura
  • p-value < 2e-16: El modelo es estadísticamente significativo

Ecuación del modelo:
Estatura = 547.55 + 1.18 × (Longitud miembro inferior)

4. Diagnóstico del modelo

4.1 Gráficos de diagnóstico

par(mfrow = c(2, 2))
plot(modelo, which = 1)
plot(modelo, which = 2)
plot(modelo, which = 3)
plot(modelo, which = 4)

par(mfrow = c(1, 1))

Interpretación de los gráficos:

  1. Residuos vs Ajustados: Los residuos se distribuyen aleatoriamente alrededor de cero sin patrones claros, sugiriendo linealidad y homocedasticidad adecuadas.

  2. Q-Q Plot: La mayoría de los puntos siguen la línea diagonal, aunque hay ligeras desviaciones en los extremos. Esto indica una normalidad aproximada de los residuos.

  3. Escala-Localización: La línea roja es relativamente horizontal, confirmando homocedasticidad (varianza constante).

  4. Distancia de Cook: No hay observaciones con valores extremadamente altos (>1.0), indicando que no hay puntos altamente influyentes que distorsionen el modelo.

4.2 Pruebas formales de supuestos

# Independencia de errores (Durbin-Watson)
cat("Prueba de Durbin-Watson (Independencia):\n")
## Prueba de Durbin-Watson (Independencia):
dwtest(modelo)
## 
##  Durbin-Watson test
## 
## data:  modelo
## DW = 1.8769, p-value = 0.1248
## alternative hypothesis: true autocorrelation is greater than 0
# Homocedasticidad (Breusch-Pagan)
cat("\nPrueba de Breusch-Pagan (Homocedasticidad):\n")
## 
## Prueba de Breusch-Pagan (Homocedasticidad):
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 0.46976, df = 1, p-value = 0.4931
# Normalidad de residuos (Shapiro-Wilk)
cat("\nPrueba de Shapiro-Wilk (Normalidad de residuos):\n")
## 
## Prueba de Shapiro-Wilk (Normalidad de residuos):
shapiro.test(residuals(modelo))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo)
## W = 0.98185, p-value = 0.0002814

Interpretación de las pruebas:

  • Durbin-Watson (DW=1.88, p=0.12): No se rechaza H0, los residuos son independientes ✓
  • Breusch-Pagan (BP=0.47, p=0.49): No se rechaza H0, hay homocedasticidad ✓
  • Shapiro-Wilk (W=0.98, p=0.0003): Se rechaza H0, los residuos no son completamente normales.

4.3 Valores influyentes

influenceIndexPlot(modelo, vars = c("Cook", "Studentized", "Hat"), id.n = 3)

Interpretación: Se identifican las observaciones 253 y 335 como potencialmente influyentes, pero sus valores de Cook son bajos (<0.5), por lo que no representan un problema grave.

4.4 Diagnóstico visual adicional

qqPlot(modelo, main = "Q-Q plot de residuos")

## [1] 253 335
plot(fitted(modelo), residuals(modelo),
     main = "Residuos vs Valores Ajustados",
     xlab = "Valores Ajustados",
     ylab = "Residuos")
abline(h = 0, col = "red")

Interpretación: El gráfico Q-Q adicional y el de residuos vs ajustados confirman que el modelo lineal simple es apropiado, con ligeras desviaciones que no comprometen seriamente los supuestos.

5. Modelos alternativos

5.1 Modelo cuadrático

# Ajustar modelo cuadrático
modelo2 <- lm(X11 ~ X18 + I(X18^2), data = Cholula)
summary(modelo2)
## 
## Call:
## lm(formula = X11 ~ X18 + I(X18^2), data = Cholula)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -113.349  -17.552   -0.885   15.239  107.259 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.135e+03  3.900e+02   2.911  0.00385 **
## X18         -1.139e-01  8.569e-01  -0.133  0.89431   
## I(X18^2)     7.101e-04  4.700e-04   1.511  0.13178   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.96 on 336 degrees of freedom
## Multiple R-squared:  0.8095, Adjusted R-squared:  0.8084 
## F-statistic:   714 on 2 and 336 DF,  p-value: < 2.2e-16
# Comparar modelos
anova(modelo, modelo2)
## Analysis of Variance Table
## 
## Model 1: X11 ~ X18
## Model 2: X11 ~ X18 + I(X18^2)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    337 228010                           
## 2    336 226471  1    1538.5 2.2825 0.1318

Interpretación: El modelo cuadrático (R²=0.8095) mejora ligeramente al modelo lineal (R²=0.8082). Sin embargo, la prueba ANOVA (p=0.13) indica que esta mejora NO es estadísticamente significativa, por lo tanto el término cuadrático no es necesario.

ggplot(Cholula, aes(x = X18, y = X11)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), se = TRUE, color = "blue") +
  labs(title = "Modelo cuadrático: Estatura vs Longitud Miembro Inferior",
       x = "Longitud del Miembro Inferior (mm)",
       y = "Estatura (mm)") +
  theme_minimal()

5.2 Modelo cúbico

# Ajustar modelo cúbico
modelo3 <- lm(X11 ~ X18 + I(X18^2) + I(X18^3), data = Cholula)
summary(modelo3)
## 
## Call:
## lm(formula = X11 ~ X18 + I(X18^2) + I(X18^3), data = Cholula)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -110.97  -17.83   -1.01   14.71  106.03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -6.308e+03  4.494e+03  -1.404   0.1614  
## X18          2.441e+01  1.478e+01   1.652   0.0995 .
## I(X18^2)    -2.616e-02  1.617e-02  -1.618   0.1066  
## I(X18^3)     9.787e-06  5.887e-06   1.662   0.0974 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.89 on 335 degrees of freedom
## Multiple R-squared:  0.8111, Adjusted R-squared:  0.8094 
## F-statistic: 479.4 on 3 and 335 DF,  p-value: < 2.2e-16
# Comparar con modelos anteriores
anova(modelo, modelo2)
## Analysis of Variance Table
## 
## Model 1: X11 ~ X18
## Model 2: X11 ~ X18 + I(X18^2)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    337 228010                           
## 2    336 226471  1    1538.5 2.2825 0.1318
anova(modelo2, modelo3)
## Analysis of Variance Table
## 
## Model 1: X11 ~ X18 + I(X18^2)
## Model 2: X11 ~ X18 + I(X18^2) + I(X18^3)
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1    336 226471                              
## 2    335 224618  1    1853.1 2.7637 0.09736 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación: El modelo cúbico (R²=0.8111) tampoco mejora significativamente (p=0.097) respecto al cuadrático.

6. Conclusiones

  1. Modelo seleccionado: El modelo lineal simple (X11 ~ X18) es el más apropiado por su simplicidad y buen ajuste.

  2. Capacidad predictiva: El modelo explica el 80.82% de la variabilidad de la estatura, lo cual sería relativamente bueno, recordar que la utilidad práctica depende de otros factores.

  3. Validez del modelo: Cumple satisfactoriamente con los supuestos de independencia y homocedasticidad. La ligera desviación de normalidad no es crítica con esta muestra grande.

  4. Aplicación forense: Este modelo puede utilizarse para estimar la estatura de individuos cuando solo se dispone de la longitud del miembro inferior, con un error estándar residual de ±26.01 mm.

# Guardar análisis completo
save(Cholula, modelo, modelo2, modelo3, file = "Cholula.RData")