Tarea de Estadística Forense 2025
## Defino mi directorio de trabajo
setwd("~/Documents/Estadistica")
##Abriendo paquete pacman
library(pacman)
## El archivo esta en formato SPSS, lo abrimos mediante la libreria haven ## MASS permite analisis discriminante
p_load(haven,dplyr,ggplot2,MASS)
p_load(haven,dplyr,ggplot2,tinytex,tidyr, GGally)
Hombro <- read_sav("Datos hombro.sav")
## Definimos como factor la variable sexoN
Hombro$sexoN <- factor(Hombro$sexoN,
levels = c(1, 2),
labels = c("Hombre", "Mujer"))
table (Hombro$sexoN) ## Frecuencias de sexo
##
## Hombre Mujer
## 50 30
# Omitimos los NA de la tabla
HombroS <- na.omit(Hombro[, c("sexoN", "LMCD","C12CD")])
# Para resumir por sexo la información de Longitud Maxima de clavicula derecha
res_lmcd <- Hombro %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(LMCD)),
media = mean(LMCD, na.rm = TRUE),
sd = sd(LMCD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_lmcd
## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Hombre 35 151. 10.2
## 2 Mujer 25 131. 7.56
# Para resumir por sexo la información de Circunferencia media de clavicula derecha
res_C12CD <- HombroS %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(C12CD)),
media = mean(C12CD, na.rm = TRUE),
sd = sd(C12CD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_C12CD
## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Hombre 35 35.8 2.34
## 2 Mujer 25 30.3 2.86
## 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()
## 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()
# 1. Comparaciones por sexo entre las variables
## Realizaremos la prueba t para comparar las 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 Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
## 14.91249 24.53489
## sample estimates:
## mean in group Hombre mean in group Mujer
## 151.0263 131.3026
t.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 Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
## 15.14104 24.30633
## sample estimates:
## mean in group Hombre mean in group Mujer
## 151.0263 131.3026
t.test(C12CD ~ sexoN, data = HombroS, var.equal = TRUE)
##
## Two Sample t-test
##
## data: C12CD by sexoN
## t = 8.1156, df = 58, p-value = 3.922e-11
## alternative hypothesis: true difference in means between group Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
## 4.112868 6.806030
## sample estimates:
## mean in group Hombre mean in group Mujer
## 35.77151 30.31206
t.test(C12CD ~ sexoN, data = HombroS, var.equal = FALSE)
##
## Welch Two Sample t-test
##
## data: C12CD by sexoN
## t = 7.8461, df = 45.117, p-value = 5.68e-10
## alternative hypothesis: true difference in means between group Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
## 4.058095 6.860803
## sample estimates:
## mean in group Hombre mean in group Mujer
## 35.77151 30.31206
# los valores son parametricos por lo que no es necesario hacer la prueba U
## 2. Construir una función discriminante para 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 <- coef(lda1) ## El factor que multiplica a LMCD
pred0 <- predict(lda1,
newdata = data.frame(LMCD = 0))
b <- pred0$x # valor de D cuando LMCD=0 es b
pred <- predict(lda1) ## Predicción de sexo para todos los valores de
# medias de la función discriminante por grupo
centroide_H <- mean(pred$x[HombroS$sexoN == "Hombre"])
centroide_M <- mean(pred$x[HombroS$sexoN == "Mujer"])
# Punto de corte (promedio de centroides, priors iguales)
cutoff <- mean(c(centroide_H, centroide_M))
## El modelo discriminante es
cat("Función discriminante: D(x) = ", round(a, 4), " * LMCD + ", round(b, 4), "\n")
## Función discriminante: D(x) = 0.1089 * LMCD + -15.5587
cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 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 = pred$class)
tabla_clas
## Predicho
## Observado Hombre Mujer
## Hombre 32 3
## Mujer 3 22
prop_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 %
## El porcentaje de clasificación correcta es 90 %
## Ejemplo de aplicación
# Valores del modelo obtenidos
a <- 0.1089 # coeficiente de LMCD
b <- -15.5587 # intercepto
cutoff <- -0.1791 # punto de corte
# Valor dado
LMCD_nuevo <- 130
# Calcular la función discriminante
D_x <- a * LMCD_nuevo + b
# Mostrar el valor calculado
cat("El valor discriminante para LMCD=130 es:", round(D_x, 4), "\n")
## El valor discriminante para LMCD=130 es: -1.4017
# Clasificación según el punto de corte
if (D_x > cutoff) {
cat("Clasificación: Hombre\n")
} else {
cat("Clasificación: Mujer\n")
}
## Clasificación: Mujer
## Se clasifica como mujer cuando LMCD= 130
## 3. Construir un modelo discriminante para las dos varables
##Ajustar Disciminante Lineal para dos variables
# Eliminar NAs de forma explícita
Hombro_sinNA <- 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 = Hombro_sinNA) ##Ajustar Disciminante Lineal
# Obtener los coeficientes a y b
a1 <- coef(lda2)[1]
a2 <- coef(lda2)[2]
lda2
## Call:
## lda(sexoN ~ LMCD + C12CD, data = Hombro_sinNA)
##
## Prior probabilities of groups:
## Hombre Mujer
## 0.5833333 0.4166667
##
## Group means:
## LMCD C12CD
## Hombre 151.0263 35.77151
## Mujer 131.3026 30.31206
##
## 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 = HombroS$sexoN, Predicho = pred2$class)
## Predicho
## Real Hombre Mujer
## Hombre 34 1
## Mujer 1 24
# Porcentaje de acierto
mean(pred2$class == HombroS$sexoN)*100
## [1] 96.66667
# 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.95782
centroide_H <- mean(pred2$x[Hombro_sinNA$sexoN == "Hombre"])
centroide_M <- mean(pred2$x[Hombro_sinNA$sexoN == "Mujer"])
cutoff <- mean(c(centroide_H, centroide_M))
cutoff
## [1] 0.2359658
cat("Función discriminante: D(x) = ", round(a1, 4), " * LMHD + ", round(a2, 4),"*ABHD +",round(b, 4), "\n")
## Función discriminante: D(x) = -0.0727 * LMHD + -0.256 *ABHD + 18.9578
cat("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 * LMHD + -0.256 *ABHD + 18.9578
## Punto de corte= 0.236 Si D> 0.236 es Hombre
ggplot(Hombro_sinNA, aes(x = LMCD, y = C12CD, color = sexoN)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Discriminante lineal: LMCDD vs C12CD")
a1 <- lda2$scaling["LMCD", 1]
a2 <- lda2$scaling["C12CD", 1]
b <- predict(lda2, newdata = data.frame(LMCD = 0, C12CD = 0))$x - cutoff
# Ecuación de la recta: a1*LMCD + a2*C12CD + b = 0
# => C12CD = -(a1/a2)*LMCD - b/a2
LD1 <- pred2$x[, 1]
ggplot(data.frame(LD1, sexoN = Hombro_sinNA$sexoN),
aes(x = LD1, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = cutoff, linetype = "dashed") +
scale_x_continuous(limits = c(-5, 5)) + # <-- escala de -5 a 5
theme_minimal() +
labs(title = "Distribución sobre la función discriminante",
x = "LD1 (puntuación discriminante)",
y = "Densidad")
## Ejemplo de aplicación
# Coeficientes del modelo ya calculado
a1 <- -0.0727 # coeficiente LMCD
a2 <- -0.256 # coeficiente C12CD
c <- 18.9578 # intercepto
cutoff <- 0.236 # punto de corte
# Valores nuevos
LMCD_nuevo <- 140
C12CD_nuevo <- 35
# Calcular la función discriminante
D_x <- a1 * LMCD_nuevo + a2 * C12CD_nuevo + c
# Mostrar resultado
cat("Valor de D(x) para LMCD=140 y C12CD=35:", round(D_x, 4), "\n")
## Valor de D(x) para LMCD=140 y C12CD=35: -0.1802
# Clasificación
if (D_x > cutoff) {
cat("Clasificación: Hombre\n")
} else {
cat("Clasificación: Mujer\n")
}
## Clasificación: Mujer
## 4. Modelo logistico para calcular % de clasificación correcta
## 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
## Hombre 3 32
## Mujer 22 3
# Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho2) * 100
## [1] 10
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., 150)
LMCD_hipotetico_glm <- 150
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=150 mm):
ggplot(Hombro_sinNA, aes(x = prob_Hombre2, fill = sexoN)) +
geom_density(alpha = 0.4) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "black") +
labs(title = "Probabilidades Predichas de ser Hombre (LMCD)",
x = "P(Hombre)",
y = "Densidad") +
theme_minimal()
## 5. Ajustar modelo logistico, calcular el % de clasificación correcta usando clasificación % correcta usando variables LMCD Y 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
## Hombre 1 34
## Mujer 24 1
# Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100
## [1] 3.333333
## 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 <- 160
C12CD_val <- 30
# 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.017
# 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 = 3.33 %