GLM Regresión Logística en R
Yolanda Larriba
2025-03-01
En esta práctica aprenderemos a ajustar el modelo GLM de Regresión Logística en R a través de distintos ejemplos que incluyen tanto datos reales como simulados. Se particulariza en la correcta interpretación de los resultados del modelo ajustado, así como en la validación y comparación de modelos, inferencias sobre los parámetros y detección de puntos influyentes y/o atípicos. El objetivo es que el estudiante adquiera un dominio básico para abordar problemas similares en la práctica.
1. Ajuste del modelo, estimación e interpretación. Datos simulados.
Se desea predecir la probabilidad de conseguir un empleo en función del nivel de educación y si una persona pertenece a un grupo sujeto a discriminación, a partir del siguiente modelo:
\[ y_i = 1.5 \cdot \text{education}_i - 2 d_i + \epsilon_i \]
donde:
- \(\text{education}_i \sim N(0,1)\) es el nivel educativo del individuo, que asumimos sigue una distribución normal estándar.
- \(d_i\) es una variable binaria que indica si el individuo pertenece (1) o no (0) a un grupo sujeto a discriminación. Se asume que la probabilidad de estar en el grupo es del 50%.
- \(\epsilon_i \sim N(0,1)\) es un término de error aleatorio que sigue una distribución normal estándar.
- Se simulan 500 casos para realizar el análisis.
Ajustar un modelo de regresión logística para estimar la probabilidad de conseguir empleo y analizar la influencia de la educación y la discriminación en este proceso.
Generación de Datos
Simularemos un conjunto de datos en el que:
education: Nivel educativo (variable continua, distribuida normal estándar).d: Grupo de pertenencia (0 = no discriminado, 1 = discriminado).y: Variable auxiliar en la escala logit (log-odds).y_proba: Probabilidad transformada a la escala de probabilidades mediante la función logística.y_dummy: Variable binaria final (0 = no consigue empleo, 1 = consigue empleo).
# Cargar paquetes necesarios
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggeffects)## Warning: package 'ggeffects' was built under R version 4.3.3# Fijar semilla para reproducibilidad
#set.seed(123)
# Número de observaciones
N <- 500
# Simulación de predictores
education <- rnorm(N, 0, 1) # Educación como variable continua
d <- rbinom(N, 1, 0.5) # Grupo discriminado (1) o no (0)
# Modelo lineal en escala logit
error <- rnorm(N) # Término de error aleatorio
y <- 1.5 * education + (-2 * d) + error # Modelo lineal en log-odds
# Transformación logística para obtener probabilidades
y_proba <- plogis(y)
# Generación de la variable binaria a partir de las probabilidades
y_dummy <- rbinom(n = N, size = 1, prob = y_proba)
# Construcción del dataframe final
df_sim <- data.frame(id = 1:N, education, d, y, y_proba, y_dummy)
# Inspección de los primeros registros
head(df_sim)## id education d y y_proba y_dummy
## 1 1 0.49702691 1 -1.4123985 0.195856025 1
## 2 2 -0.52738012 0 -1.5234206 0.178958366 1
## 3 3 0.69159700 1 -0.2728023 0.432219276 1
## 4 4 0.02129867 0 0.8791121 0.706638201 0
## 5 5 -1.76493127 1 -5.5133306 0.004016457 0
## 6 6 -0.70236515 1 -0.3511817 0.413095900 1plot(education,y_dummy)
Ajuste del modelo de regresión logística
# Ajustar el modelo GLM de regresión logística
glm_model <- glm(y_dummy ~ education + d, data = df_sim, family = binomial)
# Resumen del modelo
summary(glm_model)##
## Call:
## glm(formula = y_dummy ~ education + d, family = binomial, data = df_sim)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.005065 0.153121 0.033 0.974
## education 1.307326 0.142802 9.155 < 2e-16 ***
## d -1.649321 0.236216 -6.982 2.91e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 644.94 on 499 degrees of freedom
## Residual deviance: 488.46 on 497 degrees of freedom
## AIC: 494.46
##
## Number of Fisher Scoring iterations: 5# Transformar coeficientes a razones de odds
exp(coef(glm_model)) # Interpretación en términos de odds## (Intercept) education d
## 1.0050781 3.6962773 0.1921804# Education:Un aumento de 1 unidad en educación multiplica las probabilidades de conseguir empleo por exp(β1).
#Discriminación: Pertenecer al grupo discriminado reduce las probabilidades de empleo en un factor de exp(β2), dado que el coeficiente es negativo.Comparación con el modelo lineal
# Ajustar modelo de regresión lineal
lm_model <- lm(y_dummy ~ education + d, data = df_sim)
# Resumen del modelo
summary(lm_model)##
## Call:
## lm(formula = y_dummy ~ education + d, data = df_sim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90094 -0.32392 -0.09848 0.34877 1.02121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.49950 0.02716 18.391 < 2e-16 ***
## education 0.21591 0.01883 11.464 < 2e-16 ***
## d -0.28073 0.03663 -7.663 9.56e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4073 on 497 degrees of freedom
## Multiple R-squared: 0.2712, Adjusted R-squared: 0.2683
## F-statistic: 92.48 on 2 and 497 DF, p-value: < 2.2e-16# La regresión lineal puede producir probabilidades fuera del rango [0,1], lo cual es problemático.Predicción y visualización de resultados
# Predicción de probabilidades
df_sim$pred_prob <- predict(glm_model, type = "response")
# Comparación de probabilidades reales vs. predichas
head(df_sim[c("y_proba", "pred_prob")])## y_proba pred_prob
## 1 0.195856025 0.27002828
## 2 0.178958366 0.33528163
## 3 0.432219276 0.32297954
## 4 0.706638201 0.50822665
## 5 0.004016457 0.01886140
## 6 0.413095900 0.07159281# Gráfico de predicciones con ggpredict()
g_pred <- ggpredict(glm_model, terms = c("education[all]", "d"))
g_pred <- as.data.frame(g_pred)
# Visualización con ggplot
ggplot(g_pred, aes(x, predicted, colour = group)) +
geom_point() + geom_line() +
labs(x = "Nivel de Educación", y = "Probabilidad de Conseguir Empleo") +
theme_minimal()
2. Interpretación de los coeficientes según tipo de variable. Datos birth
Se tiene información sobre el peso de 189 recién nacidos y 10 variables clínicas de sus madres. El objetivo es describir la ocurrencia de bajo peso al nacer (\(Y\), variable respuesta dicotómica) atendiendo la raza, ser o no madre fumadora, el peso de la madre en el último ciclo menstrual, o la edad entre otros factores.
Carga de datos
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectdata("birthwt") # Cargamos los datos
dim(birthwt)## [1] 189 10head(birthwt)## low age lwt race smoke ptl ht ui ftv bwt
## 85 0 19 182 2 0 0 0 1 0 2523
## 86 0 33 155 3 0 0 0 0 3 2551
## 87 0 20 105 1 1 0 0 0 1 2557
## 88 0 21 108 1 1 0 0 1 2 2594
## 89 0 18 107 1 1 0 0 1 0 2600
## 91 0 21 124 3 0 0 0 0 0 2622attach(birthwt) # Permite trabajar con los nombres de las variables del data.frame directamente
#-----------------------------------------------------------------
# low Peso al nacer (0 = Peso $>$ 2500g,
# 1 = Peso $<$ 2500g)
# age Edad de la madre en agnos
# lwt Peso de la madre en el ultimo ciclo menstrual
# race Raza (1 = Blanca, 2 = Negra, 3 = Otra)
# smoke Fumadora durante el embarazo (1 = Si, 0 = No)
#-----------------------------------------------------------------Variable independiente dicotómica
Comenzaremos con este caso ya que es la base para los demás tipos de variables independientes. Sea \(X\) la variable independiente codificada como 0 o 1 en nuestro modelo de regresión logística, por ejemplo la variable madre fumadora: \[\ln(\frac{p}{1-p})=\beta_0+\beta_1X.\]
La diferencia en el logit para un individuo con \(X=0\) y \(X=1\) es \(\beta_1\).
Si en la ecuación anterior despejamos \(p,\) se obtiene: \[p=\frac{e^{\beta_0+\beta_1X}}{1+e^{\beta_0+\beta_1X}} \text{, y por ende, } 1-p=\frac{1}{1+e^{\beta_0+\beta_1X}}\]
Los distintos valores de las probabilidades para las 4 combinaciones entre la variable dependiente e independiente son:
De acuerdo con esta tabla, y volviendo a nuestro ejemplo, la odds, o ventaja, de que la variable respuesta \(Y\) (bajo peso al nacer) esté presente en los individuos con \(X=1\) (o madre fumadora) se define como: \[\frac{Pr(Y=1|X=1)}{Pr(Y=0|X=1)}=\frac{Pr(Y=1|X=1)}{1-Pr(Y=1|X=1)}\]
De manera análoga, la odds, o ventaja, de que la variable respuesta \(Y\) (bajo peso al nacer) esté presente en los individuos con \(X=0\) (o madre no fumadora) se define como: \[\frac{Pr(Y=1|X=0)}{Pr(Y=0|X=0)}=\frac{Pr(Y=1|X=0)}{1-Pr(Y=1|X=0)}\]
A partir de estas odds, podemos definir la odds ratio, o razón de ventajas, (OR) que es una medida de asociación entre dos variables cualitativas. En términos formales, se define como la probabilidad que una condición de salud o enfermedad (bajo peso al nacer) se presente en un grupo de población expuesto a un agente o factor de riesgo (madre fumadora), sobre la probabilidad de que ocurra en otro grupo sin exposición a dicho riesgo. La OR se define como:\[OR=\frac{Pr(Y=1|X=1)/1-Pr(Y=1|X=1)}{Pr(Y=1|X=0)/1-Pr(Y=1|X=0)}\]
Una OR = 1 implica que no existe asociación entre las variables consideradas, es decir, que la exposición al factor de riesgo no pone a riesgo de enfermar ni protege contra la enfermedad. Los valores mayores de 1 indican que la exposición pone a riesgo de enfermar, mientras que los valores inferiores a 1 indican que ésta protege contra la enfermedad.
Utilizando la tabla anterior de probabilidades se tiene: \[OR=\frac{(\frac{e^{\beta_0+\beta1}}{1+e^{\beta_0+\beta_1}})/(\frac{1}{1+e^{\beta_0+\beta_1}})}{(\frac{e^{\beta_0}}{1+e^{\beta_0}})/(\frac{1}{1+e^{\beta_0}})}=\frac{e^{\beta_0+\beta_1}}{e^{\beta_0}}=e^{\beta_1}\]
Nótese que la OR coincide con la exponencial del del coficiente \(\beta_1\) del modelo de regresión logística, y cuadra con la interpretación de los coeficientes que haremos en términos de odds o ventajas. Y recordemos por otro lado, que tal y como hemos visto, \(\beta_1\) es la diferencia en el logit para un individuo con \(X=0\) y \(X=1.\)
Por ejemplo, supongamos que se quiere estudiar la relación entre madre fumadora y bajo peso al nacer:
# Variables categóricas consideradas como tal en R deben definirse como tipo factor con f niveles (categorias)
# Se incluyen f-1 variables dummy en la matriz de diseño
# La categoría de referencia es el nivel/categoria primera/mas baja
smoke<-factor(smoke)
race<-factor(race)
# Procedimiento para ajustar el modelo en R
logistica1=glm(low~smoke+race,family=binomial(link=logit))
summary(logistica1)##
## Call:
## glm(formula = low ~ smoke + race, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8405 0.3529 -5.216 1.83e-07 ***
## smoke1 1.1160 0.3692 3.023 0.00251 **
## race2 1.0841 0.4900 2.212 0.02693 *
## race3 1.1086 0.4003 2.769 0.00562 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 219.97 on 185 degrees of freedom
## AIC: 227.97
##
## Number of Fisher Scoring iterations: 4table(low,smoke)## smoke
## low 0 1
## 0 86 44
## 1 29 30# Pr(Y=1|X=1)=30/74; Pr(Y=0|X=1)=44/74
# Pr(Y=1|X=0)=29/115; Pr(Y=0|X=0)=86/115
OR<-((30/74)/(44/74))/((29/115)/(86/115)) # OR=(86*30)/(29*44)
OR## [1] 2.021944logistica2<-glm(low~smoke,family=binomial(link=logit))
summary(logistica2)##
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
## smoke1 0.7041 0.3196 2.203 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.80 on 187 degrees of freedom
## AIC: 233.8
##
## Number of Fisher Scoring iterations: 4exp(logistica2$coefficients[2])## smoke1
## 2.021944exp(coefficients(logistica2))## (Intercept) smoke1
## 0.3372093 2.0219436Una OR=2.02 se interpreta como que la ocurrencia de bajo peso en recién nacidos es 2.02 veces mayor en los recién nacidos con madre fumdora.
Junto con el estimador de la OR , es importante utilizar los intervalos de confianza (IC) para obtener información adicional acerca del parámetro. Los IC para los coeficientes del modelo se calculan asumiendo que la distribución de los parámetros estimados es aproximadamente normal. Entonces un IC para la OR se calcula tomando exponenciales en el IC para los parámetros:\[exp\lbrace \hat{\beta_1}\pm z_{\alpha/2}\hat{se}(\hat{\beta_1})\rbrace\]
exp(coefficients(logistica2))## (Intercept) smoke1
## 0.3372093 2.0219436exp(confint(logistica2)) # IC## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.2177709 0.5070199
## smoke1 1.0818724 3.8005817Que el IC para el OR sea (1.08,3.80) quiere decir que con una confianza del 95% la ocurrencia de bajo peso al nacer en los recién nacidos de madre fumadora en la población de estudio es entre 1.08 y 3.80 veces más problable que en los recién nacidos con madre no fumadora.
Variable independiente politómica
En este caso la variable independiente tiene más de 2 categorías. En nuestro ejemplo, lo ilustraremos con la variable race.
logistica3=glm(low~race,family=binomial(link=logit))
summary(logistica3)##
## Call:
## glm(formula = low ~ race, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1550 0.2391 -4.830 1.36e-06 ***
## race2 0.8448 0.4634 1.823 0.0683 .
## race3 0.6362 0.3478 1.829 0.0674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.66 on 186 degrees of freedom
## AIC: 235.66
##
## Number of Fisher Scoring iterations: 4exp(coefficients(logistica3))## (Intercept) race2 race3
## 0.3150685 2.3275362 1.8892340Para el ajuste de este modelo, R crea dos variables dicotómicas (valores 0/1) que indican si la persona es o no de una determinada raza, a estas variables se las conoce como variables dummy y se incorporan a la matriz de diseño del GLM. Por defecto, R asume que la primera categoría de la variable (categoría de referencia) toma siempre el valor 0, y compara el resto de grupos con este asignado así los distintos valores 0/1 de las variables dummy. En nuestro caso, la categoría de referencia es race=1 (blanca). La OR de recién nacido de bajo peso para race2 (negra) es la odds de las mujeres de race=2 dividida por la odds de las mujeres de race=1 (blanca), 2.33, esto significa que estimamos que la \(OR\) de bajo peso al nacer para mujeres de race=2 es 2.33 veces el de las mujeres con race=1. Los I.C. se obtienen siguiendo el mismo procedimiento que en el caso de las variables dicotómicas.
exp(confint(logistica3))## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.1929856 0.4950281
## race2 0.9255511 5.7746995
## race3 0.9565420 3.7578709Variable independiente continua
Cuando en un modelo de regresión logística utilizamos una variable predictora continua, la interpretación del parámetro estimado depende de las unidades en que se mide esta variable. En estos casos, el parámetro \(\beta_1\) representa el cambio en el log-odds o logit cuando la variable independiente (X) cambia en una unidad, o equivalentemente, \(e^{\beta_1}\) es el cambio en la odds cuando la variable X cambia en una unidad. Si en vez de cambiar en una unidad, cambia en \(c\) unidades, entonces el cambio vendrá dado por \(e^{c\beta_1}.\)
Por ejemplo, si se desea a establecer la relación entre el peso de la madre (en libras) en la última menstruación (lwt) y el bajo peso al nacer:
logistica4=glm(low~lwt,family=binomial(link=logit))# lwt como factor?
summary(logistica4)##
## Call:
## glm(formula = low ~ lwt, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99831 0.78529 1.271 0.2036
## lwt -0.01406 0.00617 -2.279 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 228.69 on 187 degrees of freedom
## AIC: 232.69
##
## Number of Fisher Scoring iterations: 4exp(coefficients(logistica4))## (Intercept) lwt
## 2.7137035 0.9860401Por tanto, por cada por cada unidad (en libras) que aumenta el peso de la madre en su última menstruación, el cambio en la odds de ocurrencia de recién nacido con bajo peso es 0.986 veces menor (efecto protector). Si aumenta 10 unidades, la odds de recién nacido con bajo peso sería \(0.986^{10}=0.868\) veces menor. Podemos ver los resultados en un gráfico:
ajustados4=predict(logistica4,type="response")# type=response, devuelve las probabilidades; type=link logit
plot(lwt,ajustados4)
La curva es decreciente ya que la odds < 1. Hay que tener cuidado al interpretar los resultados para lwt altos ya que están basados en muy pocas observaciones.
Variables independientes categóricas y continuas
Hasta ahora hemos visto como ajustar modelos de regresión logística con una sola variable independiente, a estos modelos se les llama univariantes. Este tipo de modelos son válidos en pocas ocasiones, ya que en general, una variable independiente está asociada a otras variables independientes. Por eso, es necesario considerar un modelo multivariante para poder compreneder mejor lo que ocurre con los datos. El objetivo es ajustar el efecto de cada variable en el modelo por las otras variables independientes. Por lo tanto, en el modelo de regresión logística, cada coeficiente del modelo proporciona una estimación de la log-odds (logit) ajustando por las otras variables incluidas en el modelo.
Por ejemplo, consideremos un modelo para explicar la relación entre bajo peso al nacer con madre fumadora (\(X_1\)) y el peso de ésta en su última menstruación (\(X_2\)): \[\ln(\frac{p}{1-p})=\beta_0+\beta_1X_1+\beta_2X_2\]
logistica5=glm(low~smoke,family=binomial(link=logit))
logistica6=glm(low~smoke+lwt,family=binomial(link=logit))
summary(logistica6)##
## Call:
## glm(formula = low ~ smoke + lwt, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.62200 0.79592 0.781 0.4345
## smoke1 0.67667 0.32470 2.084 0.0372 *
## lwt -0.01332 0.00609 -2.188 0.0287 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 224.34 on 186 degrees of freedom
## AIC: 230.34
##
## Number of Fisher Scoring iterations: 4En este caso, \(e^{0.677}=1.968\) es la OR para un mismo valor de peso (lwt) de la madre en su última menstruación, es decir, la OR que esperaríamos encontrar si el peso medio en madres fumadoras y no fumadoras fuese similar. Nótese que en estos casos, la interpretación de la \(OR\) siempre debe de tener en cuenta la la presencia de otras variables incluidas en el modelo.
Concretamente, se observa que el parámetro \(\beta_1\) es muy similar en el modelo que solo considera smoke y en el modelo con smoke y lwt (logistica5 y logistica6), esto es debido a que el peso medio de madres fumadoras y no fumadoras es muy similar.
summary(logistica5)##
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
## smoke1 0.7041 0.3196 2.203 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.80 on 187 degrees of freedom
## AIC: 233.8
##
## Number of Fisher Scoring iterations: 4mean(birthwt[birthwt[,"smoke"]==0,"lwt"])#kg=lb/2.2046## [1] 130.8957mean(birthwt[birthwt[,"smoke"]==1,"lwt"])## [1] 128.1351Por tanto ajustar por el peso de la madre en la última regla no cambia el log-odds de bajo peso en recién nacido de madre fumadora/no fumadora (se observan líneas paralelas entre los valores de la variable respuesta y la covariable para cada nivel del factor, no aumenta ni disminuye la tendencia). Tal y como concretaremos posteriormente, diremos que en este modelo no hay interacción.
ajustados6=predict(logistica6,type="link")
plot(lwt,ajustados6,type="n",main="sin interaccion")
points(lwt[smoke==0],ajustados6[smoke==0],col=2)
points(lwt[smoke==1],ajustados6[smoke==1],col=4)
legend(180,0, col=c(2,4),pch=c(1,1),c("no fumadora","fumadora"))
Confusión e interacción
El término confusión se utiliza para describir una covariable que está asociada al mismo tiempo a la variable respuesta y al factor de interés. Cuando estas dos asociaciones están presentes la relación entre la variable respuesta y el factor están confundidos. Una manera de detectar una variable de confusión, siempre que en el modelo no haya interacción, es ver si los parámetros estimados para el factor de interés cambian sustancialmente al introducir la covariable.
En el caso de bajo peso de recién nacido, para saber si el peso de la madre en la última regla (lwt) es un variable de confusión con respecto a madre fumadora, vemos si el valor del parámetro asociado a fumar (smoke) cambia mucho al introducir lwt (puesto que hemos visto que en este modelo no hay interacción):
summary(logistica5)##
## Call:
## glm(formula = low ~ smoke, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
## smoke1 0.7041 0.3196 2.203 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.80 on 187 degrees of freedom
## AIC: 233.8
##
## Number of Fisher Scoring iterations: 4summary(logistica6)##
## Call:
## glm(formula = low ~ smoke + lwt, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.62200 0.79592 0.781 0.4345
## smoke1 0.67667 0.32470 2.084 0.0372 *
## lwt -0.01332 0.00609 -2.188 0.0287 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 224.34 on 186 degrees of freedom
## AIC: 230.34
##
## Number of Fisher Scoring iterations: 4En este caso, hay muy poca diferencia, por lo que podemos considerar que no hay confusión. En concreto, la estimación del parámetro pasa de 0.704 a 0.677 por lo que no parece que lwt sea una variable de confusión (un decremento/incremento de más de 25% sería un indicio de variable de confusión).
Por otro lado, cuando la asociación entre una covariable y la variable respuesta es la misma para cada nivel del factor, se dice que no hay interacción entre la covariable y el factor. Gráficamente, la ausencia de interacción da lugar a un modelo con dos rectas paralelas, una para cada nivel del factor. Cuando la interacción está presente, la asociación entre el factor y la variable respuesta depende del nivel de la covariable, es decir la covariable modifica el efecto del factor. En este caso las rectas no son paralelas, y se dice que la covariable tiene un efecto modificador. Siendo imposible estimar la OR para el factor sin especificar el valor de la covariable para el que se hará dicha comparación.
En nuestro ejemplo, ¿existe interacción entre el factor madre fumadora y el peso de éstas en su última regla al estudiar la asociación de ambas con la ocurrencia de recién nacido de bajo peso?
logistica7=glm(low~smoke+lwt+smoke:lwt,family=binomial(link=logit))
ajustados7=predict(logistica7,type="link")
summary(logistica7)##
## Call:
## glm(formula = low ~ smoke + lwt + smoke:lwt, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.93234 1.29209 1.496 0.1348
## smoke1 -1.51089 1.61737 -0.934 0.3502
## lwt -0.02389 0.01039 -2.299 0.0215 *
## smoke1:lwt 0.01757 0.01279 1.373 0.1697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 222.37 on 185 degrees of freedom
## AIC: 230.37
##
## Number of Fisher Scoring iterations: 4plot(lwt,ajustados7,type="n",main="con interaccion")
points(lwt[smoke==0],ajustados7[smoke==0],col=4,pch=1)
points(lwt[smoke==1],ajustados7[smoke==1],col=2,pch=1)
legend(180,0, col=c(2,4,3),pch=c(1,1,2),c("no fumadora","fumadora"))
Ya vimos que la asociación entre smoke y el recién nacido de bajo peso era la misma para cada nivel del factor. De hecho, al incluir un término de interacción en el modelo, el parámetro asociado a smoke deja de ser significativo, y lo más importante, la interacción no es significativa. En el gráfico vemos que las rectas a las que da lugar el modelo con interacción no son paralelas, pero más importante el coeficiente asociado a la interacción no es significativo y no debe considerarse.
En la práctica, cualquier cambio significativo en el coeficiente del factor de riesgo, sugiere que la variable es de confusión. Si esto ocurre y la interacción no es estadísticamente significativa, la variable ha de salir del modelo. Por otra parte, una variable es un efecto modificador, solo cuando el término de interacción añadido al modelo es estadísticamente significativo. En cualquier caso, ambas conclusiones hay que comprobarlas mediante tests de hipótesis que comparen los modelos.
Interpretación de la OR en presencia de interacción
Cuando existe interacción entre un factor de riesgo y otra variable, el parámetro estimado para el factor de riesgo depende del valor de la variable que interactúa con éste. Luego, no es posible obtener la OR simplemente exponenciando el valor del parámetro, los pasos a seguir son:
- Escribir la ecuación del logit para los dos niveles del factor de riesgo.
- Calcular la diferencia entre los logits.
- Exponenciar el valor obtenido.
Por ejemplo, en un modelo que contiene solo dos variables y su interacción. Llamamos \(F\) al factor, \(X\) a la covariable y \(F\times X\) a la interacción. El logit para \(F=f\) y \(X=x\) es: \[\ln(\frac{p(f,x)}{1-p(f,x)})=\beta_0+\beta_1f+\beta_2x+\beta_3f\times x,\] si queremos la OR para los niveles de F \(f_1\) versus \(f_0\): \[ \begin{align} \ln(\frac{p(f_1,x)}{1-p(f_1,x)})&=\beta_0+\beta_1f_1+\beta_2x+\beta_3f_1\times x\\ \ln(\frac{p(f_0,x)}{1-p(f_0,x)})&=\beta_0+\beta_1f_0+\beta_2x+\beta_3f_0\times x\\ \ln(\frac{p(f_1,x)}{1-p(f_1,x)})-\ln(\frac{p(f_0,x)}{1-p(f_0,x)})&=\beta_1(f_1-f_0)+\beta_3 x(f_1-f_0)\\ OR&=exp\lbrace \beta_1(f_1-f_0)+\beta_3 x(f_1-f_0) \rbrace \end{align} \] En nuestro modelo con interacción anterior (que asumiremos puntualmente significativa) entre el peso de la madre en su última menstruación y ser madre fumadora o no, la ventaja de recién nacido de bajo peso de madre fumadora frente a no fumadora para una madre de peso en la última menstruación \(x\) se calcula como:
summary(logistica7) # Asumiremos que sola interaccion es significativa para ilustrarlo##
## Call:
## glm(formula = low ~ smoke + lwt + smoke:lwt, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.93234 1.29209 1.496 0.1348
## smoke1 -1.51089 1.61737 -0.934 0.3502
## lwt -0.02389 0.01039 -2.299 0.0215 *
## smoke1:lwt 0.01757 0.01279 1.373 0.1697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 222.37 on 185 degrees of freedom
## AIC: 230.37
##
## Number of Fisher Scoring iterations: 4# OR se calcula:
#exp(-1.51089+0.01757*x)#x es un valor del peso de la madre3.2.7.- Interpretación de los valores ajustados
En general, en regresión logística, los parámetros y las OR son el centro de interés, aunque en ocasiones los valores ajustados pueden ser también importantes. Calcular los IC para el logit es sencillo, ya que R da los errores estándar. Para calcular el IC para las probabilidades sólo hay que utilizar la relación entre el logit y la probabilidad: \[\frac{e^{IC logit}}{1+e^{IC logit}}\] Por ejemplo, supongamos que nos interesa estudiar el efecto del peso de la mujeres fumadoras en su última menstruación sobre la ocurrencia de recién nacidos de bajo peso:
logistica6=glm(low~smoke+lwt,family=binomial(link=logit))
ajustados6=predict(logistica6, se.fit=TRUE)
names(ajustados6)## [1] "fit" "se.fit" "residual.scale"L.inf=with(ajustados6,exp(fit-1.96*se.fit)/(1+exp(fit-1.96*se.fit))) # Evalua un expresión de R con la función pasada como argumento
L.sup=with(ajustados6,exp(fit+1.96*se.fit)/(1+exp(fit+1.96*se.fit)))
with(ajustados6,plot(lwt[smoke==1],(exp(fit)/(1+exp(fit)))[smoke==1],
ylim=c(0,1),ylab="probabilidad"))
points(lwt[smoke==1],L.inf[smoke==1],col=2)
points(lwt[smoke==1],L.sup[smoke==1],col=2)
legend(150,0.9, c("probabilidad","I.C."),col=c(1,2),pch=c(1,1))
abline(v=mean(lwt[smoke==1]),lty=2)
abline(v=135)
# ¿Y en las no fumadoras?El IC es más ancho en los extremos ya que hay menos observaciones y por lo tanto las estimaciones son menos precisas. La línea vertical punteda corresponde al peso medio de las mujeres fumadoras, donde el IC es más estrecho. Por ejemplo, para un peso de 130 libras (línea vertical sin puntear), el valor estimado es 0.45 y el IC [0.3,0.55].
3. Selección de modelo y predicción. Datos Cleveland
Se quiere estudiar la probabilidad de padecer endermedad cardiaca a partir de una colección de variables sociodemográficas y biomédicas relacionadas con el estado de salud. ### Lectura de datos
getwd()## [1] "C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM"datos <- read.table(file="C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM/Datos/Cleveland.txt",sep=",")
colnames(datos) <- c('age', 'sex', 'cp', 'trestbps', 'chol',
'fbs', 'restecg', 'thalach', 'exang',
'oldpeak', 'slope', 'ca', 'thal', 'target')
# Exploremos los datos
str(datos)## 'data.frame': 303 obs. of 14 variables:
## $ age : int 63 67 67 37 41 56 62 57 63 53 ...
## $ sex : int 1 1 1 1 0 1 0 0 1 1 ...
## $ cp : int 1 4 4 3 2 2 4 4 4 4 ...
## $ trestbps: int 145 160 120 130 130 120 140 120 130 140 ...
## $ chol : int 233 286 229 250 204 236 268 354 254 203 ...
## $ fbs : int 1 0 0 0 0 0 0 0 0 1 ...
## $ restecg : int 2 2 2 0 2 0 2 0 2 2 ...
## $ thalach : int 150 108 129 187 172 178 160 163 147 155 ...
## $ exang : int 0 1 1 0 0 0 0 1 0 1 ...
## $ oldpeak : num 2.3 1.5 2.6 3.5 1.4 0.8 3.6 0.6 1.4 3.1 ...
## $ slope : int 3 2 2 3 1 1 3 1 2 3 ...
## $ ca : int 0 3 2 0 0 0 2 0 1 0 ...
## $ thal : int 6 3 7 3 3 3 3 3 7 7 ...
## $ target : int 0 2 1 0 0 0 3 0 2 1 ...# Se deben tener en cuenta algunas modificaciones
# Convertimos a factor las variables cualitativas
library(dplyr)
datos <- datos %>% mutate(sex = as.factor(sex),
cp = as.factor(cp),
fbs = as.factor(fbs),
restecg = as.factor(restecg),
exang=as.factor(exang),
slope=as.factor(slope),
thal=as.factor(thal))
# Crear la variable Y a partir de la variable target.
# Y=1=presence y Y=0=absence
datos <- datos %>% mutate(y=case_when(target == 0 ~ "absence",
target != 0 ~ "presence"))
# De nuevo hay que convertirlo a factor
datos <- datos %>% mutate(y = as.factor(y))
# Para ver una parte de los datos
datos %>% head()## age sex cp trestbps chol fbs restecg thalach exang oldpeak slope ca thal
## 1 63 1 1 145 233 1 2 150 0 2.3 3 0 6
## 2 67 1 4 160 286 0 2 108 1 1.5 2 3 3
## 3 67 1 4 120 229 0 2 129 1 2.6 2 2 7
## 4 37 1 3 130 250 0 0 187 0 3.5 3 0 3
## 5 41 0 2 130 204 0 2 172 0 1.4 1 0 3
## 6 56 1 2 120 236 0 0 178 0 0.8 1 0 3
## target y
## 1 0 absence
## 2 2 presence
## 3 1 presence
## 4 0 absence
## 5 0 absence
## 6 0 absenceSelección del modelo con backward y forward
# Seleccion FORWARD -------------------------------------------------------
empty_model <- glm(y ~ 1, data=datos, family=binomial(link='logit'))
# Vamos a crear una formula para el modelo full
horizonte <- formula(y ~ age + sex + cp + trestbps + chol + fbs +
restecg + thalach + exang + oldpeak + slope)
library(MASS)
mod_forw <- stepAIC(empty_model, trace=TRUE, direction="forward", scope=horizonte)## Start: AIC=419.98
## y ~ 1
##
## Df Deviance AIC
## + cp 3 331.86 339.86
## + oldpeak 1 357.20 361.20
## + exang 1 359.54 363.54
## + thalach 1 360.90 364.90
## + slope 2 370.76 376.76
## + sex 1 393.93 397.93
## + age 1 402.54 406.54
## + restecg 2 407.84 413.84
## + trestbps 1 411.03 415.03
## + chol 1 415.78 419.78
## <none> 417.98 419.98
## + fbs 1 417.79 421.79
##
## Step: AIC=339.86
## y ~ cp
##
## Df Deviance AIC
## + oldpeak 1 294.12 304.12
## + slope 2 303.41 315.41
## + thalach 1 307.69 317.69
## + sex 1 309.30 319.30
## + exang 1 314.25 324.25
## + age 1 322.47 332.47
## + trestbps 1 324.93 334.93
## + restecg 2 327.12 339.12
## <none> 331.86 339.86
## + fbs 1 330.65 340.65
## + chol 1 330.79 340.79
##
## Step: AIC=304.12
## y ~ cp + oldpeak
##
## Df Deviance AIC
## + sex 1 274.45 286.45
## + thalach 1 282.37 294.37
## + exang 1 283.26 295.26
## + slope 2 283.76 297.76
## + age 1 288.71 300.71
## + trestbps 1 290.13 302.13
## <none> 294.12 304.12
## + restecg 2 290.61 304.61
## + fbs 1 292.67 304.67
## + chol 1 293.39 305.39
##
## Step: AIC=286.46
## y ~ cp + oldpeak + sex
##
## Df Deviance AIC
## + thalach 1 260.52 274.52
## + slope 2 260.05 276.05
## + age 1 264.51 278.51
## + exang 1 265.21 279.21
## + trestbps 1 268.29 282.29
## + chol 1 270.78 284.78
## <none> 274.45 286.45
## + restecg 2 270.61 286.61
## + fbs 1 273.80 287.80
##
## Step: AIC=274.52
## y ~ cp + oldpeak + sex + thalach
##
## Df Deviance AIC
## + trestbps 1 253.32 269.32
## + exang 1 255.08 271.08
## + slope 2 253.31 271.31
## + chol 1 255.34 271.34
## + age 1 256.59 272.59
## <none> 260.52 274.52
## + restecg 2 257.12 275.12
## + fbs 1 259.98 275.98
##
## Step: AIC=269.32
## y ~ cp + oldpeak + sex + thalach + trestbps
##
## Df Deviance AIC
## + slope 2 246.07 266.07
## + chol 1 249.01 267.01
## + exang 1 249.20 267.20
## <none> 253.32 269.32
## + age 1 251.60 269.60
## + restecg 2 250.56 270.56
## + fbs 1 253.25 271.25
##
## Step: AIC=266.07
## y ~ cp + oldpeak + sex + thalach + trestbps + slope
##
## Df Deviance AIC
## + chol 1 242.59 264.59
## + exang 1 243.08 265.08
## <none> 246.07 266.07
## + age 1 244.52 266.52
## + restecg 2 243.78 267.79
## + fbs 1 246.00 268.00
##
## Step: AIC=264.59
## y ~ cp + oldpeak + sex + thalach + trestbps + slope + chol
##
## Df Deviance AIC
## + exang 1 239.91 263.92
## <none> 242.59 264.59
## + age 1 241.58 265.58
## + fbs 1 242.53 266.52
## + restecg 2 241.19 267.19
##
## Step: AIC=263.92
## y ~ cp + oldpeak + sex + thalach + trestbps + slope + chol +
## exang
##
## Df Deviance AIC
## <none> 239.91 263.92
## + age 1 238.75 264.75
## + fbs 1 239.90 265.90
## + restecg 2 238.26 266.26mod_forw$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## y ~ 1
##
## Final Model:
## y ~ cp + oldpeak + sex + thalach + trestbps + slope + chol +
## exang
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 302 417.9821 419.9821
## 2 + cp 3 86.117585 299 331.8646 339.8646
## 3 + oldpeak 1 37.746319 298 294.1182 304.1182
## 4 + sex 1 19.662743 297 274.4555 286.4555
## 5 + thalach 1 13.932037 296 260.5235 274.5235
## 6 + trestbps 1 7.201410 295 253.3220 269.3220
## 7 + slope 2 7.255186 293 246.0669 266.0669
## 8 + chol 1 3.477537 292 242.5893 264.5893
## 9 + exang 1 2.674197 291 239.9151 263.9151summary(mod_forw)##
## Call:
## glm(formula = y ~ cp + oldpeak + sex + thalach + trestbps + slope +
## chol + exang, family = binomial(link = "logit"), data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.991830 2.144868 -2.327 0.019947 *
## cp2 0.837095 0.697391 1.200 0.230014
## cp3 0.318344 0.614118 0.518 0.604196
## cp4 2.313529 0.606605 3.814 0.000137 ***
## oldpeak 0.604296 0.197372 3.062 0.002201 **
## sex1 1.983987 0.418265 4.743 2.1e-06 ***
## thalach -0.023069 0.009281 -2.486 0.012936 *
## trestbps 0.021988 0.009529 2.308 0.021024 *
## slope2 0.782971 0.401647 1.949 0.051248 .
## slope3 -0.334461 0.860390 -0.389 0.697475
## chol 0.006059 0.003358 1.805 0.071152 .
## exang1 0.633063 0.385469 1.642 0.100524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 417.98 on 302 degrees of freedom
## Residual deviance: 239.92 on 291 degrees of freedom
## AIC: 263.92
##
## Number of Fisher Scoring iterations: 5# Seleccion BACKWARD ------------------------------------------------------
full_model <- glm(horizonte, data=datos, family=binomial(link='logit'))
mod_back <- stepAIC(full_model, trace=TRUE, direction="backward")## Start: AIC=269.27
## y ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach +
## exang + oldpeak + slope
##
## Df Deviance AIC
## - restecg 2 238.75 266.75
## - fbs 1 237.27 267.27
## - age 1 238.24 268.24
## - chol 1 239.11 269.11
## <none> 237.27 269.27
## - exang 1 240.28 270.28
## - slope 2 242.29 270.29
## - trestbps 1 240.86 270.86
## - thalach 1 241.15 271.15
## - oldpeak 1 247.55 277.55
## - sex 1 263.94 293.94
## - cp 3 270.15 296.15
##
## Step: AIC=266.75
## y ~ age + sex + cp + trestbps + chol + fbs + thalach + exang +
## oldpeak + slope
##
## Df Deviance AIC
## - fbs 1 238.75 264.75
## - age 1 239.90 265.90
## <none> 238.75 266.75
## - chol 1 241.34 267.34
## - exang 1 241.55 267.55
## - slope 2 244.05 268.05
## - trestbps 1 242.59 268.59
## - thalach 1 242.69 268.69
## - oldpeak 1 249.01 275.01
## - sex 1 266.09 292.09
## - cp 3 271.55 293.55
##
## Step: AIC=264.75
## y ~ age + sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
## Df Deviance AIC
## - age 1 239.91 263.92
## <none> 238.75 264.75
## - chol 1 241.34 265.34
## - exang 1 241.58 265.58
## - slope 2 244.06 266.06
## - trestbps 1 242.69 266.69
## - thalach 1 242.69 266.69
## - oldpeak 1 249.04 273.04
## - sex 1 266.57 290.57
## - cp 3 271.77 291.77
##
## Step: AIC=263.92
## y ~ sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
## Df Deviance AIC
## <none> 239.91 263.92
## - exang 1 242.59 264.59
## - chol 1 243.08 265.08
## - slope 2 245.35 265.35
## - trestbps 1 245.38 267.38
## - thalach 1 246.42 268.42
## - oldpeak 1 250.46 272.46
## - sex 1 266.98 288.98
## - cp 3 272.65 290.65mod_back$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## y ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach +
## exang + oldpeak + slope
##
## Final Model:
## y ~ sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 287 237.2655 269.2655
## 2 - restecg 2 1.4842382886 289 238.7498 266.7498
## 3 - fbs 1 0.0007679681 290 238.7505 264.7505
## 4 - age 1 1.1646061219 291 239.9151 263.9151summary(mod_back)##
## Call:
## glm(formula = y ~ sex + cp + trestbps + chol + thalach + exang +
## oldpeak + slope, family = binomial(link = "logit"), data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.991830 2.144868 -2.327 0.019947 *
## sex1 1.983987 0.418265 4.743 2.1e-06 ***
## cp2 0.837095 0.697391 1.200 0.230014
## cp3 0.318344 0.614118 0.518 0.604196
## cp4 2.313529 0.606605 3.814 0.000137 ***
## trestbps 0.021988 0.009529 2.308 0.021024 *
## chol 0.006059 0.003358 1.805 0.071152 .
## thalach -0.023069 0.009281 -2.486 0.012936 *
## exang1 0.633063 0.385469 1.642 0.100524
## oldpeak 0.604296 0.197372 3.062 0.002201 **
## slope2 0.782971 0.401647 1.949 0.051248 .
## slope3 -0.334461 0.860390 -0.389 0.697475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 417.98 on 302 degrees of freedom
## Residual deviance: 239.92 on 291 degrees of freedom
## AIC: 263.92
##
## Number of Fisher Scoring iterations: 5# Seleccion BOTH
mod_both <- stepAIC(full_model, trace=TRUE, direction="both")## Start: AIC=269.27
## y ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach +
## exang + oldpeak + slope
##
## Df Deviance AIC
## - restecg 2 238.75 266.75
## - fbs 1 237.27 267.27
## - age 1 238.24 268.24
## - chol 1 239.11 269.11
## <none> 237.27 269.27
## - exang 1 240.28 270.28
## - slope 2 242.29 270.29
## - trestbps 1 240.86 270.86
## - thalach 1 241.15 271.15
## - oldpeak 1 247.55 277.55
## - sex 1 263.94 293.94
## - cp 3 270.15 296.15
##
## Step: AIC=266.75
## y ~ age + sex + cp + trestbps + chol + fbs + thalach + exang +
## oldpeak + slope
##
## Df Deviance AIC
## - fbs 1 238.75 264.75
## - age 1 239.90 265.90
## <none> 238.75 266.75
## - chol 1 241.34 267.34
## - exang 1 241.55 267.55
## - slope 2 244.05 268.05
## - trestbps 1 242.59 268.59
## - thalach 1 242.69 268.69
## + restecg 2 237.27 269.27
## - oldpeak 1 249.01 275.01
## - sex 1 266.09 292.09
## - cp 3 271.55 293.55
##
## Step: AIC=264.75
## y ~ age + sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
## Df Deviance AIC
## - age 1 239.91 263.92
## <none> 238.75 264.75
## - chol 1 241.34 265.34
## - exang 1 241.58 265.58
## - slope 2 244.06 266.06
## - trestbps 1 242.69 266.69
## - thalach 1 242.69 266.69
## + fbs 1 238.75 266.75
## + restecg 2 237.27 267.27
## - oldpeak 1 249.04 273.04
## - sex 1 266.57 290.57
## - cp 3 271.77 291.77
##
## Step: AIC=263.92
## y ~ sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
## Df Deviance AIC
## <none> 239.91 263.92
## - exang 1 242.59 264.59
## + age 1 238.75 264.75
## - chol 1 243.08 265.08
## - slope 2 245.35 265.35
## + fbs 1 239.90 265.90
## + restecg 2 238.26 266.26
## - trestbps 1 245.38 267.38
## - thalach 1 246.42 268.42
## - oldpeak 1 250.46 272.46
## - sex 1 266.98 288.98
## - cp 3 272.65 290.65mod_both$anova## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## y ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach +
## exang + oldpeak + slope
##
## Final Model:
## y ~ sex + cp + trestbps + chol + thalach + exang + oldpeak +
## slope
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 287 237.2655 269.2655
## 2 - restecg 2 1.4842382886 289 238.7498 266.7498
## 3 - fbs 1 0.0007679681 290 238.7505 264.7505
## 4 - age 1 1.1646061219 291 239.9151 263.9151summary(mod_both)##
## Call:
## glm(formula = y ~ sex + cp + trestbps + chol + thalach + exang +
## oldpeak + slope, family = binomial(link = "logit"), data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.991830 2.144868 -2.327 0.019947 *
## sex1 1.983987 0.418265 4.743 2.1e-06 ***
## cp2 0.837095 0.697391 1.200 0.230014
## cp3 0.318344 0.614118 0.518 0.604196
## cp4 2.313529 0.606605 3.814 0.000137 ***
## trestbps 0.021988 0.009529 2.308 0.021024 *
## chol 0.006059 0.003358 1.805 0.071152 .
## thalach -0.023069 0.009281 -2.486 0.012936 *
## exang1 0.633063 0.385469 1.642 0.100524
## oldpeak 0.604296 0.197372 3.062 0.002201 **
## slope2 0.782971 0.401647 1.949 0.051248 .
## slope3 -0.334461 0.860390 -0.389 0.697475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 417.98 on 302 degrees of freedom
## Residual deviance: 239.92 on 291 degrees of freedom
## AIC: 263.92
##
## Number of Fisher Scoring iterations: 5# Nota: los modelos finales son el mismo, no tiene por quéOtra opción, a partir de cualquier modelo:
library(RcmdrMisc)## Loading required package: car## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## Loading required package: sandwichglm_intermedio <- glm(y ~age + sex + chol , data=datos, family=binomial(link='logit'))
modelo <- stepwise(glm_intermedio, direction='backward', criterion='AIC')##
## Direction: backward
## Criterion: AIC
##
## Start: AIC=376.59
## y ~ age + sex + chol
##
## Df Deviance AIC
## <none> 368.59 376.59
## - chol 1 372.31 378.31
## - age 1 386.96 392.96
## - sex 1 402.03 408.03modelo <- stepwise(glm_intermedio, direction='forward', criterion='AIC')##
## Direction: forward
## Criterion: AIC
##
## Start: AIC=419.98
## y ~ 1
##
## Df Deviance AIC
## + sex 1 393.93 397.93
## + age 1 402.54 406.54
## + chol 1 415.78 419.78
## <none> 417.98 419.98
##
## Step: AIC=397.93
## y ~ sex
##
## Df Deviance AIC
## + age 1 372.31 378.31
## + chol 1 386.96 392.96
## <none> 393.93 397.93
##
## Step: AIC=378.31
## y ~ sex + age
##
## Df Deviance AIC
## + chol 1 368.59 376.59
## <none> 372.31 378.31
##
## Step: AIC=376.59
## y ~ sex + age + cholmodelo <- stepwise(glm_intermedio, direction='forward/backward', criterion='AIC') # No tiene porque llegarse al mismo resultado##
## Direction: forward/backward
## Criterion: AIC
##
## Start: AIC=419.98
## y ~ 1
##
## Df Deviance AIC
## + sex 1 393.93 397.93
## + age 1 402.54 406.54
## + chol 1 415.78 419.78
## <none> 417.98 419.98
##
## Step: AIC=397.93
## y ~ sex
##
## Df Deviance AIC
## + age 1 372.31 378.31
## + chol 1 386.96 392.96
## <none> 393.93 397.93
## - sex 1 417.98 419.98
##
## Step: AIC=378.31
## y ~ sex + age
##
## Df Deviance AIC
## + chol 1 368.59 376.59
## <none> 372.31 378.31
## - age 1 393.93 397.93
## - sex 1 402.54 406.54
##
## Step: AIC=376.59
## y ~ sex + age + chol
##
## Df Deviance AIC
## <none> 368.59 376.59
## - chol 1 372.31 378.31
## - age 1 386.96 392.96
## - sex 1 402.03 408.03Valores predichos y matriz de confusión
Una de las principales aplicaciones de un modelo de regresión
logística es que puede ser empleada como una herramienta de
clasificación. La siguiente tabla resulta de cruzar la variable
respuesta con una variable dicotómica cuyos valores resultan de las
probabilidades estimadas en el modelo. Para definir esta variable hemos
de decidir un punto de corte z y comparar las probabilidades estimadas
con z, si son mayores, la variable vale 1 y si son menores 0. Una
posibilidad estándar es considerar z=0.5 , aunque podría ser cualquier
otro valor. 
En ningún caso, las tablas deben usarse para comparar modelos ya que depende en gran parte de la distribución de las probabilidades en la muestra en la que están basadas, el mismo modelo evaluado en dos muestras diferentes puede dar lugar a clasificaciones distintas. Si el objetivo del análisis no es la clasificación de sujetos, entonces las tablas de clasificación son sólo el complemento de otras medidas de ajuste.
Una de las medidas más utilizadas para valorar la calidad de la predicción el el área bajo la curva ROC (Receiver Operating Characteristic) para ello debemos definir primeramente algunos conceptos:
- Sensitividad: proporción de verdaderos 1s estimados como 1s (la probabilidad de predecir correctamente un 1): Ss = a/(a + c).
- Especificidad es la proporción de verdaderos 0s estimados como 0s (la probabilidad de predecir correctamente un 0): Sp = d/(b + d).
- Tasa de falsos positivos: proporción de verdaderos 0s estimados como 1s (la probabilidad de predecir incorrectamente un 0): F+ = b/(b + d).
- Tasa de falsos negativos es la proporción de verdaderos 1s estimados como 0s (la probabilidad de predecir incorrectamente un 1): F- = c/(a + c).
Cada vez que elegimos un punto de corte la sensitividad y especificidad cambian, si nuestro objetivo es buscar el punto de corte óptimo desde el punto de vista de la clasificación, entonces buscaremos aquel que maximiza la sensitividad y especificidad. La curva ROC es un gráfico de sensitividad frente a 1-especificidad (tasa de falsos positivos) para todos los puntos de corte. El área por debajo de esa curva puede tomar valores entre 0 y 1, y nos proporciona una medida de la habilidad del modelo para discriminar entre aquellos individuos que presentan la respuesta de interés frente a los que no la presentan.
En general:
- Si ROC <0.5 el modelo no ayuda a discriminar.
- Si 0.6 <ROC<0.8 el modelo discrimina de forma adecuada.
- Si 0.8 <ROC <0.9 el modelo discrimina de forma excelente.
- Si ROC >0.9 el modelo discrimina de forma excepcional.
prob_hat <- predict(mod_back, newdata=datos, type="response")
cut_off <- 0.5 # Este threshold puede variar, se podría construir curva ROC
y_hat <- ifelse(prob_hat >= cut_off, "presence", "absence")
y_hat <- as.factor(y_hat)
# Manual
conf_mat <- table(Prediccion=y_hat, Real=datos$y)
addmargins(conf_mat)## Real
## Prediccion absence presence Sum
## absence 144 33 177
## presence 20 106 126
## Sum 164 139 303# Con el paquete caret
caret::confusionMatrix(data=y_hat, reference=datos$y)## Confusion Matrix and Statistics
##
## Reference
## Prediction absence presence
## absence 144 33
## presence 20 106
##
## Accuracy : 0.8251
## 95% CI : (0.7775, 0.8661)
## No Information Rate : 0.5413
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.6452
##
## Mcnemar's Test P-Value : 0.09929
##
## Sensitivity : 0.8780
## Specificity : 0.7626
## Pos Pred Value : 0.8136
## Neg Pred Value : 0.8413
## Prevalence : 0.5413
## Detection Rate : 0.4752
## Detection Prevalence : 0.5842
## Balanced Accuracy : 0.8203
##
## 'Positive' Class : absence
## # Calculando otras medidas ------------------------------------------------
# McFadden’s R2 ranges from 0 to just under 1.
# Values close to 0 indicate that the model has no predictive power.
# values of 0.2 to 0.4 for McFadden’s R2 represent EXCELLENT fit
pscl::pR2(mod_back)["McFadden"]## fitting null model for pseudo-r2## McFadden
## 0.4260158# The importance of each predictor variable in the model by using the
# varImp function from the caret package
caret::varImp(mod_back)## Overall
## sex1 4.7433722
## cp2 1.2003239
## cp3 0.5183755
## cp4 3.8138983
## trestbps 2.3075581
## chol 1.8045067
## thalach 2.4855314
## exang1 1.6423187
## oldpeak 3.0617170
## slope2 1.9494005
## slope3 0.38873174. Ejercicio propuesto. Datos Challenger
El 28 de enero de 1986, el transbordador espacial Challenger explotó 73 segundos después del despegue, causando la muerte de los siete astronautas a bordo. Uno de los factores clave en el desastre fue el posible fallo de los O-rings en temperaturas frías.
Los ingenieros advirtieron que las temperaturas bajas podrían aumentar el riesgo de fallo de los O-rings. Sin embargo, la NASA ignoró la advertencia y el Challenger fue lanzado a 31°F, la temperatura más baja registrada en un lanzamiento de transbordador.
Explorarar la relación entre la temperatura y la probabilidad de fallo de los O-rings utilizando regresión logística. ¿Qué ocurre a 31°F?
Carga de Datos y Exploración Inicial
# Cargar librerías necesarias
library(tidyverse)
library(ggplot2)
library(ResourceSelection) # Para el test de Hosmer-Lemeshow## ResourceSelection 0.3-6 2023-06-27challenger_data <- read.table(file="C:/Users/Yolanda/OneDrive - Universidad de Valladolid/Docencia/MASTER25/Practicas/GLM/Datos/Challenger.txt",sep="\t",header=TRUE)
# Los datos contienen 23 lanzamientos previos al Challenger con información sobre la temperatura y fallos en los O-rings.
str(challenger_data)## 'data.frame': 23 obs. of 9 variables:
## $ flight : chr "1" "2" "3" "5" ...
## $ date : chr "12/04/81" "12/11/81" "22/03/82" "11/11/82" ...
## $ nfails.field : int 0 1 0 0 0 0 0 0 1 1 ...
## $ nfails.nozzle: int 0 0 0 0 2 0 0 0 1 1 ...
## $ fail.field : int 0 1 0 0 0 0 0 0 1 1 ...
## $ fail.nozzle : int 0 0 0 0 1 0 0 0 1 1 ...
## $ temp : int 66 70 69 68 67 72 73 70 57 63 ...
## $ pres.field : int 50 50 50 50 50 50 100 100 100 200 ...
## $ pres.nozzle : int 50 50 50 50 50 50 50 100 100 200 ...# Visualización inicial con eje X desde 30°F y línea en 31°F
ggplot(challenger_data, aes(x = temp, y = fail.field)) +
geom_point() +
geom_vline(xintercept = 31, linetype = "dashed", color = "gray") + # Línea en 31°F
scale_x_continuous(limits = c(30, max(challenger_data$temp))) + # Eje X desde 30°F
labs(title = "fallo de O-rings vs. Temperatura",
x = "Temperatura (°F)", y = "O-ring falló (1 = sí, 0 = no)") +
theme_minimal()
Ajuste del modelo de regresión logística
# Ajustar el modelo logístico
chall_glm <- glm(fail.field ~ temp, data = challenger_data, family = binomial)
# Resumen del modelo
summary(chall_glm)##
## Call:
## glm(formula = fail.field ~ temp, family = binomial, data = challenger_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.0429 7.3786 2.039 0.0415 *
## temp -0.2322 0.1082 -2.145 0.0320 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28.267 on 22 degrees of freedom
## Residual deviance: 20.315 on 21 degrees of freedom
## AIC: 24.315
##
## Number of Fisher Scoring iterations: 5# Obtener razones de odds
exp(coef(chall_glm))## (Intercept) temp
## 3.412315e+06 7.928171e-01#Cada 1°F de aumento reduce las odds de fallo en un factor de exp(β1) ≈ 0.79 (disminuye en un 21%).
#Cada 1°F de reducción aumenta las odds de fallo en un factor de exp(-β1) ≈ 1.26 (incremento del 26%).Predicción
# Crear datos para predicción
temp_seq <- data.frame(temp = seq(min(challenger_data$temp), max(challenger_data$temp), length.out = 100))
# Predicciones con el modelo
temp_seq$prob_fail <- predict(chall_glm, newdata = temp_seq, type = "response")
# Gráfico de regresión logística
ggplot(challenger_data, aes(x = temp, y = fail.field)) +
geom_point() +
geom_line(data = temp_seq, aes(x = temp, y = prob_fail), color = "red") +
geom_vline(xintercept = 31, linetype = "dashed", color = "gray") + # Temperatura del Challenger
labs(title = "Probabilidad de Fallo de O-rings vs. Temperatura",
x = "Temperatura (°F)", y = "Probabilidad de Fallo") +
theme_minimal()
# El punto de 31°F (línea gris) muestra que la probabilidad de fallo era extremadamente alta.Diagnóstico del modelo
En regresión logística, el test Homer-Lemeshow calculado sobre los datos una vez que las observaciones se han segmentado en grupos basados en probabilidades predichas similares. Este test examina si las proporciones observadas de eventos son similares a las probabilidades predichas de ocurrencia en subgrupos del conjunto de datos, y lo hace con una prueba de chi cuadrado de Pearson.
La hipótesis nula sostiene que el modelo se ajusta a los datos y en el siguiente ejemplo rechazamos H0. No obstante, siempre preferimos el enfoque de la deviance.
#library(ResourceSelection)
# Test de bondad de ajuste de Hosmer-Lemeshow
#hoslem.test(challenger_data$fail_binary, fitted(chall_glm), g = 5)Probabilidad de fallo para Challenger
El Challenger fue lanzado a 31°F. Calculamos la probabilidad de fallo con nuestro modelo.
# Predicción para 31°F
prob_31F <- predict(chall_glm, newdata = data.frame(temp = 31), type = "response")
prob_31F## 1
## 0.9996088# Probabilidad de fallo ≈ 99.96%.
# La decisión de lanzar a esta temperatura era extremadamente arriesgada.
# Si se hubiera esperado hasta 53°F, la probabilidad de fallo hubiera bajado drásticamente.
# Comparación con 53°F
prob_53F <- predict(chall_glm, newdata = data.frame(temp = 53), type = "response")
prob_53F## 1
## 0.9392478# Factor de reducción
exp(coef(chall_glm)["temp"] * (53 - 31))## temp
## 0.006050706#Esperar hasta 53°F habría reducido las odds de fallo en un 99.4%.5. Ejercicio propuesto. Datos titanic
El transatlántico británico RMS Titanic se hundió después de chocar con un iceberg en su viaje inaugural. La tragedia se saldó con la muerte de 1502 de los 2224 pasajeros y tripulantes. En el siguiente enlace https://web.stanford.edu/class/archive/cs/cs109/cs109.1166/stuff/titanic.csv puedes descargar los datos relativos a los pasajeros, como por ejemplo si sobrevivieron, la clase en la que viajaban, el número de hijos a bordo etc. La descripción de las variables aparece a continuación:
- survival - Survival (0 = No; 1 = Yes)
- class - Passenger Class (1 = 1st; 2 = 2nd; 3 = 3rd)
- name - Name
- sex - Sex
- age - Age
- sibsp - Number of Siblings/Spouses Aboard
- parch - Number of Parents/Children Aboard
- ticket - Ticket Number
- fare - Passenger Fare
- cabin - Cabin
- embarked - Port of Embarkation (C = Cherbourg; Q = Queenstown; S = Southampton)
- boat - Lifeboat (if survived)
- body - Body number (if did not survive and body was recovered)
Si bien hubo algún elemento de suerte involucrado en la supervivencia, parece que algunos grupos de personas tenían más probabilidades de sobrevivir que otros. Analice varios modelos de regressión logística y seleccione uno de ellos para modelizar la supervivencia de los pasajeros. Para el modelo seleccionado incluya la bondad del ajuste, diagnosis del modelo y las principales conclusiones que se puedan extraer.
# Cargar datos
titanic <- read.csv("https://web.stanford.edu/class/archive/cs/cs109/cs109.1166/stuff/titanic.csv")
titanic <- as.data.frame(titanic)
str(titanic)## 'data.frame': 887 obs. of 8 variables:
## $ Survived : int 0 1 1 1 0 0 0 0 1 1 ...
## $ Pclass : int 3 1 3 1 3 3 1 3 3 2 ...
## $ Name : chr "Mr. Owen Harris Braund" "Mrs. John Bradley (Florence Briggs Thayer) Cumings" "Miss. Laina Heikkinen" "Mrs. Jacques Heath (Lily May Peel) Futrelle" ...
## $ Sex : chr "male" "female" "female" "female" ...
## $ Age : num 22 38 26 35 35 27 54 2 27 14 ...
## $ Siblings.Spouses.Aboard: int 1 1 0 1 0 0 0 3 0 1 ...
## $ Parents.Children.Aboard: int 0 0 0 0 0 0 0 1 2 0 ...
## $ Fare : num 7.25 71.28 7.92 53.1 8.05 ...titanic %>% summary()## Survived Pclass Name Sex
## Min. :0.0000 Min. :1.000 Length:887 Length:887
## 1st Qu.:0.0000 1st Qu.:2.000 Class :character Class :character
## Median :0.0000 Median :3.000 Mode :character Mode :character
## Mean :0.3856 Mean :2.306
## 3rd Qu.:1.0000 3rd Qu.:3.000
## Max. :1.0000 Max. :3.000
## Age Siblings.Spouses.Aboard Parents.Children.Aboard
## Min. : 0.42 Min. :0.0000 Min. :0.0000
## 1st Qu.:20.25 1st Qu.:0.0000 1st Qu.:0.0000
## Median :28.00 Median :0.0000 Median :0.0000
## Mean :29.47 Mean :0.5254 Mean :0.3833
## 3rd Qu.:38.00 3rd Qu.:1.0000 3rd Qu.:0.0000
## Max. :80.00 Max. :8.0000 Max. :6.0000
## Fare
## Min. : 0.000
## 1st Qu.: 7.925
## Median : 14.454
## Mean : 32.305
## 3rd Qu.: 31.137
## Max. :512.329