Análisis Discriminante y Regresión Logística en una muestra mexicana
Ruth Sua Cruz Méndez
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 Nicolas Tolentino y San Lorenzo Tezonco.
- Abriendo paquetes estableciendo directorio de trabajo y leyendo archivo
setwd("~/Metodos_cuantitativos_sua")
##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 50Podemos observar que la muestra consta de 80 individuos de los cuales 50 son hombres y 30 mujeres
Estadística descriptiva
Empezemos con la estadística 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 Longuitud Máxima del húmero derecho. También se vislumbra que los hombres tienen valores mayores que las mujeres
Ahora veremos el caso de la Altura biomecánica del húmero
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 en este caso de la tabla anterior observamos que en ambos sexos hay datos perdidos en la variable Altura biomecánica del húmero. También se vislumbra que los hombres tienen valores mayores que las mujeres
REalizaremos una grafica comparativa de longuitud máxima del húmero derercho 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()`).
calcular la 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.427154Realizaremos la prueba t para comparar las medias
## Realizaremos la 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 0cuando las varianas son iguales =Ho cuando son difererntes Ha Pueba de homogeneidad de varianzas
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.8811p-valuehombre=0.8811 es mayor a 0.05 por lo que se asume distribucion nomrmalp-valuemujer=0.9149 es mayor a 0.05 por lo que se asume distribucion normal Ambos grupos se asume distribucion normal por lo que se realiza prueba t-student
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.2842770.4642 es mayor a 0.05, no se rechaza Ho, se asume igualdad de las varianzas Asumiendo varianzas iguales, pruebas t
P-VALUE 1.156e-11 P<0.05 por lo que se rechaza H0, las medias de hombre y mujeres son significativas
#Análisis discriminante
En este apartado aplicaremos la técnica de analisis discriminante que nos permitirá clasificar a un individuo como Hombre o Mujer, dependiendo de los valores discriminantes obtenidos, un paso intermedio déspues de estimar la función discriminante es obtener el punto de corte.
##Ajustar Disciminante 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 dos variables
##Ajustar Disciminante 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 logistica
Aplicaremos la regresión logistica 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 *LMHDAquí aplicaremos la regresión lineal en 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