Práctica 2. Construcción del modelo de regresión logístico tomando en cuenta las variables Longitud Máxima de la Clavícula Derecha y Circunferencia en punto medio de la Clavícula Derecha
Olympia Citlallin González Pulido
2025-12-10
Primero defino mi directorio de trabajo
setwd("D:/Métodos cuantitativos 2/Olympia")Instalo los paquetes necesarios para abrir el archivo .sav y poder analizar la información
##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")Convertir variable sexo a factor
## 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 50Empezamos con la estadística descriptiva, para resumir por sexo la información de la variable Longitud Máxima de la Clavícula 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 Mujer 25 131. 7.56
## 2 Hombre 35 151. 10.2Analizamos lo mismo, ahora con la variable Circunferencia en Punto Medio de la Clavícula Derecha
res_c12cd <- Hombro %>%
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 Mujer 29 30.3 2.83
## 2 Hombre 46 35.9 2.33##Realizamos representaciones visaules a traves de graficas de las estadísticas descriptivas
## Gráfica comparativa lmcd
ggplot(Hombro, aes(x = LMCD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Gráfica 1. Longitud máxima de la clavícula derecha por sexo",
x = "Longitud máxima de la clavícula derecha (mm)",
y = "Densidad"
) +
theme_minimal()## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_density()`).
## Gráfica comparativa c12cd
ggplot(Hombro, aes(x = C12CD, fill = sexoN)) +
geom_density(alpha = 0.5) +
labs(
title = "Gráfica 1. Circunferencia del punto medio de la clavícula derecha por sexo",
x = "Circunferencia del punto medio de la clavícula derecha (mm)",
y = "Densidad"
) +
theme_minimal()## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_density()`).
## Gráfica de caja o box plot lmcd
ggplot(Hombro, aes(x = sexoN, y = LMCD, fill = sexoN)) +
geom_boxplot(alpha = 0.7) +
labs(
title = "Longitud máxima de la clavícula derecha",
x = "Sexo",
y = " "
) +
theme_minimal()## Warning: Removed 20 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Gráfica de caja o box plot c12cd
ggplot(Hombro, aes(x = sexoN, y = C12CD, fill = sexoN)) +
geom_boxplot(alpha = 0.7) +
labs(
title = "Circunferencia del punto medio de la clavícula derecha",
x = "Sexo",
y = " "
) +
theme_minimal()## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
##1. Realizar comparaciones por sexo de las dos variables (Prueba t O U-de Mann Whitney)
Como primer paso, se realiza Prueba de normalidad Shapiro-Wilk para saber si tiene o no distribución normal, H0: existe distribución normal, vs. HA: no existe distribución normal se cumple el supuesto de normalidad si p valor es mayor a .05
# Supuesto de normalidad, prueba de Shapiro-Wilk lmcd
by(Hombro$LMCD, Hombro$sexoN, shapiro.test)## Hombro$sexoN: Mujer
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97687, p-value = 0.8169
##
## ------------------------------------------------------------
## Hombro$sexoN: Hombre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.94899, p-value = 0.1055#Resultados:
p-value = 0.8169 Mujer, se distribuye normal
p-value = 0.1055 Hombre, se distribuye normalPara ambos sexos, la variable Longitud Máxima de la Clavícula Derecha se distribuye normal, porque p valor es mayor que 0.05 en ambos sexos
#Prueba de Normalidad para la variable C12CD:# Supuesto de normalidad, prueba de Shapiro-Wilk c12cd
by(Hombro$C12CD, Hombro$sexoN, shapiro.test)## Hombro$sexoN: Mujer
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95394, p-value = 0.2313
##
## ------------------------------------------------------------
## Hombro$sexoN: Hombre
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.95329, p-value = 0.06285 p-value = 0.2313 Mujer, se distribuye normal
p-value = 0.06285 Hombre, se distribuye normalPara ambos sexos, la variable Circunferencia del Punto Medio de la Clavícula Derecha se distribuye normal, porque p valor es mayor que 0.05 en ambos sexos
Entonces para ambas variables tenemos, por sexo, distribución normal, por lo que se aplicará el estadístico de prueba t de student.
#Procederemos primero con la variable LMCD, y después con la variable C12CD
## Realizaremos la prueba t para comparar las medias
t.test(LMCD ~ sexoN, data = Hombro, 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.0263 t.test(LMCD ~ sexoN, data = Hombro, 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 ## la alternativa no paramétrica
wilcox.test(LMCD ~ 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: LMCD by sexoN
## W = 43.5, p-value = 3.584e-09
## alternative hypothesis: true location shift is not equal to 0 ## Realizaremos la prueba t para comparar las medias
t.test(C12CD ~ sexoN, data = Hombro, var.equal = TRUE)##
## Two Sample t-test
##
## data: C12CD by sexoN
## t = -9.4101, df = 73, p-value = 3.183e-14
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -6.849137 -4.454997
## sample estimates:
## mean in group Mujer mean in group Hombre
## 30.28626 35.93832 t.test(C12CD ~ sexoN, data = Hombro, var.equal = FALSE)##
## Welch Two Sample t-test
##
## data: C12CD by sexoN
## t = -9.0054, df = 51.263, p-value = 3.866e-12
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -6.911930 -4.392204
## sample estimates:
## mean in group Mujer mean in group Hombre
## 30.28626 35.93832 ## la alternativa no paramétrica
wilcox.test(C12CD ~ 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: C12CD by sexoN
## W = 76, p-value = 1.154e-10
## alternative hypothesis: true location shift is not equal to 0#Resultados Prueba t para LMCD H0: Existe igualdad entre las medias HA: Existen diferencias entre las medias p-value = 2.767e-11, p valor es menor a 0.05,
Se rechaza HO porque en las medias el p valor es <0.05, por lo que existen diferencias significativas entre ambas variables.
#Resultados Prueba t para C12CD H0: Existe igualdad entre las medias HA: Existen diferencias entre las medias p-value = 3.183e-14, p valor es menor a 0.05,
Se rechaza HO porque en las medias el p valor es <0.05, por lo que existen diferencias significativas entre ambas variables.
##Ahora se realizara la prueba de homogeneidad de varianza Igual, primero para la variable LMCD
# Supuesto de homogeneidad de varianza
var.test(LMCD ~ sexoN, data = Hombro)##
## F test to compare two variances
##
## data: LMCD by sexoN
## F = 0.55358, num df = 24, denom df = 34, p-value = 0.1344
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.266830 1.206637
## sample estimates:
## ratio of variances
## 0.5535798Hipotesis: H0: Igualdad de varianzas HA: Diferencia entre las varianzas Se rechaza H0 si p valor < 0.05 p-value LMCD= 0.1344 > 0.05 Por lo que se asumen diferencias entre las varianzas
p-value C12CD= 0.2424 > 0.05 Por lo que se asumen diferencias entre las varianzas
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", "LMCD")])
# Ajustar el modelo con la base depurada
## El modelo discriminante es D=a*LMCD+b
lda1 <- lda(sexoN ~ LMCD, data = Hombro_sinNA)
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[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), " * 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 # Tabla cruzada
tabla_clas <- table(Observado = Hombro_sinNA$sexoN,
Predicho = pred$class)
tabla_clas## Predicho
## Observado Mujer Hombre
## Mujer 22 3
## Hombre 3 32 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 %##Resultados deln modelo con punto de corte ###Función discriminante: D(x) = 0.1089 * LMCD + -15.5587 ###Punto de corte= -0.1791 Si D> -0.1791 es Hombre> Porcentaje de acierto: 90 %
Ejemplo: Supongamos que se encuentra una clavícula Derecha cuya Longitud Máxima es de 119.5mm. Usaremos el modelo para predecir.
prede1 <- predict(lda1, newdata = data.frame(LMCD = 119.5))
prede1$x## LD1
## 1 -2.539385if(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 -2.539385 < -0.1790721 es mujer##Resultado LD1: -2.539385 Debido a que -2.539385 < -0.1790721 es mujer
Ahora 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", "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:
## 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 = Hombro_sinNA$sexoN, Predicho = pred2$class)## Predicho
## Real Mujer Hombre
## Mujer 24 1
## Hombre 1 34 # Porcentaje de acierto
mean(pred2$class == Hombro_sinNA$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), " * LMCD + ", round(a2, 4),"*C12CD +",round(b, 4), "\n")## Función discriminante: D(x) = 0.0727 * LMCD + 0.256 *C12CD + -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 ggplot(Hombro_sinNA, 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
# Ecuación de la recta: a1*LMCD + a2*C12CD + b = 0
# => C12CD = -(a1/a2)*LMCD - b/a2
ggplot(Hombro_sinNA, aes(x = LMCD, y = C12CD, 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")
#Resultados Función discriminante: D(x) = 0.0727 * LMCD + 0.256 *C12CD +
-18.9578 Punto de corte= -0.236 Si D> -0.236 es Hombre Porcentaje de
acierto: 96.66667
Ejemplo: Supongamos que se encuentra una clavícula Derecha cuya Longitud Máxima es de 119.5mm y cuya circunferencia es de 31mm. Usaremos el modelo para predecir.
prede2 <- predict(lda2, newdata = data.frame(LMCD = 119.5, C12CD = 31))
prede2$x## LD1
## 1 -2.333787if(prede2$x>cutoff) cat("Debido a que", prede2$x,">",cutoff, "es hombre") else cat("Debido a que", prede2$x, "<", cutoff,"es mujer")## Debido a que -2.333787 < -0.2359658 es mujer#Resultados LD1: -2.333787, Debido a que -2.333787 < -0.2359658 es mujer,
## Regresión Logística
Aplicaremos la Regresión Logística para determinar la probabilidad de hombre o mujer, usando la clavícula (LMCD y C12CD). Aplicaremos primero la Regresión con una variable (LMCD)
### Ahora aplicaremos Regresión Logística
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMCD")])
# 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_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 3 22
## Hombre 32 3 # Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho) * 100## [1] 90 ## 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) LMCD
## -41.7269684 0.2994501 cat("Ecuación logística:\nlogit(p) = ",
round(coef(modelo1)[1], 4), " + ",
round(coef(modelo1)[2], 4), "*LMCD\n")## Ecuación logística:
## logit(p) = -41.727 + 0.2995 *LMCD##Resultado (Intercept) LMCD -41.7269684 0.2994501 Ecuación logística: logit(p) = -41.727 + 0.2995 * LMCD
Ejemplo: Ejemplo: Supongamos que se encuentra una clavícula Derecha cuya Longitud Máxima es de 123.5mm. Usaremos el modelo para predecir.
prede3 <- predict(modelo1, newdata = data.frame(LMCD = 123.5))
prede3## 1
## -4.744881if(prede3>cutoff) cat("Debido a que", prede3,">",cutoff, "es hombre") else cat("Debido a que", prede3, "<", cutoff,"es mujer")## Debido a que -4.744881 < -0.2359658 es mujer##Resultados Debido a que 0.008621123 > -0.2359658, es hombre
Ahora vamos a ajustar un Discriminante lineal con 2 variables
### Ahora aplicaremos Regresión Logística
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMCD","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##Resultados Ecuación logística: logit(p) = -69.0885 + 0.3067 LMCD + 0.8022 C12CD
Ejemplo: Ejemplo: Supongamos que se encuentra una clavícula Derecha cuya Longitud Máxima es de 123.5mm y cuya circunferencia es de 35mm. Usaremos el modelo para predecir.
prede4 <- predict(modelo2,
newdata = data.frame(LMCD = 123.5, C12CD = 35),
type = "response") # ESTO ES IMPORTANTE
if(prede4>cutoff) cat("Debido a que", prede4,">",
cutoff, "es hombre") else cat("Debido a que",
prede4, "<",
cutoff,"es mujer")## Debido a que 0.04191591 > -0.2359658 es hombre##Resultados Debido a que 0.04191591 > -0.2359658 es hombre