Análisis Discriminante y Regresión Logística en una muestra mexicana
Olympia Citlallin González Pulido
2025-11-07
Introducción
En esta práctica se hará uso de una base de datos contemporánea procedente de dos panteones de la Ciudad de México: San Nicolás Tolentino, y San Lorenzo Tezonco.
- Abriendo paquetes estableciendo directorio de trabajo y leyendo archivo
setwd("~/Olympia")
##Abriendo paquete pacman
library(pacman)
## El archivo esta en formato SPSS, lo abrimos mediante la libreria haven
p_load(haven,dplyr,ggplot2,MASS)
p_load(tinytex)
Hombro <- read_sav("Datos hombro.sav")## Definimos como factor la variable sexoN
Hombro$sexoN <- factor(Hombro$sexoN,
levels = c(2, 1),
labels = c("Mujer", "Hombre"))
table (Hombro$sexoN) ## Frecuencias de sexo##
## Mujer Hombre
## 30 50table (Hombro$sexoN) ## Frecuencias de sexo##
## Mujer Hombre
## 30 50Podemos observar que la muestra consta de 80 individuos 50 son hombres y 30 son mujeres.
Estadística descriptiva
Empezemos con la estadistica descriptiva, para resumir por sexo la información de Longitud Máxima del húmero derecho
res_dmhd <- Hombro %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(LMHD)),
media = mean(LMHD, na.rm = TRUE),
sd = sd(LMHD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_dmhd## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Mujer 26 283. 15.7
## 2 Hombre 43 312. 13.9De la tabla anterior observamos que en ambos sexos hay datos perdidos en la variable Longitud Máxima del Húmero Derecho. También se observa que los hombres tienen valores mayores que las mujeres.
Ahora veremos el caso de la Altura biomecánica del húmero derecho
res_abhd <- Hombro %>%
group_by(sexoN) %>%
summarise(
n = sum(!is.na(ABHD)),
media = mean(ABHD, na.rm = TRUE),
sd = sd(ABHD, na.rm = TRUE)
) %>%
mutate(across(c(media, sd), ~round(.x, 2)))
res_abhd## # A tibble: 2 × 4
## sexoN n media sd
## <fct> <int> <dbl> <dbl>
## 1 Mujer 26 279. 14.7
## 2 Hombre 41 308. 13.6Observamos que…
Realizaremos una gráfica comparativa de Longitud Máxima del Húmero Derecho por Sexo
ggplot(Hombro, aes(x = LMHD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Gráfica 1. Longitud máxima del húmero derecho por sexo",
x = "Longitud máxima del húmero derecho (mm)",
y = "Densidad"
) +
theme_minimal()## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_density()`).
Otra gráfica, ahora de caja.
ggplot(Hombro, aes(x = sexoN, y = LMHD, fill = sexoN)) +
geom_boxplot(alpha = 0.7) +
labs(
title = "Longitud máxima del húmero derecho",
x = "Sexo",
y = " "
) +
theme_minimal()## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
Cálculo d de Cohen
p_load(effsize)
# Comparar entre hombres y mujeres
cohen.d(Hombro$LMHD, Hombro$sexoN,na.rm = TRUE) ##
## Cohen's d
##
## d estimate: -2.031314 (large)
## 95 percent confidence interval:
## lower upper
## -2.635474 -1.427154Observamos que…
Prueba t para comparar las medias
t.test(LMHD ~ sexoN, data = Hombro, var.equal = TRUE)##
## Two Sample t-test
##
## data: LMHD by sexoN
## t = -8.1766, df = 67, p-value = 1.156e-11
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -36.86178 -22.39620
## sample estimates:
## mean in group Mujer mean in group Hombre
## 282.5998 312.2288 t.test(LMHD ~ sexoN, data = Hombro, var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: LMHD by sexoN
## t = -7.9259, df = 47.778, p-value = 2.941e-10
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -37.14616 -22.11183
## sample estimates:
## mean in group Mujer mean in group Hombre
## 282.5998 312.2288 ## la alternativa no paramétrica
wilcox.test(LMHD ~ sexoN, data = Hombro)## 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: LMHD by sexoN
## W = 81, p-value = 3.344e-09
## alternative hypothesis: true location shift is not equal to 0Se observa que…
Supuesto de normalidad, prueba de Shapiro-Wilk
by(Hombro$LMHD, Hombro$sexoN, shapiro.test)## Hombro$sexoN: Mujer
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.98209, p-value = 0.9149
##
## ------------------------------------------------------------
## Hombro$sexoN: Hombre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.98631, p-value = 0.8811Se observa que… ¿Hay diferencias significativas entre Hombres y Mujeres en la variable LMHD? # Ho:media de Hombres = media de Mujeres vs. Ha: media de Hombres es diferente a media de Mujeres # ¿Hay distribución normal? Prueba de Shapiro-Wilks by(Hombro\(LMHD, Hombro\)sexoN, shapiro.test) #p-v Hombres: 0.8811 >0.05, por lo que se asume Distribución normal #p-v Mujeres: 0.9149 >0.05, por lo que se asume Distribución normal #En ambos grupos se asume Dn, por lo que realiza la prueba t-student
Ahora se analiza el Supuesto de homogeneidad de varianza
var.test(LMHD ~ sexoN, data = Hombro)##
## F test to compare two variances
##
## data: LMHD by sexoN
## F = 1.2843, num df = 25, denom df = 42, p-value = 0.4642
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6498778 2.7084863
## sample estimates:
## ratio of variances
## 1.284277Cuando las variables son diferentes: #Prueba de homogeneidad de varianzas var.test(LMHD ~ sexoN, data = Hombro) #p-v 0.4642 es> a 0.05, no se rechaza Ho, se asume igualdad de varianzas #Haremos prueba de igualdad de varianzas, Prueba t t.test(LMHD ~ sexoN, data = Hombro, var.equal = TRUE) #p-v 1.156e-11 (0.0000000000001156) #p-v<0.05, por lo que se rechaza Ho
Análisis discriminante
En este apartado aplicaremos la técnica de Análisis Discriminante que nos permitirá clasificar a un individuo como Hombre o Mujer, dependiendo de los valores discriminantes obtenidos. Un paso intermedio después de estimar la función discriminante es obtener el punto de corte.
##Ajustar Discriminante Lineal
# Eliminar NAs de forma explícita
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMHD")])
# Ajustar el modelo con la base depurada
## El modelo discriminante es D=a*LMHD+b
lda1 <- lda(sexoN ~ LMHD, data = Hombro_sinNA)
a <- coef(lda1) ## El factor que multiplica a LMHD
pred0 <- predict(lda1,
newdata = data.frame(LMHD = 0))
b <- pred0$x # valor de D cuando LMHD=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[Hombro_sinNA$sexoN == "Hombre"])
centroide_M <- mean(pred$x[Hombro_sinNA$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), " * LMHD + ", round(b, 4), "\n")## Función discriminante: D(x) = 0.0686 * LMHD + -20.6405 cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 4), "es Hombre")## Punto de corte= -0.2502 Si D> -0.2502 es Hombre # Tabla cruzada
tabla_clas <- table(Observado = Hombro_sinNA$sexoN,
Predicho = pred$class)
tabla_clas## Predicho
## Observado Mujer Hombre
## Mujer 21 5
## Hombre 5 38 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 85.5 %Ejemplo: Supongamos que se encuentra un húmero derecho cuya Longitud Máxima es de 314.5mm. Usaremos el modelo para predecir.
prede1 <- predict(lda1, newdata = data.frame(LMHD = 314.5))
prede1$x## LD1
## 1 0.9211348if(prede1$x>cutoff) cat("Debido a que", prede1$x,">",cutoff, "es hombre") else cat("Debido a que", prede1$x, "<", cutoff,"es mujer")## Debido a que 0.9211348 > -0.2502343 es hombreAhora vamos a ajustar un Discriminante lineal con 2 variables
##Ajustar Discriminante Lineal para dos variables
# Eliminar NAs de forma explícita
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMHD","ABHD")])
# Ajustar el modelo con la base depurada
## El modelo discriminante que queremos ajustar es D=a*LMHD+b*ABHD+c
lda2 <- lda(sexoN ~ LMHD+ABHD, data = Hombro_sinNA) ##Ajustar Disciminante Lineal
# Obtener los coeficientes a y b
a1 <- coef(lda2)[1]
a2 <- coef(lda2)[2]
lda2## Call:
## lda(sexoN ~ LMHD + ABHD, data = Hombro_sinNA)
##
## Prior probabilities of groups:
## Mujer Hombre
## 0.3880597 0.6119403
##
## Group means:
## LMHD ABHD
## Mujer 282.5998 278.9212
## Hombre 312.2643 308.3958
##
## Coefficients of linear discriminants:
## LD1
## LMHD -0.1352937
## ABHD 0.2092925 # Predicción de valores mediante modelo
pred2 <- predict(lda2)
# Matriz de confusión
table(Real = Hombro_sinNA$sexoN, Predicho = pred2$class)## Predicho
## Real Mujer Hombre
## Mujer 21 5
## Hombre 3 38 # Porcentaje de acierto
mean(pred2$class == Hombro_sinNA$sexoN)*100## [1] 88.0597 # valor de D cuando LMHD=0 y ABHD=0 (intercepto)
b <- predict(lda2,
newdata = data.frame(LMHD = 0, ABHD = 0))$x
b## LD1
## 1 -21.46112 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.2412749 cat("Función discriminante: D(x) = ", round(a1, 4), " * LMHD + ", round(a2, 4),"*ABHD +",round(b, 4), "\n")## Función discriminante: D(x) = -0.1353 * LMHD + 0.2093 *ABHD + -21.4611 cat("Punto de corte=", round(cutoff, 4), "Si D>",round(cutoff, 4), "es Hombre")## Punto de corte= -0.2413 Si D> -0.2413 es Hombre ggplot(Hombro_sinNA, aes(x = LMHD, y = ABHD, color = sexoN)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Discriminante lineal: LMHD vs ABHD")
a1 <- lda2$scaling["LMHD", 1]
a2 <- lda2$scaling["ABHD", 1]
b <- predict(lda2, newdata = data.frame(LMHD = 0, ABHD = 0))$x - cutoff
# Ecuación de la recta: a1*LMHD + a2*ABHD + b = 0
# => ABHD = -(a1/a2)*LMHD - b/a2
ggplot(Hombro_sinNA, aes(x = LMHD, y = ABHD, color = sexoN)) +
geom_point(size = 3) +
geom_abline(slope = -(a1/a2), intercept = -b/a2, color = "black", linetype = "dashed") +
theme_minimal() +
labs(title = "Frontera de clasificación LDA")
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")
## Regresión Logística
Aplicaremos la Regresión Logística para determinar la probabilidad de hombre o mujer, usando los huesos del hombro (LMHD y ABHD). Aplicaremos primero la Regresión con una variable (LMHD)
### Ahora aplicaremos Regresión Logística
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMHD")])
# Ajustar modelo logístico binario
modelo1 <- glm(sexoN ~ LMHD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
summary(modelo1)##
## Call:
## glm(formula = sexoN ~ LMHD, family = binomial(link = "logit"),
## data = Hombro_sinNA)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -42.25632 10.21782 -4.136 3.54e-05 ***
## LMHD 0.14368 0.03449 4.166 3.10e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 91.422 on 68 degrees of freedom
## Residual deviance: 45.765 on 67 degrees of freedom
## AIC: 49.765
##
## Number of Fisher Scoring iterations: 6 # Probabilidades de ser "Hombre"
Hombro_sinNA$prob_Hombre <- predict(modelo1, 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 5 21
## Hombre 38 5 # Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100## [1] 85.50725 ## 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(modelo1)## (Intercept) LMHD
## -42.2563211 0.1436795 cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo1)[1], 4), " + ",
round(coef(modelo1)[2], 4), "*LMHD\n")## Ecuación logística:
## logit(p) = -42.2563 + 0.1437 *LMHDAhora aplicaremos la Regresión logística con dos variables
### Ahora aplicaremos Regresión Logística
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMHD","ABHD")])
# Ajustar modelo logístico binario
modelo2 <- glm(sexoN ~ LMHD + ABHD,
data = Hombro_sinNA,
family = binomial(link = "logit"))
summary(modelo2)##
## Call:
## glm(formula = sexoN ~ LMHD + ABHD, family = binomial(link = "logit"),
## data = Hombro_sinNA)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -50.0490 12.6492 -3.957 7.6e-05 ***
## LMHD -0.4774 0.3048 -1.566 0.1173
## ABHD 0.6555 0.3294 1.990 0.0466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 89.495 on 66 degrees of freedom
## Residual deviance: 39.992 on 64 degrees of freedom
## AIC: 45.992
##
## Number of Fisher Scoring iterations: 6 # 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 5 21
## Hombre 38 3 # Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100## [1] 88.0597 ## 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) LMHD ABHD
## -50.0490482 -0.4773712 0.6555261 cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo2)[1], 4), " + ",
round(coef(modelo2)[2], 4), "*LMHD + ",
round(coef(modelo2)[3], 4), "*ABHD\n")## Ecuación logística:
## logit(p) = -50.049 + -0.4774 *LMHD + 0.6555 *ABHD