Análisis Discriminante y Regresión Logística en una muestra mexicana
Nohemi Angel
2025-11-16
R Markdown
Es una herramienta que combina código de programación (principalmente de R, pero también de otros lenguajes como Python) con texto narrativo para crear documentos reproducibles. Permite integrar código, sus resultados (como tablas y gráficos) y explicaciones en un único archivo, el cual se puede exportar fácilmente a múltiples formatos como PDF, HTML o Word.
INTRODUCCION
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("~/Blanca")
##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 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 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….
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()`).
Análisis discriminante
En este apartado aplicaremos la técnica de análisis discriminante que nos permitirá clasificara 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 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 discirminante 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 - cutoffRegresió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 *LMHDAquí vamos a aplicar la regresión lineal 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