{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, # mostrar código message = FALSE, # ocultar mensajes de paquetes warning = FALSE # ocultar warnings )
library(haven) library(dplyr) library(ggplot2) library(MASS)
{r cargar-datos} ############################ # Punto 1: Cargar datos # ############################
Hombro <- haven::read_sav(“Datos hombro.sav”) # Leemos el archivo SPSS y lo guardamos en el objeto Hombro
str(Hombro) # Mostramos la estructura de la base: tipo de objeto, número de casos y variables
names(Hombro) # Mostramos los nombres de todas las variables de Hombro
{r preparar-hombro2} ########################################################### # Punto 2: Revisar y transformar la variable de sexo (sexoN) ###########################################################
table(Hombro$sexoN, useNA = “ifany”) # Contamos cuántos códigos 1 y 2 hay en sexoN y verificamos si existen NA
Hombro\(sexoN <- haven::as_factor(Hombro\)sexoN) # Convertimos la variable etiquetada de SPSS en un factor de R con etiquetas
table(Hombro$sexoN, useNA = “ifany”) # Verificamos cómo quedan las frecuencias después de convertir a factor
levels(Hombro$sexoN) # Mostramos el orden actual de las categorías de sexoN
Hombro\(sexoN <- relevel(Hombro\)sexoN, ref = “Mujer”) # Fijamos “Mujer” como categoría de referencia en el factor sexoN
levels(Hombro$sexoN) # Confirmamos el nuevo orden de las categorías de sexoN
Hombro2 <- Hombro %>% dplyr::select(sexoN, LMCD, C12CD) %>% # Nos quedamos solo con sexoN, LMCD y C12CD filter(!is.na(sexoN), # Retenemos solo casos con sexoN no faltante !is.na(LMCD), # Retenemos solo casos con LMCD no faltante !is.na(C12CD)) # Retenemos solo casos con C12CD no faltante
str(Hombro2) # Comprobamos cuántas observaciones y variables tiene el nuevo objeto Hombro2
table(Hombro2$sexoN) # Vemos cuántos hombres y mujeres quedan en el subconjunto Hombro2
{r comparar-por-sexo} ######################################## # Punto 4: Comparaciones por sexo # ########################################
comparar_por_sexo <- function(var_nombre) { x <- Hombro2[[var_nombre]] grupo_H <- x[Hombro2\(sexoN == "Hombre"] grupo_M <- x[Hombro2\)sexoN == “Mujer”]
p_H <- shapiro.test(grupo_H)\(p.value p_M <- shapiro.test(grupo_M)\)p.value
cat(“:”, var_nombre, “”) cat(“Shapiro H =”, round(p_H,4), ” / Shapiro M =“, round(p_M,4),”“)
if (p_H > 0.05 & p_M > 0.05) { cat(“Usamos t de Student”) test <- t.test(x ~ sexoN, data = Hombro2) } else { cat(“Usamos Mann-Whitney”) test <- wilcox.test(x ~ sexoN, data = Hombro2) } print(test) }
comparar_por_sexo(“LMCD”) # Comparamos LMCD entre hombres y mujeres usando la función comparar_por_sexo
comparar_por_sexo(“C12CD”) # Comparamos C12CD entre hombres y mujeres con la misma función
{r lda-lmcd} ############################# # Punto 5: LDA con solo LMCD #############################
disc1 <- lda(sexoN ~ LMCD, data = Hombro2) # Ajustamos una función discriminante lineal donde LMCD clasifica el sexo
disc1 # Mostramos medias por grupo y coeficientes del modelo discriminante
pred1 <- predict(disc1) # Obtenemos clasificaciones y probabilidades posteriores para cada individuo
pred1 # Mostramos las predicciones (clase y posterior probabilities)
tab_disc1 <- table(Real = Hombro2\(sexoN, Predicho = pred1\)class) # Construimos la matriz de confusión entre sexo real y sexo predicho
tab_disc1 # Visualizamos la matriz de confusión
acc_disc1 <- mean(pred1\(class == Hombro2\)sexoN) * 100 # Calculamos el porcentaje de casos correctamente clasificados por el modelo
acc_disc1 # Mostramos el porcentaje de clasificación correcta del modelo con LMCD
nuevo_LMCD <- data.frame(LMCD = 140) # Definimos un ejemplo hipotético con LMCD igual a 140
pred_nuevo1 <- predict(disc1, nuevo_LMCD) # Obtenemos la clasificación y probabilidades para el caso hipotético
pred_nuevo1 # Mostramos sexo asignado y probabilidades posteriores del ejemplo
{r lda-dos-variables} ######################################### # Punto 6: LDA con LMCD y C12CD # #########################################
disc2 <- lda(sexoN ~ LMCD + C12CD, data = Hombro2) # Ajustamos una función discriminante con las dos variables LMCD y C12CD
disc2 # Inspeccionamos coeficientes y medias por grupo del modelo con dos predictores
pred2 <- predict(disc2) # Calculamos clasificaciones y probabilidades posteriores para todos los casos
pred2 # Mostramos el objeto de predicciones (clase y probabilidades)
tab_disc2 <- table(Real = Hombro2\(sexoN, Predicho = pred2\)class) # Construimos la matriz de confusión del modelo con LMCD y C12CD
tab_disc2 # Visualizamos la matriz de confusión del segundo modelo
acc_disc2 <- mean(pred2\(class == Hombro2\)sexoN) * 100 # Obtenemos el porcentaje de clasificación correcta del modelo con dos variables
acc_disc2 # Mostramos el porcentaje de clasificación correcta
nuevo_2 <- data.frame( LMCD = 140, # Valor hipotético de LMCD C12CD = 40 # Valor hipotético de C12CD ) # Definimos un caso hipotético con ambas medidas
pred_nuevo2 <- predict(disc2, nuevo_2) # Obtenemos la clasificación discriminante para el individuo hipotético
pred_nuevo2 # Mostramos sexo asignado y probabilidades posteriores del caso hipotético
{r logistico-lmcd} ############################################# # Punto 7: Modelo logístico con solo LMCD # #############################################
mod_log1 <- glm(sexoN ~ LMCD, data = Hombro2, family = binomial(link = “logit”)) # Ajustamos un modelo de regresión logística donde LMCD predice la probabilidad de ser Hombre
summary(mod_log1) # Revisamos coeficientes, errores estándar y significancia del modelo logístico
Hombro2$prob_Hombre1 <- predict(mod_log1, type = “response”) # Calculamos la probabilidad estimada de ser Hombre para cada individuo
Hombro2\(predicho_log1 <- ifelse(Hombro2\)prob_Hombre1 >= 0.5, “Hombre”, “Mujer”) # Asignamos la etiqueta Hombre o Mujer usando un punto de corte de 0.5
Hombro2\(predicho_log1 <- factor(Hombro2\)predicho_log1, levels = levels(Hombro2$sexoN)) # Convertimos la clasificación en factor con los mismos niveles que el sexo real
tab_log1 <- table(Real = Hombro2\(sexoN, Predicho = Hombro2\)predicho_log1) # Construimos la matriz de confusión del modelo logístico con LMCD
tab_log1 # Visualizamos la matriz de confusión
acc_log1 <- mean(Hombro2\(sexoN == Hombro2\)predicho_log1) * 100 # Calculamos el porcentaje de acierto del modelo logístico
acc_log1 # Mostramos el porcentaje de casos correctamente clasificados
nuevo_LMCD <- data.frame(LMCD = 150) # Definimos un nuevo valor hipotético de LMCD en 150
prob_hombre_LMCD <- predict(mod_log1, nuevo_LMCD, type = “response”) # Obtenemos la probabilidad de ser Hombre para ese valor hipotético
prob_hombre_LMCD # Mostramos la probabilidad numérica estimada
{r grafica-log1} ggplot(Hombro2, aes(x = prob_Hombre1, fill = sexoN)) + geom_density(alpha = 0.4) + geom_vline(xintercept = 0.5, linetype = “dashed”) + theme_minimal() + labs(title = “Probabilidades predichas de ser Hombre (modelo LMCD)”, x = “P(Hombre)”, y = “Densidad”) # Graficamos las distribuciones de P(Hombre) por sexo real # y observamos cómo funciona el punto de corte 0.5 en la separación de grupos
{r logistico-dos-variables} ############################################# # Punto 8: Modelo logístico con LMCD + C12CD #############################################
mod_log2 <- glm(sexoN ~ LMCD + C12CD, data = Hombro2, family = binomial(link = “logit”)) # Ajustamos un modelo logístico donde LMCD y C12CD predicen conjuntamente el sexo
summary(mod_log2) # Revisamos la significancia y los coeficientes del modelo con dos predictores
Hombro2$prob_Hombre2 <- predict(mod_log2, type = “response”) # Calculamos la probabilidad de ser Hombre para cada caso usando ambos predictores
Hombro2\(predicho_log2 <- ifelse(Hombro2\)prob_Hombre2 >= 0.5, “Hombre”, “Mujer”) # Asignamos la clasificación Hombre/Mujer según el punto de corte 0.5
Hombro2\(predicho_log2 <- factor(Hombro2\)predicho_log2, levels = levels(Hombro2$sexoN)) # Garantizamos que la variable predicha tenga los mismos niveles que sexoN
tab_log2 <- table(Real = Hombro2\(sexoN, Predicho = Hombro2\)predicho_log2) # Armamos la matriz de confusión del modelo logístico con dos variables
tab_log2 # Visualizamos la matriz de confusión resultante
acc_log2 <- mean(Hombro2\(sexoN == Hombro2\)predicho_log2) * 100 # Calculamos el porcentaje de acierto del modelo logístico con LMCD y C12CD
acc_log2 # Mostramos el valor final de ese porcentaje
nuevo_vals <- data.frame( LMCD = 150, C12CD = 38 ) # Creamos un caso hipotético con LMCD = 150 y C12CD = 38
prob_hombre_2 <- predict(mod_log2, nuevo_vals, type = “response”) # Obtenemos la probabilidad estimada de ser Hombre para ese individuo hipotético
prob_hombre_2 # Mostramos la probabilidad numérica asociada al caso con ambas medidas