Punto 1: Cargar datos y explorar variables
############################
# Punto 1: Cargar datos #
############################
Hombro <- haven::read_sav("Datos hombro.sav")
# Leemos el archivo SPSS y lo guardamos en el objeto Hombro
str(Hombro)
## tibble [80 × 16] (S3: tbl_df/tbl/data.frame)
## $ FOLIO: chr [1:80] "028" "042" "053" "054" ...
## ..- attr(*, "format.spss")= chr "A3"
## ..- attr(*, "display_width")= int 9
## $ sexoN: dbl+lbl [1:80] 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## ..@ label : chr "Sexo"
## ..@ format.spss: chr "F8.0"
## ..@ labels : Named num [1:2] 1 2
## .. ..- attr(*, "names")= chr [1:2] "Hombre" "Mujer"
## $ LMHD : num [1:80] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "Longitud máxima húmero derecho"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ LMHI : num [1:80] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "Longitud máxima húmero izquierdo"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ ABHD : num [1:80] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "Altura biomecánica húmero derecho"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ ABHI : num [1:80] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "label")= chr "Altura biomecánica húmero izquierdo"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ LMCD : num [1:80] 148 158 NA NA NA NA NA NA 135 NA ...
## ..- attr(*, "label")= chr "Longitud máxima clavícula derecha"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ LMCI : num [1:80] 152 159 NA NA NA ...
## ..- attr(*, "label")= chr "Longitud máxima clavícula izquierdo"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ C12CD: num [1:80] 34.4 35 37 33 NA ...
## ..- attr(*, "label")= chr "Circunferencia en punto medio clavícula derecha"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
## $ C12CI: num [1:80] 36 35 35 35 NA ...
## ..- attr(*, "label")= chr "Circunferencia en punto medio clavícula izquierdo"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
## $ DAPCD: num [1:80] 11.5 13 11 12 NA 12.5 NA 11 12 10 ...
## ..- attr(*, "label")= chr "Diámetro antero-posterior clavícula derecha"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
## $ DAPCI: num [1:80] 13 12.5 12 12 NA 14 12 NA 11 NA ...
## ..- attr(*, "label")= chr "Diámetro antero-posterior clavícula izquierdo"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
## $ AlEED: num [1:80] NA NA 146 NA NA ...
## ..- attr(*, "label")= chr "Altura escapular escápula derecha"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ AlEEI: num [1:80] NA NA 147 NA NA ...
## ..- attr(*, "label")= chr "Altura escapular escápula izquierdo"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ AnEED: num [1:80] 103 NA 103 NA 101 ...
## ..- attr(*, "label")= chr "Anchura escapular escápula derecha"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
## $ AnEEI: num [1:80] 104 NA 104 NA 101 ...
## ..- attr(*, "label")= chr "Anchura escapular escápula izquierdo"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
# Mostramos la estructura de la base: tipo de objeto, número de casos y variables
names(Hombro)
## [1] "FOLIO" "sexoN" "LMHD" "LMHI" "ABHD" "ABHI" "LMCD" "LMCI" "C12CD"
## [10] "C12CI" "DAPCD" "DAPCI" "AlEED" "AlEEI" "AnEED" "AnEEI"
# Mostramos los nombres de todas las variables de Hombro
Punto 2: Transformar sexo y construir Hombro2
###########################################################
# Punto 2: Revisar y transformar la variable de sexo (sexoN)
###########################################################
table(Hombro$sexoN, useNA = "ifany")
##
## 1 2
## 50 30
# 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")
##
## Hombre Mujer
## 50 30
# Verificamos cómo quedan las frecuencias después de convertir a factor
levels(Hombro$sexoN)
## [1] "Hombre" "Mujer"
# 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)
## [1] "Mujer" "Hombre"
# Confirmamos el nuevo orden de las categorías de sexoN
######################################################
# Punto 3: Construir subconjunto con LMCD y C12CD #
######################################################
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)
## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
## $ sexoN: Factor w/ 2 levels "Mujer","Hombre": 2 2 1 1 1 1 1 1 1 1 ...
## $ LMCD : num [1:60] 148 158 135 125 116 ...
## ..- attr(*, "label")= chr "Longitud máxima clavícula derecha"
## ..- attr(*, "format.spss")= chr "F5.1"
## ..- attr(*, "display_width")= int 12
## $ C12CD: num [1:60] 34.4 35 33 28.4 30 ...
## ..- attr(*, "label")= chr "Circunferencia en punto medio clavícula derecha"
## ..- attr(*, "format.spss")= chr "F4.1"
## ..- attr(*, "display_width")= int 12
# Comprobamos cuántas observaciones y variables tiene el nuevo objeto Hombro2
table(Hombro2$sexoN)
##
## Mujer Hombre
## 25 35
# Vemos cuántos hombres y mujeres quedan en el subconjunto Hombro2
Punto 3: Comparaciones por sexo (LMCD y C12CD)
########################################
# 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("\nVariable:", var_nombre, "\n")
cat("Shapiro H =", round(p_H,4),
" / Shapiro M =", round(p_M,4), "\n")
if (p_H > 0.05 & p_M > 0.05) {
cat("Usamos t de Student\n")
test <- t.test(x ~ sexoN, data = Hombro2)
} else {
cat("Usamos Mann-Whitney\n")
test <- wilcox.test(x ~ sexoN, data = Hombro2)
}
print(test)
}
comparar_por_sexo("LMCD")
##
## Variable: LMCD
## Shapiro H = 0.1055 / Shapiro M = 0.8169
## Usamos t de Student
##
## Welch Two Sample t-test
##
## data: x 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
# Comparamos LMCD entre hombres y mujeres usando la función comparar_por_sexo
comparar_por_sexo("C12CD")
##
## Variable: C12CD
## Shapiro H = 0.2534 / Shapiro M = 0.2948
## Usamos t de Student
##
## Welch Two Sample t-test
##
## data: x by sexoN
## t = -7.8461, df = 45.117, p-value = 5.68e-10
## alternative hypothesis: true difference in means between group Mujer and group Hombre is not equal to 0
## 95 percent confidence interval:
## -6.860803 -4.058095
## sample estimates:
## mean in group Mujer mean in group Hombre
## 30.31206 35.77151
# Comparamos C12CD entre hombres y mujeres con la misma función
Punto 4: LDA con solo 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
## Call:
## lda(sexoN ~ LMCD, data = Hombro2)
##
## Prior probabilities of groups:
## Mujer Hombre
## 0.4166667 0.5833333
##
## Group means:
## LMCD
## Mujer 131.3026
## Hombre 151.0263
##
## Coefficients of linear discriminants:
## LD1
## LMCD 0.1089485
# Mostramos medias por grupo y coeficientes del modelo discriminante
pred1 <- predict(disc1)
# Obtenemos clasificaciones y probabilidades posteriores para cada individuo
pred1
## $class
## [1] Hombre Hombre Mujer Mujer Mujer Mujer Mujer Mujer Mujer Mujer
## [11] Mujer Mujer Hombre Mujer Hombre Mujer Mujer Mujer Mujer Mujer
## [21] Mujer Mujer Mujer Hombre Mujer Hombre Mujer Mujer Hombre Hombre
## [31] Hombre Hombre Mujer Hombre Hombre Hombre Hombre Mujer Hombre Hombre
## [41] Hombre Hombre Hombre Hombre Hombre Hombre Mujer Hombre Hombre Hombre
## [51] Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre
## Levels: Mujer Hombre
##
## $posterior
## Mujer Hombre
## 1 0.1260029038 0.873997096
## 2 0.0136815848 0.986318415
## 3 0.7515192045 0.248480795
## 4 0.9712929261 0.028707074
## 5 0.9956905066 0.004309493
## 6 0.8966228697 0.103377130
## 7 0.9861741417 0.013825858
## 8 0.9951350049 0.004864995
## 9 0.9488811825 0.051118817
## 10 0.8852580772 0.114741923
## 11 0.6544142881 0.345585712
## 12 0.8284898988 0.171510101
## 13 0.2688805653 0.731119435
## 14 0.9567570026 0.043242997
## 15 0.3700322955 0.629967704
## 16 0.9105832491 0.089416751
## 17 0.9876826201 0.012317380
## 18 0.7926291698 0.207370830
## 19 0.9427330877 0.057266912
## 20 0.9427330877 0.057266912
## 21 0.9069802971 0.093019703
## 22 0.6803884961 0.319611504
## 23 0.8284898988 0.171510101
## 24 0.4260556107 0.573944389
## 25 0.5997405749 0.400259425
## 26 0.3173030151 0.682696985
## 27 0.6544142881 0.345585712
## 28 0.7515192045 0.248480795
## 29 0.0215416906 0.978458309
## 30 0.4260556107 0.573944389
## 31 0.1871645030 0.812835497
## 32 0.1023954120 0.897604588
## 33 0.5424646191 0.457535381
## 34 0.1260029038 0.873997096
## 35 0.1260029038 0.873997096
## 36 0.3633345989 0.636665401
## 37 0.0341757375 0.965824263
## 38 0.8444896642 0.155510336
## 39 0.0125723447 0.987427655
## 40 0.1871645030 0.812835497
## 41 0.0428050787 0.957194921
## 42 0.0827919587 0.917208041
## 43 0.1260029038 0.873997096
## 44 0.2464999368 0.753500063
## 45 0.0281356116 0.971864388
## 46 0.0002305177 0.999769482
## 47 0.8284898988 0.171510101
## 48 0.2384359256 0.761564074
## 49 0.0341757375 0.965824263
## 50 0.0428050787 0.957194921
## 51 0.0827919587 0.917208041
## 52 0.0243009347 0.975699065
## 53 0.3704483262 0.629551674
## 54 0.0341757375 0.965824263
## 55 0.1023954120 0.897604588
## 56 0.0003681246 0.999631875
## 57 0.0108568362 0.989143164
## 58 0.0796619614 0.920338039
## 59 0.0006607779 0.999339222
## 60 0.0013328728 0.998667127
##
## $x
## LD1
## 1 0.56564625
## 2 1.65513096
## 3 -0.85068387
## 4 -1.97441778
## 5 -2.86843900
## 6 -1.34095199
## 7 -2.32148823
## 8 -2.81175635
## 9 -1.69503615
## 10 -1.28647776
## 11 -0.63278693
## 12 -1.06858081
## 13 0.12985237
## 14 -1.77674587
## 15 -0.08804458
## 16 -1.41565452
## 17 -2.37596246
## 18 -0.95963234
## 19 -1.63915983
## 20 -1.63915983
## 21 -1.39542623
## 22 -0.68726117
## 23 -1.06858081
## 24 -0.19699305
## 25 -0.52383846
## 26 0.02090390
## 27 -0.63278693
## 28 -0.85068387
## 29 1.44016152
## 30 -0.19699305
## 31 0.34774931
## 32 0.67459472
## 33 -0.41488999
## 34 0.56564625
## 35 0.56564625
## 36 -0.07462270
## 37 1.21933708
## 38 -1.12305505
## 39 1.69500093
## 40 0.34774931
## 41 1.11038861
## 42 0.78354319
## 43 0.56564625
## 44 0.18432660
## 45 1.31274182
## 46 3.56172920
## 47 -1.06858081
## 48 0.20475893
## 49 1.21933708
## 50 1.11038861
## 51 0.78354319
## 52 1.38275978
## 53 -0.08887492
## 54 1.21933708
## 55 0.67459472
## 56 3.34383226
## 57 1.76407943
## 58 0.80306301
## 59 3.07146108
## 60 2.74461567
# 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
## Predicho
## Real Mujer Hombre
## Mujer 22 3
## Hombre 3 32
# 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
## [1] 90
# 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
## $class
## [1] Hombre
## Levels: Mujer Hombre
##
## $posterior
## Mujer Hombre
## 1 0.4840437 0.5159563
##
## $x
## LD1
## 1 -0.3059415
# Mostramos sexo asignado y probabilidades posteriores del ejemplo
Punto 5: LDA con LMCD y C12CD
#########################################
# 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
## Call:
## lda(sexoN ~ LMCD + C12CD, data = Hombro2)
##
## 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
# Inspeccionamos coeficientes y medias por grupo del modelo con dos predictores
pred2 <- predict(disc2)
# Calculamos clasificaciones y probabilidades posteriores para todos los casos
pred2
## $class
## [1] Hombre Hombre Mujer Mujer Mujer Mujer Mujer Mujer Mujer Mujer
## [11] Mujer Mujer Hombre Mujer Mujer Mujer Mujer Mujer Mujer Mujer
## [21] Mujer Mujer Mujer Hombre Mujer Hombre Mujer Hombre Hombre Hombre
## [31] Hombre Hombre Mujer Hombre Hombre Mujer Hombre Mujer Hombre Hombre
## [41] Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre
## [51] Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre Hombre
## Levels: Mujer Hombre
##
## $posterior
## Mujer Hombre
## 1 6.132754e-02 0.9386724621
## 2 5.368094e-03 0.9946319056
## 3 7.237080e-01 0.2762919870
## 4 9.983476e-01 0.0016524175
## 5 9.990426e-01 0.0009573886
## 6 9.990549e-01 0.0009450619
## 7 9.886497e-01 0.0113503097
## 8 9.997498e-01 0.0002502171
## 9 9.560486e-01 0.0439514311
## 10 9.908592e-01 0.0091407631
## 11 9.848640e-01 0.0151360139
## 12 8.908575e-01 0.1091425370
## 13 4.585656e-01 0.5414343566
## 14 9.995850e-01 0.0004150006
## 15 9.709280e-01 0.0290720472
## 16 9.811070e-01 0.0188930072
## 17 9.988241e-01 0.0011758575
## 18 9.960018e-01 0.0039982119
## 19 9.067925e-01 0.0932074839
## 20 5.691369e-01 0.4308631185
## 21 9.380314e-01 0.0619686090
## 22 9.953473e-01 0.0046527259
## 23 7.334613e-01 0.2665386703
## 24 7.966774e-02 0.9203322634
## 25 5.854784e-01 0.4145216371
## 26 3.353530e-01 0.6646469552
## 27 7.817730e-01 0.2182269746
## 28 1.596810e-01 0.8403190133
## 29 2.498618e-03 0.9975013818
## 30 4.323372e-01 0.5676627841
## 31 2.138794e-01 0.7861205970
## 32 3.327774e-02 0.9667222606
## 33 9.773277e-01 0.0226723436
## 34 4.344422e-01 0.5655577659
## 35 2.007508e-02 0.9799249212
## 36 8.417076e-01 0.1582923998
## 37 6.765708e-04 0.9993234292
## 38 7.530960e-01 0.2469039912
## 39 3.560689e-03 0.9964393110
## 40 6.000781e-02 0.9399921885
## 41 1.488325e-02 0.9851167468
## 42 2.830352e-03 0.9971696477
## 43 4.783939e-03 0.9952160615
## 44 4.763575e-03 0.9952364246
## 45 5.973455e-04 0.9994026545
## 46 8.096146e-06 0.9999919039
## 47 3.100537e-01 0.6899462701
## 48 7.718711e-02 0.9228128918
## 49 2.877081e-03 0.9971229187
## 50 7.265162e-03 0.9927348379
## 51 1.339026e-02 0.9866097436
## 52 8.949111e-03 0.9910508885
## 53 9.108027e-04 0.9990891973
## 54 1.214763e-02 0.9878523651
## 55 1.640106e-02 0.9835989359
## 56 2.219213e-04 0.9997780787
## 57 1.837757e-02 0.9816224338
## 58 4.261755e-02 0.9573824498
## 59 3.712447e-04 0.9996287553
## 60 3.791787e-05 0.9999620821
##
## $x
## LD1
## 1 0.60870665
## 2 1.48936478
## 3 -0.69486072
## 4 -2.61637258
## 5 -2.80936708
## 6 -2.81394801
## 7 -1.93238668
## 8 -3.28351358
## 9 -1.44242456
## 10 -2.00963432
## 11 -1.82938149
## 12 -1.09626073
## 13 -0.29612770
## 14 -3.10477443
## 15 -1.59384186
## 16 -1.74973118
## 17 -2.73670037
## 18 -2.30348809
## 19 -1.15825980
## 20 -0.45308863
## 21 -1.31438054
## 22 -2.24971498
## 23 -0.71228054
## 24 0.50933929
## 25 -0.47674091
## 26 -0.11320788
## 27 -0.80543431
## 28 0.23166256
## 29 1.76045493
## 30 -0.25862109
## 31 0.10491193
## 32 0.83500534
## 33 -1.68396828
## 34 -0.26164845
## 35 1.01828553
## 36 -0.94491780
## 37 2.22248554
## 38 -0.74863384
## 39 1.63498421
## 40 0.61688552
## 41 1.12583176
## 42 1.71631153
## 43 1.53025911
## 44 1.53177279
## 45 2.26649642
## 46 3.78567755
## 47 -0.07231355
## 48 0.52146110
## 49 1.71051195
## 50 1.38181855
## 51 1.16369874
## 52 1.30759827
## 53 2.11741253
## 54 1.19853836
## 55 1.09099213
## 56 2.61631716
## 57 1.05009780
## 58 0.74421330
## 59 2.43455064
## 60 3.24037801
# 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
## Predicho
## Real Mujer Hombre
## Mujer 24 1
## Hombre 1 34
# 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
## [1] 96.66667
# 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
## $class
## [1] Hombre
## Levels: Mujer Hombre
##
## $posterior
## Mujer Hombre
## 1 0.005821304 0.9941787
##
## $x
## LD1
## 1 1.46058
# Mostramos sexo asignado y probabilidades posteriores del caso hipotético
Punto 6: Modelo logístico con solo 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)
##
## Call:
## glm(formula = sexoN ~ LMCD, family = binomial(link = "logit"),
## data = Hombro2)
##
## 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
# 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
## Predicho
## Real Mujer Hombre
## Mujer 22 3
## Hombre 3 32
# 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
## [1] 90
# 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
## 1
## 0.960477
# Mostramos la probabilidad numérica estimada
Gráfica de probabilidades predichas (modelo LMCD)
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
Punto 7: Modelo logístico con LMCD + C12CD
#############################################
# 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)
##
## Call:
## glm(formula = sexoN ~ LMCD + C12CD, family = binomial(link = "logit"),
## data = Hombro2)
##
## 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
# 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
## Predicho
## Real Mujer Hombre
## Mujer 24 1
## Hombre 1 34
# 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
## [1] 96.66667
# 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
## 1
## 0.9993926
# Mostramos la probabilidad numérica asociada al caso con ambas medidas