MeCoBi: Modelos de Regresion Logística
D. S. Fernandez del Viso
Septiembre 2019
Regresión Logística
En una regresión logística la variable respuesta (dependiente) es una variable binaria (dicótoma en términos generales), y las variables predictoras (independientes) pueden ser binarias, categóricas o continuas.
Logistic Regression vs Linear Regression 
Como la variable respuesta solo puede tener valores entre 0 o 1, y las variables independientes pueden ser continuas, con cualquier valor real, para formular un modelo de regresión logística se hace una transformación de probabilidades a ‘razones de probabilidades’ (en inglés odd ratios (ORs)), o proporción de casos favorables a desfavorables.
La transformación tiene esta forma:
\[OR_i = \frac{\pi}{1-\pi} \qquad (1)\] donde \(\pi\) es la probabilidad de que el evento dependiente (\(Y\)) ocurra; luego se calcula el logaritmo del \(OR\), el cual se denomina logit:
\[logit(\pi) = log\frac{\pi}{1-\pi} \qquad (2)\]
El modelo para la regresión, con dos variables independientes (\(x_1\), \(x_2\)), es:
\[logit(\pi)=\beta_0 + x_1\beta_1 +x_2\beta_2\qquad (3)\]
EJEMPLO 1: Predicción de alta presión en mujeres
Vamos a explorar un modelo de diagnóstico de alta presión, a partir de datos de mujeres (n = 189), en relación a su edad, estado de menopausia, y el índice de masa corporal. Los datos están codificados como variables binarias:
| Variables | Nombre | Código |
|---|---|---|
| Edad (años) | age | 0:≤50, 1>50 |
| Alta presión | dxhigh | 0:no, 1:sí |
| BMI (kg/m2) | bmi | 0:≤25, 1:>25 |
| Menopausia | menop | 0:pre-, 1:menopausia |
Para este modelo utilizaremos la función glm (general linear model):
#leemos los datos:
library(readxl)
hbp <- read_excel("mod_empiricos.xlsx",
sheet = "hbp")Modelo logístico con glm
#modelo logístico:
highbp <- glm(dxhigh ~ age + bmi + menop, data = hbp,
family=binomial(logit))
#resumen de los resultados:
summary(highbp)##
## Call:
## glm(formula = dxhigh ~ age + bmi + menop, family = binomial(logit),
## data = hbp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0392 -0.8084 -0.6580 1.3222 2.0781
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0367 0.4550 -4.476 7.61e-06 ***
## age 0.6376 0.4920 1.296 0.195
## bmi 0.6165 0.4620 1.334 0.182
## menop 0.4484 0.4887 0.918 0.359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 218.39 on 188 degrees of freedom
## Residual deviance: 206.60 on 185 degrees of freedom
## AIC: 214.6
##
## Number of Fisher Scoring iterations: 4#obtener los coeficientes como ORs:
exp(cbind(OR = coef(highbp), confint(highbp)))## OR 2.5 % 97.5 %
## (Intercept) 0.1304608 0.0491431 0.2984743
## age 1.8919266 0.7125247 4.9772745
## bmi 1.8523471 0.7836975 4.9120329
## menop 1.5658815 0.6016810 4.1448987EJEMPLO 2: Variables predictoras de diabetes tipo 2, en una población de nativos Pima
En este caso las variables independientes no son binarias (0,1) sino datos continuos. No se calculan los \(OR\) (ni logit), y el resultado son sólo los coeficientes de la regresión.
Datos
# paquetes
library(tidyverse)
library(caret)
library(mlbench)
# cargar datos (estos son directamente del paquete mlbench, en su caso debe usar readxl o similar)
data("PimaIndiansDiabetes2", package = "mlbench")
# inspeccionar los datos
sample_n(PimaIndiansDiabetes2, 6)## pregnant glucose pressure triceps insulin mass pedigree age diabetes
## 1 1 106 76 NA NA 37.5 0.197 26 neg
## 2 1 116 78 29 180 36.1 0.496 25 neg
## 3 1 88 62 24 44 29.9 0.422 23 neg
## 4 4 184 78 39 277 37.0 0.264 31 pos
## 5 1 143 86 30 330 30.1 0.892 23 neg
## 6 3 173 82 48 465 38.4 2.137 25 posDividiendo los datos en “training” y “test”
# Split the data into training and test set
set.seed(123)
training.samples <- PimaIndiansDiabetes2$diabetes %>%
createDataPartition(p = 0.8, list = FALSE)
train.data <- PimaIndiansDiabetes2[training.samples, ]
test.data <- PimaIndiansDiabetes2[-training.samples, ]Modelo usando glm
library(MASS)
# Fit the model
model <- glm(diabetes ~., data = train.data, family = binomial)
# Summarize the final selected model
summary(model)##
## Call:
## glm(formula = diabetes ~ ., family = binomial, data = train.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6751 -0.6666 -0.3588 0.6581 2.5650
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.196902 1.369247 -7.447 9.54e-14 ***
## pregnant 0.082032 0.061219 1.340 0.18025
## glucose 0.036517 0.006381 5.723 1.05e-08 ***
## pressure 0.001333 0.013053 0.102 0.91863
## triceps 0.008425 0.018649 0.452 0.65145
## insulin -0.001219 0.001441 -0.845 0.39784
## mass 0.081434 0.030448 2.675 0.00748 **
## pedigree 1.186528 0.470790 2.520 0.01173 *
## age 0.030886 0.020337 1.519 0.12884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 402.38 on 314 degrees of freedom
## Residual deviance: 280.83 on 306 degrees of freedom
## (300 observations deleted due to missingness)
## AIC: 298.83
##
## Number of Fisher Scoring iterations: 5Gráfica logística diabetes vs glucosa
library(ggplot2)
#pasar datos de diabetes "pos" y "neg" a 1s y 0s
diabetes01 <- ifelse(PimaIndiansDiabetes2$diabetes == "pos", 1, 0)
#gráfica con curva logística
ggplot(PimaIndiansDiabetes2, aes(x=glucose, y=diabetes01, na.rm = TRUE)) +
geom_point() +
geom_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE)
Predicciones usando el 20% de los datos
# Make predictions
probabilities <- model %>% predict(test.data, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, "pos", "neg")
# Prediction accuracy
observed.classes <- test.data$diabetes
mean(predicted.classes == observed.classes, na.rm = TRUE)## [1] 0.8181818