Practica Estadística Forense
Itzama Ortiz
2025-11-21
Introducción
En esta practica se aplican diversas técnicas de análisis multivariado (Análisis Discriminante Lineal y Regresión Logística) sobre la base de datos de hombro (Datos hombro.sav), utilizando las variables Longitud Máxima de la Clavícula Derecha (\(\text{LMCD}\)) y Circunferencia en el Punto Medio de la Clavícula Derecha (\(\text{C12CD}\)) con el objetivo de evaluar su capacidad para la determinación del sexo.
1. Comparación por Sexo (Prueba t o U de Mann-Whitney)
Se evalúa si existen diferencias estadísticamente significativas entre hombres y mujeres para las variables \(\text{LMCD}\) y \(\text{C12CD}\), seleccionando la prueba adecuada según el supuesto de normalidad (Prueba de Shapiro-Wilk).
## 1. Comparación por sexo
# Función para realizar Shapiro-Wilk y elegir la prueba de comparación
# Supuesto de normalidad, prueba de Shapiro-Wilk hipótesis nula los datos estan distribución normal
shapiro.test(HombroS$LMCD)##
## Shapiro-Wilk normality test
##
## data: HombroS$LMCD
## W = 0.98259, p-value = 0.5472shapiro.test(HombroS$C12CD)##
## Shapiro-Wilk normality test
##
## data: HombroS$C12CD
## W = 0.96789, p-value = 0.1149## Al realizar el Shapiro-Wilk para el supuesto de normalidad de los datos p-value > 0.5
## LMCD p-value= 0.5472, C12CD p-value= 0.1149
## LMCD Cumple con la hipotesis de normalidad por lo que se realiza la prueba T
# Se realiza prueba T-test para comparar la medias
t.test(LMCD ~ sexoN, data = HombroS, var.equal = TRUE)##
## Two Sample t-test
##
## data: LMCD by sexoN
## t = -8.2061, df = 58, p-value = 2.767e-11
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -24.53489 -14.91249
## sample estimates:
## mean in group Mujer mean in group Hombre
## 131.3026 151.0263t.test(LMCD ~ sexoN, data = HombroS, var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: LMCD by sexoN
## t = -8.6158, df = 57.876, p-value = 5.838e-12
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -24.30633 -15.14104
## sample estimates:
## mean in group Mujer mean in group Hombre
## 131.3026 151.0263## C12CD No cumple con el supuesto de normalidad por lo que se realiza la prueba U
## la alternativa no paramétrica prueba u de mann-whitney
wilcox.test(C12CD ~ sexoN, data = HombroS)## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: C12CD by sexoN
## W = 60.5, p-value = 1.477e-08
## alternative hypothesis: true location shift is not equal to 0## data: C12CD by sexoN W = 60.5, p-value = 1.477e-08 alternative hypothesis: true location shift is not equal to 0
## Gráfica comparativa
ggplot(HombroS, aes(x = LMCD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Gráfica 1. Longitud maxima clavicula derecha por sexo",
x = "Longitud clavicula derecha por sexo (mm)",
y = "Densidad"
) +
theme_minimal()
ggplot(HombroS, aes(x = C12CD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Gráfica 2. Circunferencia en punto medio de la clavicula derecha",
x = "Circunferencia en punto medio de la clavicula derecha (mm)",
y = "Densidad"
) +
theme_minimal()
Interpretación: Se debe interpretar si los valores p de las pruebas de
comparación son menores a \(\alpha =
0.05\). Si lo son, se rechaza la hipótesis nula (\(\text{H}_{0}\)) de que no existen
diferencias por sexo y se concluye que la media de la variable es
significativamente diferente entre hombres y mujeres.
Conclusión: Se espera rechazar la hipótesis nula (\(\text{H}_{0}\))5, lo que confirma que las medias de \(\text{LMCD}\) es significativamente diferente entre hombres y mujeres, validando su uso en los modelos de determinación sexual.
2. Función Discriminante con LMCD (LDA)
Se ajusta un modelo de Análisis Discriminante Lineal (\(\text{LDA}\)) para clasificar el sexo utilizando únicamente la variable \(\text{LMCD}\).
## Ajustar Discriminante Lineal
# Eliminar NAs de forma explícita
HombroS <- na.omit(Hombro[, c("sexoN", "LMCD")])
# Ajustar el modelo con la base depurada
## El modelo discriminante es D=a*LMCD+b
lda1 <- lda(sexoN ~ LMCD, data = HombroS)
a_LMCD <- coef(lda1) ## El factor que multiplica a LMCD
pred0 <- predict(lda1,
newdata = data.frame(LMCD = 0))
b_LMCD <- pred0$x # valor de D cuando LMCD=0 es b
pred1 <- predict(lda1) ## Predicción de sexo para todos los valores de
# medias de la función discriminante por grupo
centroide_H <- mean(pred1$x[HombroS$sexoN == "Hombre"])
centroide_M <- mean(pred1$x[HombroS$sexoN == "Mujer"])
# Punto de corte (promedio de centroides, priors iguales)
cutoff_lmcd <- mean(c(centroide_H, centroide_M))
## El modelo discriminante es
cat("Función discriminante: D(x) = ", round(a_LMCD, 4), " * LMCD + ", round(b_LMCD, 4), "\n")## Función discriminante: D(x) = 0.1089 * LMCD + -15.5587cat("Punto de corte=", round(cutoff_lmcd, 4), "Si D>",round(cutoff_lmcd, 4), "es Hombre")## Punto de corte= -0.1791 Si D> -0.1791 es Hombre## Función discriminante: D(x) = 0.1089 * LMCD + -15.5587
## Punto de corte= -0.1791 Si D> -0.1791 es Hombre
# Tabla cruzada
tabla_clas <- table(Observado = HombroS$sexoN,
Predicho = pred1$class)
tabla_clas## Predicho
## Observado Mujer Hombre
## Mujer 22 3
## Hombre 3 32prop_clas <- sum(diag(tabla_clas)) / sum(tabla_clas)*100
cat("El porcentaje de clasificación correcta es ",round(prop_clas,1),"%")## El porcentaje de clasificación correcta es 90 %## Resultados: El % de clasificación correcta de LMCD es de 90%, con un punto de corte de -0.1791 Sí D> -0.1791 es Hombre.Función discriminante: D(x) = 0.1089 * LMCD + -15.5587
## Ejemplo de aplicación (LMCD hipotético, e.g., 135mm)
LMCD_hipotetico <- 135
D_ejemplo_lmcd <- (a_LMCD[1] * LMCD_hipotetico) + b_LMCD[1]
clasificacion_lmcd <- ifelse(D_ejemplo_lmcd > cutoff_lmcd, "Hombre", "Mujer")
cat(paste0("**Ejemplo de Aplicación** (LMCD=", LMCD_hipotetico, " mm):\n"))## **Ejemplo de Aplicación** (LMCD=135 mm):cat(paste0("D = ", round(D_ejemplo_lmcd, 4), ". Clasificación: ", clasificacion_lmcd, "\n"))## D = -0.8507. Clasificación: Mujer# Crea el dataframe para el gráfico
data_scores_lmcd <- data.frame(D_score = pred1$x, sexoN = HombroS$sexoN)# Gráfico para LDA (histograma de puntuaciones con corte)
ggplot(data_scores_lmcd, aes(x = pred1$x, fill = sexoN)) +
geom_density(alpha = 0.45, color = NA) +
geom_vline(xintercept = cutoff_lmcd,
linetype = "dotted",
color = "red",
linewidth = 1) +
labs(title = "Distribución de la Puntuación Discriminante (LMCD)",
subtitle = "Separada por sexo",
x = "Puntuación D",
y = "Densidad") +
scale_fill_brewer(palette = "Set2") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.title = element_blank()
)
Conclusión: El modelo multifactorial \((\text{LMCD} + \text{C12CD})\) debe ofrecer un mayor porcentaje de clasificación correcta que un modelo unifactorial, demostrando el valor de predicción conjunto de la longitud y la circunferencia en punto medio de la clavícula.
3. Función Discriminante con LMCD y C12CD (LDA Múltiple)
Se ajusta un modelo \(\text{LDA}\) para clasificar el sexo utilizando las dos variables \(\text{LMCD}\) y \(\text{C12CD}\) simultáneamente.
# Eliminar NAs de forma explícita
HombroNA <- na.omit(Hombro[, c("sexoN", "LMCD","C12CD")])
# Ajustar el modelo con la base depurada
## El modelo discriminante que queremos ajustar es D=a*LMCD+b*C12CD+c
lda2 <- lda(sexoN ~ LMCD+C12CD, data = HombroNA) ##Ajustar Disciminante Lineal
# Obtener los coeficientes a y b
a1 <- coef(lda2)[1]
a2 <- coef(lda2)[2]
lda2## Call:
## lda(sexoN ~ LMCD + C12CD, data = HombroNA)
##
## Prior probabilities of groups:
## Mujer Hombre
## 0.4166667 0.5833333
##
## Group means:
## LMCD C12CD
## Mujer 131.3026 30.31206
## Hombre 151.0263 35.77151
##
## Coefficients of linear discriminants:
## LD1
## LMCD 0.0727066
## C12CD 0.2559868# Predicción de valores mediante modelo
pred2 <- predict(lda2)
# Matriz de confusión
table(Real = HombroNA$sexoN, Predicho = pred2$class)## Predicho
## Real Mujer Hombre
## Mujer 24 1
## Hombre 1 34# Porcentaje de acierto
mean(pred2$class == HombroNA$sexoN)*100## [1] 96.66667## Porcentaje de acierto de clasificación 96.66%
# valor de D cuando LMCD=0 y C12CD=0 (intercepto)
b <- predict(lda2,
newdata = data.frame(LMCD = 0, C12CD = 0))$x
b## LD1
## 1 -18.95782centroide_H <- mean(pred2$x[HombroNA$sexoN == "Hombre"])
centroide_M <- mean(pred2$x[HombroNA$sexoN == "Mujer"])
cutoff <- mean(c(centroide_H, centroide_M))
cutoff## [1] -0.2359658cat("Función discriminante: D(x) = ", round(a1, 4), " * LMCD + ", round(a2, 4),"*C12CD +",round(b, 4), "\n")## Función discriminante: D(x) = 0.0727 * LMCD + 0.256 *C12CD + -18.9578cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 4), "es Hombre")## Punto de corte= -0.236 Si D> -0.236 es Hombre## Función discriminante: D(x) = 0.0727 * LMCD + 0.256 *C12CD + -18.9578
## Punto de corte= -0.236 Si D> -0.236 es Hombre
## Grafico de dispersión
ggplot(HombroNA, aes(x = LMCD, y = C12CD, color = sexoN)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Discriminante lineal: LMCD vs C12CD")
a1 <- lda2$scaling["LMCD", 1]
a2 <- lda2$scaling["C12CD", 1]
b <- predict(lda2, newdata = data.frame(LMCD = 0, C12CD = 0))$x - cutoff
## Grafico LDA
ggplot(HombroNA, aes(x = LMCD, y = C12CD, color = sexoN)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Discriminante lineal: LMCD vs C12CD")
a1 <- lda2$scaling["LMCD", 1]
a2 <- lda2$scaling["C12CD", 1]
b <- predict(lda2, newdata = data.frame(LMCD = 0, C12CD = 0))$x - cutoff
## Ejemplo de aplicación
# 1. Base sin NA
HombroNA <- na.omit(Hombro[, c("sexoN", "LMCD", "C12CD")])
# 2. Ajustar modelo LDA
lda2 <- lda(sexoN ~ LMCD + C12CD, data = HombroNA)
# 3. Predicciones del modelo
pred2 <- predict(lda2)
# 4. Porcentaje de clasificación correcta
porcentaje_clasificacion <- mean(pred2$class == HombroNA$sexoN) * 100
cat("Porcentaje de clasificación correcta =",
round(porcentaje_clasificacion, 2), "%\n\n")## Porcentaje de clasificación correcta = 96.67 %# 5. Obtener coeficientes de la función discriminante
a1 <- lda2$scaling["LMCD", 1]
a2 <- lda2$scaling["C12CD", 1]
# 6. Intercepto (valor de D con LMCD=0, C12CD=0)
b_intercepto <- predict(lda2,
newdata = data.frame(LMCD = 0, C12CD = 0))$x
# 7. Puntos de corte usando centroides
centroide_H <- mean(pred2$x[HombroNA$sexoN == "Hombre"])
centroide_M <- mean(pred2$x[HombroNA$sexoN == "Mujer"])
cutoff <- mean(c(centroide_H, centroide_M))
cat("Punto de corte =", round(cutoff, 4), "\n\n")## Punto de corte = -0.236cat("Función discriminante:\n")## Función discriminante:cat("D(x) = ", round(a1, 4), "* LMCD +",
round(a2, 4), "* C12CD +",
round(b_intercepto, 4), "\n\n")## D(x) = 0.0727 * LMCD + 0.256 * C12CD + -18.9578# 8. Calcular D(x) para LMCD = 135 y C12CD = 20
LMCD_val <- 135
C12CD_val <- 20
D_x <- a1 * LMCD_val + a2 * C12CD_val + b_intercepto
cat("D(135, 20) =", round(D_x, 4), "\n")## D(135, 20) = -4.0227# 9. Clasificación usando el punto de corte
clasificacion <- ifelse(D_x > cutoff, "Hombre", "Mujer")
cat("Clasificación del caso nuevo:", clasificacion, "\n")## Clasificación del caso nuevo: Mujer## Clasificación del caso MUJER4. Modelo Logístico con LMCD (GLM)
Se ajusta un modelo de Regresión Logística (\(\text{GLM}\)) para estimar la probabilidad de ser hombre en función de \(\text{LMCD}\).
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMCD","C12CD")])
# Ajustar modelo logístico binario
modelo1 <- glm(sexoN ~ LMCD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
summary(modelo1)##
## Call:
## glm(formula = sexoN ~ LMCD, family = binomial(link = "logit"),
## data = Hombro_sinNA)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -41.72697 10.79320 -3.866 0.000111 ***
## LMCD 0.29945 0.07717 3.880 0.000104 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 81.503 on 59 degrees of freedom
## Residual deviance: 32.430 on 58 degrees of freedom
## AIC: 36.43
##
## Number of Fisher Scoring iterations: 6 # Probabilidades de ser "Hombre"
Hombro_sinNA$prob_Hombre2 <- predict(modelo1, type = "response")
# Clasificación usando 0.5 como punto de corte
Hombro_sinNA$predicho2 <- ifelse(Hombro_sinNA$prob_Hombre2 >= 0.5, "Hombre", "Mujer")
# Matriz de confusión
table(Real = Hombro_sinNA$sexoN, Predicho2 = Hombro_sinNA$predicho2)## Predicho2
## Real Hombre Mujer
## Mujer 3 22
## Hombre 32 3 # Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho2) * 100## [1] 90 ggplot(Hombro_sinNA, aes(x = prob_Hombre2, fill = sexoN)) +
geom_density(alpha = 0.45, linewidth = 1) +
geom_vline(xintercept = 0.5, color = "red", linetype = "dashed", linewidth = 1) +
scale_fill_manual(values = c("Hombre" = "#1f77b4", "Mujer" = "#ff7f0e")) +
theme_minimal(base_size = 14) +
labs(
title = "Distribución de las probabilidades predichas de ser Hombre",
subtitle = "Modelo logístico con LMCD y C12CD",
x = "Probabilidad predicha: P(Hombre)",
y = "Densidad",
fill = "Sexo real"
) +
theme(
plot.title = element_text(face = "bold"),
legend.position = "top"
)
coef(modelo1)## (Intercept) LMCD
## -41.7269684 0.2994501 # Ajustar el modelo logístico
glm_lmcd <- glm(sexoN ~ LMCD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
# Coeficientes
coef_glm_lmcd <- coef(modelo1)
# Ejemplo de aplicación (LMCD hipotético, e.g., 140)
LMCD_hipotetico_glm <- 140
logit_ejemplo_glm_lmcd <- coef_glm_lmcd[1] + (coef_glm_lmcd[2] * LMCD_hipotetico_glm)
prob_ejemplo_glm_lmcd <- exp(logit_ejemplo_glm_lmcd) / (1 + exp(logit_ejemplo_glm_lmcd))
clasificacion_glm_lmcd <- ifelse(prob_ejemplo_glm_lmcd >= 0.5, "Hombre", "Mujer")
cat(paste0("**Ejemplo de Aplicación** (LMCD=", LMCD_hipotetico_glm, " mm):\n"))## **Ejemplo de Aplicación** (LMCD=140 mm):Conclusión: Los modelos \(\text{GLM}\) proveen una estimación de probabilidad, lo que es útil en casos forenses. El modelo que combina ambas variables debe ser el más preciso.
5. Modelo Logístico con LMCD y C12CD (GLM Múltiple)
Se ajusta un modelo de Regresión Logística (\(\text{GLM}\)) utilizando las dos variables \(\text{LMCD}\) y \(\text{C12CD}\).
# Ajustar modelo logístico binario
modelo2 <- glm(sexoN ~ LMCD + C12CD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
summary(modelo2)##
## Call:
## glm(formula = sexoN ~ LMCD + C12CD, family = binomial(link = "logit"),
## data = Hombro_sinNA)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -69.0885 22.9591 -3.009 0.00262 **
## LMCD 0.3067 0.1117 2.747 0.00602 **
## C12CD 0.8022 0.3405 2.356 0.01847 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 81.503 on 59 degrees of freedom
## Residual deviance: 18.731 on 57 degrees of freedom
## AIC: 24.731
##
## Number of Fisher Scoring iterations: 8 # Probabilidades de ser "Hombre"
Hombro_sinNA$prob_Hombre <- predict(modelo2, type = "response")
# Clasificación usando 0.5 como punto de corte
Hombro_sinNA$predicho <- ifelse(Hombro_sinNA$prob_Hombre >= 0.5, "Hombre", "Mujer")
# Matriz de confusión
table(Real = Hombro_sinNA$sexoN, Predicho = Hombro_sinNA$predicho)## Predicho
## Real Hombre Mujer
## Mujer 1 24
## Hombre 34 1 # Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100## [1] 96.66667 ## Gráfica de probabilidades predichas
ggplot(Hombro_sinNA, aes(x = prob_Hombre, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
theme_minimal() +
labs(title = "Probabilidades predichas de ser Hombre",
x = "P(Hombre)", y = "Densidad")
coef(modelo2)## (Intercept) LMCD C12CD
## -69.0884686 0.3067252 0.8022468 cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo2)[1], 4), " + ",
round(coef(modelo2)[2], 4), "*LMCD + ",
round(coef(modelo2)[3], 4), "*C12CD\n")## Ecuación logística:
## logit(p) = -69.0885 + 0.3067 *LMCD + 0.8022 *C12CD# Ejemplo de aplicación
# 1. Extraer coeficientes del modelo ya ajustado
b0 <- coef(modelo2)[1]
b1 <- coef(modelo2)[2]
b2 <- coef(modelo2)[3]
# 2. Valores específicos
LMCD_val <- 155
C12CD_val <- 25
# 3. Calcular logit y probabilidad
logit_p <- b0 + b1 * LMCD_val + b2 * C12CD_val
probabilidad <- exp(logit_p) / (1 + exp(logit_p))
cat("Probabilidad predicha de ser Hombre =", round(probabilidad, 4), "\n")## Probabilidad predicha de ser Hombre = 0.1839 # 4. Clasificación según punto de corte = 0.5
clasificacion <- ifelse(probabilidad >= 0.5, "Hombre", "Mujer")
cat("Clasificación predicha =", clasificacion, "\n\n")## Clasificación predicha = Mujer # 5. Porcentaje de clasificación correcta del modelo completo
Hombro_sinNA$prob_Hombre <- predict(modelo2, type = "response")
Hombro_sinNA$predicho <- ifelse(Hombro_sinNA$prob_Hombre >= 0.5, "Hombre", "Mujer")
porcentaje_clasificacion <- mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100
cat("Porcentaje de clasificación correcta del modelo =",
round(porcentaje_clasificacion, 2), "%\n")## Porcentaje de clasificación correcta del modelo = 96.67 %