1)              Un estudio quiere establecer un modelo que permita calcular la probabilidad de obtener una matrícula de honor al final del bachillerato en función de la nota que se ha obtenido en matemáticas. La variable matrícula está codificada como 0 si no se tiene matrícula y 1 si se tiene.

 

a)   Genere la ecuación logística. Interprete los odds. Grafique el modelo.

b)   Realice la prueba estadística para saber si haydiferencia de los residuos del modelo y el modelo nulo.

c)    A un nivel de confianza es del 95%, ¿cuál es el intervalo de confianza?

d)   compare de clasificación predicha y observaciones.

matricula <- as.factor(c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1,
1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1,
0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1)) 
matematicas <- c(41, 53, 54, 47, 57, 51,
42, 45, 54, 52, 51, 51, 71, 57, 50, 43, 51, 60, 62, 57, 35, 75, 45, 57, 45, 46,
66, 57, 49, 49, 57, 64, 63, 57, 50, 58, 75, 68, 44, 40, 41, 62, 57, 43, 48, 63,
39, 70, 63, 59, 61, 38, 61, 49, 73, 44, 42, 39, 55, 52, 45, 61, 39, 41, 50, 40,
60, 47, 59, 49, 46, 58, 71, 58, 46, 43, 54, 56, 46, 54, 57, 54, 71, 48, 40, 64,
51, 39, 40, 61, 66, 49, 65, 52, 46, 61, 72, 71, 40, 69, 64, 56, 49, 54, 53, 66,
67, 40, 46, 69, 40, 41, 57, 58, 57, 37, 55, 62, 64, 40, 50, 46, 53, 52, 45, 56,
45, 54, 56, 41, 54, 72, 56, 47, 49, 60, 54, 55, 33, 49, 43, 50, 52, 48, 58, 43,
41, 43, 46, 44, 43, 61, 40, 49, 56, 61, 50, 51, 42, 67, 53, 50, 51, 72, 48, 40,
53, 39, 63, 51, 45, 39, 42, 62, 44, 65, 63, 54, 45, 60, 49, 48, 57, 55, 66, 64,
55, 42, 56, 53, 41, 42, 53, 42, 60, 52, 38, 57, 58, 65) 
datos<-data.frame(matricula, matematicas)
head(datos,4)
##   matricula matematicas
## 1         0          41
## 2         0          53
## 3         0          54
## 4         0          47

Representación de las observaciones

library(ggplot2) 
## Warning: package 'ggplot2' was built under R version 4.2.2
table(datos$matricula)
## 
##   0   1 
## 151  49
ggplot(data = datos, aes(x = matricula, y = matematicas, color = matricula)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.1) + theme_bw() + theme(legend.position = "null")

Parece existir una diferencia entre la nota de las personas con matrícula y sin matrícula.

modelo <- glm(matricula ~ matematicas, data = datos, family = "binomial")
summary(modelo)
## 
## Call:
## glm(formula = matricula ~ matematicas, family = "binomial", data = datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0332  -0.6785  -0.3506  -0.1565   2.6143  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.79394    1.48174  -6.610 3.85e-11 ***
## matematicas  0.15634    0.02561   6.105 1.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 222.71  on 199  degrees of freedom
## Residual deviance: 167.07  on 198  degrees of freedom
## AIC: 171.07
## 
## Number of Fisher Scoring iterations: 5
#para hallar los odd
#El coeficiente estimado para la intersección es el valor esperado del logaritmo de odds de que un estudiante obtenga matrícula teniendo un 0 en en matemáticas. Como es de esperar, los odds son muy bajos
exp(-9.793394)
## [1] 5.581913e-05

lo que se corresponde con una probabilidad de obtener matrícula de 0.00

Acorde al modelo, el logaritmo de los odds de que un estudiante tenga matrícula está positivamente relacionado con la puntuación obtenida en matemáticas (coeficiente de regresión = 0.1563404). Esto significa que, por cada unidad que se incrementa la variable matemáticas, se espera que el logaritmo de odds de la variable matrícula se incremente en promedio 0.1563404 unidades.

exp(0.15634)
## [1] 1.169224

se obtiene que, por cada unidad que se incrementa la variable matemáticas, los odds de obtener matrícula se incremente en promedio 1.169 unidades. No hay que confundir esto último con que la probabilidad de matrícula se incremente un 1.169 %.

Además del valor de las estimaciones de los coeficientes parciales de correlación del modelo, es conveniente calcular sus correspondientes intervalos de confianza.

confint(object = modelo, level = 0.95 )
## Waiting for profiling to be done...
##                   2.5 %     97.5 %
## (Intercept) -12.9375208 -7.0938806
## matematicas   0.1093783  0.2103937

Gráfico del modelo

Dado que un modelo logístico modela el logaritmo de ODDs, estas son las unidades en las que se devuelven las predicciones. Es necesario convertirlas de nuevo en probabilidad mediante la función logit. En R, la función predict() puede devolver directamente las probabilidades en lugar de los logODDs si se indica el argumento type=“response”. Sin embargo, si se quieren calcular intervalos de confianza y que estos no se salgan del rango [0, 1] es necesario emplear los logODDs y una vez que se les ha sustraído o sumado el margen de error (Z x SE) se transforman en probabilidades.

# MEDIANTE BASE GRAPHICS SIN INTERVALOS DE CONFIANZA 
# Codificación 0,1 de la variable respuesta 
datos$matricula <- as.character(datos$matricula) 
datos$matricula <- as.numeric(datos$matricula) 
plot(matricula ~ matematicas, datos, col = "darkblue", main = "Modelo regresión logística", ylab = "P(matrícula=1|matemáticas)", xlab = "matemáticas", pch = "I") 

#type="response" devuelve las predicciones en forma de probabilidad en lugar de en #log_ODDs
#curve(predict(modelo, data.frame(matematicas = x), type ="response"), col = "firebrick", lwd = 2.5, add = TRUE)
matricula <- as.character(datos$matricula) 
datos$matricula <- as.numeric(datos$matricula) 
# Se crea un vector con nuevos valores interpolados en el rango de observaciones. 
nuevos_puntos <- seq(from = min(datos$matematicas), to = max(datos$matematicas), by = 0.5) 
# Predicciones de los nuevos puntos según el modelo. 
# Si se indica se.fit = TRUE se devuelve el error estándar de cada predicción 

# junto con el valor de la predicción (fit). 
predicciones <- predict(modelo, data.frame(matematicas = nuevos_puntos), se.fit = TRUE) 
# Mediante la función logit se transforman los log_ODDs a probabilidades. 
predicciones_logit <- exp(predicciones$fit) / (1 + exp(predicciones$fit)) # Se calcula el límite inferior y superior del IC del 95% sustrayendo e # incrementando el logODDs de cada predicción 1.95*SE. Una vez calculados los 
# logODDs del intervalo se transforman en probabilidades con la función logit. 

limite_inferior <- predicciones$fit - 1.96 * predicciones$se.fit

limite_inferior_logit <- exp(limite_inferior) / (1 + exp(limite_inferior)) 

limite_superior <- predicciones$fit + 1.96 * predicciones$se.fit

limite_superior_logit <- exp(limite_superior) / (1 + exp(limite_superior)) 
# Se crea un data frame con los nuevos puntos y sus predicciones
datos_curva <- data.frame(matematicas = nuevos_puntos, probabilidad_matricula = predicciones_logit, limite_inferior_logit = limite_inferior_logit, limite_superior_logit =limite_superior_logit)
ggplot(datos, aes(x = matematicas, y = matricula)) + geom_point(aes(color = as.factor(matricula)), shape = "I", size = 3) + geom_line(data = datos_curva, aes(y = probabilidad_matricula), color = "firebrick") + geom_line(data = datos_curva, aes(y = limite_inferior_logit), linetype = "dashed") + geom_line(data = datos_curva, aes(y = limite_superior_logit), linetype = "dashed") + theme_bw() + labs(title = "Modelo regresión logística matrícula ~ nota matemáticas", y = "P(matrícula = 1 | matemáticas)", y = "matemáticas") + theme(legend.position = "null") + theme(plot.title = element_text(hjust = 0.5))

Evaluación del modelo

A la hora de evaluar la validez y calidad de un modelo de regresión logística, se analiza tanto el modelo en su conjunto como los predictores que lo forman.

# Diferencia de residuos # En R, un objeto glm almacena la "deviance" del modelo, así como la "deviance" # del modelo nulo.
dif_residuos <- modelo$null.deviance - modelo$deviance 
# Grados libertad 
df <- modelo$df.null - modelo$df.residual 
# p-value 4
p_value <- pchisq(q = dif_residuos,df = df, lower.tail = FALSE) 

paste("Diferencia de residuos:", round(dif_residuos, 4))
## [1] "Diferencia de residuos: 55.6368"
paste("Grados de libertad:", df)
## [1] "Grados de libertad: 1"
paste("p-value:", p_value)
## [1] "p-value: 8.71759108087094e-14"
# El mismo cálculo se puede obtener directamente con: 
anova(modelo, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: matricula
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          199     222.71              
## matematicas  1   55.637       198     167.07 8.718e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo comienza con una desviación total (Deviance) de 222.71 en 199 grados de libertad (Df) cuando se usa solo el valor NULL (ningún predictor). Luego, se agrega el predictor “matematicas” al modelo y se muestra que la desviación disminuye a 167.07, con un cambio en la desviación de 55.637. El número de grados de libertad residuales (Resid. Df) también disminuye a 198.

El valor de “Pr(> Chi)” muestra la probabilidad de que el cambio en la desviación sea debido solo al azar. En este caso, el valor de “Pr(> Chi)” es muy pequeño (8.718e-14) y se indica con tres asteriscos “***”. Esto sugiere que la adición de “matematicas” al modelo mejora significativamente el ajuste del modelo a los datos y que es poco probable que el cambio en la desviación sea debido solo al azar.

En resumen, la tabla de análisis de desviación muestra cómo cada término predictor afecta el modelo y puede ayudar a determinar qué variables son importantes para predecir la variable de respuesta. En este caso, la adición de “matematicas” al modelo parece ser una adición importante y significativa.

• McFadden’s: 𝑅𝑀𝑐𝐹2=1–𝑙𝑛𝐿^(𝑚𝑜𝑑𝑒𝑙𝑜)𝑙𝑛𝐿^(𝑚𝑜𝑑𝑒𝑙𝑜 𝑛𝑢𝑙𝑜), siendo 𝐿^ el valor de likelihood de cada modelo. La idea de esta fórmula es que, 𝑙𝑛(𝐿^), tiene un significado análogo a la suma de cuadrados de la regresión lineal. De ahí que se le denomine 𝑝𝑠𝑒𝑢𝑑𝑜𝑅2.

• Otra opción bastante extendida es el test de Hosmer-Lemeshow. Este test examina mediante un Pearson chi-square test si las proporciones de eventos observados son similares a las probabilidades predichas por el modelo, haciendo subgrupos.

Comparación de clasificación predicha y observaciones

library(vcd) 
## Warning: package 'vcd' was built under R version 4.2.2
## Loading required package: grid
predicciones <- ifelse(test = modelo$fitted.values > 0.5, yes = 1, no = 0) 
matriz_confusion <- table(modelo$model$matricula, predicciones, dnn = c("observaciones", "predicciones"))
matriz_confusion
##              predicciones
## observaciones   0   1
##             0 140  11
##             1  27  22
mosaic(matriz_confusion, shade = T, colorize = T, gp = gpar(fill = matrix(c("green3", "red2", "red2", "green3"), 2, 2)))

Esta tabla muestra una matriz de confusión para las predicciones del modelo. La matriz de confusión es una tabla que se utiliza para evaluar la precisión de un modelo de clasificación, en este caso, para predecir la variable binaria “matricula” (0 o 1).

Las filas de la tabla representan los valores observados de la variable “matricula”, mientras que las columnas representan las predicciones del modelo. En este caso, la tabla muestra que:

La diagonal principal de la tabla representa las predicciones correctas del modelo (los verdaderos positivos y verdaderos negativos), mientras que las celdas fuera de la diagonal principal representan las predicciones incorrectas (los falsos positivos y falsos negativos).

La matriz de confusión puede ser útil para calcular métricas de evaluación del modelo, como la precisión, la sensibilidad, la especificidad y el valor predictivo positivo y negativo.

#El modelo es xapaz de clasificar correctamente
p=(140+22)/(140+22+27+11)
p
## [1] 0.81

0.81(81%) de las observaciones cuando se emplean los datos de entrenamiento. No hay que olvidar que este es el error de entrenamiento, por lo que no es generalizable a nuevas observaciones. Para obtener una estimación más realista hay que calcular el error de test.

Conclusión