Introducción

#en este ejercicio ajustaremos modelos de regresión logistica que nos permitan estimar el sexo de un individuo a traves de los huesos considerando 2 variables LMDH Y ABDH https://rpubs.com/YumeWalker/1365227

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

## Defino mi directorio de trabajo
setwd("~/Itza")
##Abriendo paquete pacman
library(pacman)
## El archivo esta en formato SPSS, lo abrimos mediante la libreria haven
p_load(haven,dplyr,ggplot2,MASS,tinytex)
Hombro <- read_sav("Datos hombro.sav")
table (Hombro$sexoN)  ## Frecuencias de sexo
## 
##  1  2 
## 50 30
## Cambiamos las categorías para predecir p(Hombre)
Hombro$sexoN <- factor(Hombro$sexoN,
                       levels = c(2, 1),
                       labels = c("Mujer", "Hombre"))
table (Hombro$sexoN)  ## Frecuencias de sexo
## 
##  Mujer Hombre 
##     30     50
## Omitir datos perdidos
Hombro_sinNA <- na.omit(Hombro[, c("sexoN", "LMHD","ABHD")])
# Ajustar modelo logístico binario
modelo1 <- glm(sexoN ~ LMHD + ABHD,
                  data = Hombro_sinNA,
                  family = binomial(link = "logit"))
summary(modelo1)
## 
## 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(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     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(modelo1)
## (Intercept)        LMHD        ABHD 
## -50.0490482  -0.4773712   0.6555261
cat("Ecuación logística:\nlogit(p) = ",
    round(coef(modelo1)[1], 4), " + ",
    round(coef(modelo1)[2], 4), "*LMHD + ",
    round(coef(modelo1)[3], 4), "*ABHD\n")
## Ecuación logística:
## logit(p) =  -50.049  +  -0.4774 *LMHD +  0.6555 *ABHD
modelo2 <- glm(sexoN ~ ABHD,
               data = Hombro_sinNA,
               family = binomial(link = "logit"))
summary(modelo2)
## 
## Call:
## glm(formula = sexoN ~ ABHD, family = binomial(link = "logit"), 
##     data = Hombro_sinNA)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -44.80276   10.96751  -4.085 4.41e-05 ***
## ABHD          0.15406    0.03746   4.112 3.92e-05 ***
## ---
## 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: 42.819  on 65  degrees of freedom
## AIC: 46.819
## 
## Number of Fisher Scoring iterations: 6
# Probabilidades de ser "Hombre"
Hombro_sinNA$prob_Hombre2 <- predict(modelo2, type = "response")

# Clasificación usando 0.5 como punto de corte
Hombro_sinNA$predicho2 <- ifelse(Hombro_sinNA$prob_Hombre2 >= 0.5, "Hombre", "Mujer")

# Matriz de confusión
table(Real = Hombro_sinNA$sexoN, Predicho2 = Hombro_sinNA$predicho2)
##         Predicho2
## Real     Hombre Mujer
##   Mujer       5    21
##   Hombre     36     5
# Porcentaje de acierto
mean(Hombro_sinNA$sexoN == Hombro_sinNA$predicho2) * 100
## [1] 85.07463
## 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        ABHD 
## -50.0490482  -0.4773712   0.6555261
cat("Ecuación logística:\nlogit(p) = ",
    round(coef(modelo1)[1], 4), " + ",
    round(coef(modelo1)[2], 4), "*LMHD + ",
    round(coef(modelo1)[3], 4), "*ABHD\n")
## Ecuación logística:
## logit(p) =  -50.049  +  -0.4774 *LMHD +  0.6555 *ABHD
summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.