Evaluación de supuestos en regresión logística

1 Introducción

La regresión logística es uno de los métodos más conocidos para predecir el valor de una variable respuesta \(Y\) de tipo categórico, en función de un conjunto de covariables o variables predictoras. Algunos ejemplos del uso de este tipo de regresión son:

• \(Y=0\) ó \(Y=1\).

• \(Y=\) no sobrevive o \(Y=\) sobrevive

• \(Y=\) no paga el crédito o \(Y=\) si paga el crédito

El objetivo de este tipo de modelos es estimar la probabilidad del evento de interés:

\[P(Y=1|X=x)\]

Dado un vector \(X=x\) de covariables o variables predictoras, considerando que \(Y\) tiene solo dos posible valores, se tiene que \(Y \sim Bernoulli(p)\), donde \(P(Y=1)=p\) y \(P(Y=0)=1-p\).

Con un modelo de regresión logística es posible modelar \(P(Y=1|X=x)\). Usando la definción de la función logit se tiene que:

\[P(Y=1|X=x)=logit^-1(\beta^Tx)=\frac{exp(\beta^Tx)}{1+exp(\beta^Tx)}\]

2 Regresión logística aplicada a Titanic

Los datos a analizar corresponden al conjunto de datos del Titanic, los cuales contienen información acerca de los pasajeros que viajaban en la embarcación. Los campos de la base de datos se describen a continuación:

# Retirar columnas innecesarias
titanic <- titanic[,-c(3)]

# Renombrar columnas por facilidad
colnames(titanic) <- c("Sobreviviente", "Clase", "Sexo", "Edad",
                       "Hermanos", "Padres", "Tarifa")

# Codificar las columnas
titanic$Sobreviviente = factor(titanic$Sobreviviente)
titanic$Clase = factor(titanic$Clase)
titanic$Sexo = factor(titanic$Sexo)

• Sobreviviente: es un variable categórica que representa si el pasajero sobrevive o no.

• Clase: hace referencia a la clase en la que viajaba el pasajero, posee tres categorías; siendo 1 la clase más alta, es decir los pasajeros más adinerados.

• Sexo: representa el sexo de los pasajeros, femenino o masculino.

• Edad: variable numérica que contiene la información de la edad de los pasajeros.

• Tarifa: es una variable numérica con el precio que cada pasajero pagó para viajar en la embarcación.

• Hermanos: es una variable numérica discreta que contiene información sobre el número de hermanos o cónyuges a bordo.

• Padres: es una variable numérica discreta que contiene información sobre el número de padres e hijos a bordo.

2.1 Análisis descriptivo

La primera gráfica a analizar es el comportamiento de sobrevivencia y muerte.

ggplot(data = titanic, aes(Sobreviviente, fill = Sobreviviente)) +
  geom_bar(position="dodge") +
  ggtitle("Sobrevivientes") +
  ylab("Cantidad") + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position='top')

La base de datos utilizada en este análisis contiene información de \(887\) pasajeros de los cuales \(545\) fallecieron como se observa en la gráfica anterior y los restantes \(342\) lograron sobrevivir. Adicionalmente se concluye que no existe un problema grave de desbalanceo de clases que pueda perjudicar los modelos que se realizan en las secciones posteriores.

Las siguientes gráficas representan la catidad de pasajeros que sobrevivieron y muerion por clase y por sexo.

# Sobrevivientes por Clase
clase <- ggplot(data = titanic, aes(Clase, fill = Sobreviviente)) +
  geom_bar(position="dodge") +
  ggtitle("Sobrevivientes por Clase") +
  ylab("Cantidad") + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position='top')

# Sobrevivientes por Sexo
sexo <- ggplot(data = titanic, aes(Sexo, fill = Sobreviviente)) + 
  geom_bar(position = "dodge") + 
  ylab("Cantidad") + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Sexo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme(legend.position='top')

plot_grid(clase, sexo)

De la gráfica izquierda es posible identificar con claridad que la clase a la que pertenecia cada pasajero es un factor fundamental para la sobreviviencia o no de las personas, a la clase 3 pertenecian \(487\) pasajeros con bajos recursos económico y por ende presenta la mayor cantidad de personas fallecidas \(368\), mientras que en la clase 1 se ubicaban \(216\) pasajeros de los cuales lograron sobrevivir \(136\) y fallecieron \(80\) personas. Finalmente, para la clase 2 la cantidad de pasajeros fallecidos y sobrevivientes es muy similar, \(97\) y \(87\) respectivamente.

La gráfica de la derecha presenta la sobreviviencia de los pasajeros por sexo, de esta figura también es posible concluir con claridad que el sexo es un factor determinante en la sobrevivencia de las personas a bordo. De las \(314\) mujeres en la embarcación lograron sobrevivir \(233\), mientras que de los \(573\) hombres a bordo únicamente sobreviven \(109\).

# Sobrevivientes por edad
edad <- ggplot(data = titanic, aes(Sobreviviente, Edad)) + 
  geom_boxplot(aes(fill=Sobreviviente), alpha=0.9) + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Edad") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))+
  theme(legend.position='top')

#Sobrevivientes por Tarifa
tarifa <- ggplot(data = titanic, aes(Sobreviviente, Tarifa)) + 
  geom_boxplot(aes(fill=Sobreviviente), alpha=0.9) + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Tarifa") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))+
  theme(legend.position='top')
plot_grid(edad, tarifa)

Las gráficas anteriores presentan el comportamiento de la sobrevivencia por las variables numéricas edad y tarifa, incialmente con la variable edad se observan algunos valores atípicos, pero en general la edad no es una varible influyente en la sobreviviencia o no de los pasajeros. Respecto a la tarifa se observa una leve diferencia en la mediana de los sobrevivientes y fallecidos, sin embargo no es un diferencia grande, también posee algunos puntos atípicos, pero no son de gran relevancia.

A continuación se presenta un análisis exploratiorio de las variables categóricas en función de las variables numéricas que se tienen disponibles.

clase_edad <- ggplot(data = titanic, aes(Clase, Edad)) + 
  geom_boxplot(aes(fill=Sobreviviente), alpha=0.9) + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Edad y Clase") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))+
  theme(legend.position='top')
sexo_edad <- ggplot(data = titanic, aes(Sexo, Edad)) + 
  geom_boxplot(aes(fill=Sobreviviente), alpha=0.9) + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Edad y Sexo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))+
  theme(legend.position='top')
plot_grid(clase_edad, sexo_edad)

De la figura izquierda se observa que las personas que viajaban en la clase 1 tenian una edad más avanzada que las personas de la clase 2 y 3, pero dentro de cada clase las edades de sobrevivencia y fallecimiento son muy similares, la variabilidad más alta se presenta en la clase 1. Las edades de los psajeros que sobrevivieron estaban al rededor de los \(36\), \(28\) y \(22\) años para las clases 1, 2 y 3 respectivamente. Adicional a esto se presentan varios valores atípicos en las edades de las personas que no sobrevivieron y que pertenecían a la clase 3.

La figura de la derecha permite evidenciar que tanto mujeres como hombres sobrevivientes tenian edades muy similares, sin embargo, no hay gran diferencia entre las edades de cada sexo ya que los boxplot se traslapan y sus medianas no estan muy alejadas.

Teniendo en cuenta que la clase es un factor determinante para que cada pasajero lograra sobrevivir se realizan gráficos para analizar la sobreviviencia por sexo para cada una de las clases.

ggplot(data = titanic, aes(Sexo, fill = Sobreviviente)) + 
  geom_bar(position = "dodge") +
  facet_wrap(~Clase) + 
  ylab("Cantidad") + 
  scale_fill_manual(values=c("#87ceeb", "#ffd700")) +
  ggtitle("Sobrevivientes por Clase y Sexo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme(legend.position='top')

De la gráfica anterior se observa que en todas las clases la mayor parte de los sobrevivientes eran mujeres, algo de esperarse ya que eran prioridad al momento de abandonar la embarcación, algo también de resaltar es que los hombres fallecidos en la clase 3 son alrededor de 3 veces más de los que muerieron en las demás clases.

2.2 Ajuste del modelo

De manera inicial se realiza una partición para el conjunto de datos, tomando un 80% de las observaciones como datos de entrenamiento y un 20% restante como datos de prueba.Se realiza un proceso de selección de variables y se determina que las variables y no son significativas. El resumen para el modelo completo, ajustado con el conjunto de datos de entrenamiento se presenta a continuación:

#  Partición del conjunto de datos
split <- sample(c(rep(0, 0.8 * nrow(titanic)), rep(1, 0.2 * nrow(titanic))))
train <- titanic[split == 0, ]  
test <- titanic[split == 1, ]  

# Ajuste del modelo completo
fit <- glm(Sobreviviente ~ Clase + Sexo + Edad + Hermanos, data = train,
           family=binomial(logit))
summary(fit)

## 
## Call:
## glm(formula = Sobreviviente ~ Clase + Sexo + Edad + Hermanos, 
##     family = binomial(logit), data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7766  -0.6345  -0.4085   0.6645   2.4152  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.361206   0.463203   9.415  < 2e-16 ***
## Clase2      -1.315094   0.296503  -4.435 9.19e-06 ***
## Clase3      -2.414706   0.278987  -8.655  < 2e-16 ***
## Sexomale    -2.615949   0.214723 -12.183  < 2e-16 ***
## Edad        -0.048702   0.008662  -5.622 1.88e-08 ***
## Hermanos    -0.430459   0.120876  -3.561 0.000369 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 954.03  on 709  degrees of freedom
## Residual deviance: 649.55  on 704  degrees of freedom
## AIC: 661.55
## 
## Number of Fisher Scoring iterations: 5

A partir del resumen para el modelo completo es posible que todas las variables utilizadas como predictoras son significativas para predecir la probabilidad de sobrevivencia de los pasajeros.

Además, se presentan algunas de las probabilidades estimadas con el conjunto de datos de entrenamiento:

# Probabilidades estimadas
train$Probabilidades_estimadas <- fit$fitted.values

kable(train[c(1:4,6), c("Sobreviviente", "Clase", "Sexo", "Edad", "Hermanos", "Padres",
             "Probabilidades_estimadas")], digits = 2)

	Sobreviviente	Clase	Sexo	Edad	Hermanos	Padres	Probabilidades_estimadas
1	0	3	male	22	1	0	0.10
2	1	1	female	38	1	0	0.89
4	1	1	female	35	1	0	0.90
5	0	3	male	35	0	0	0.09
10	1	2	female	14	1	0	0.87

Por último, es posible observar la matriz de confusión del modelo, junto con las medidas de evaluación para modelos de clasificación:

confusion <- function(mod, data=NULL, cutoff=0.5) {
if (is.null(data)) y <- mod$y
else {
name_y <- deparse(formula(mod)[2])
name_y <- substr(x=name_y, start=1, stop=nchar(name_y)-2)
y <- data[[name_y]]
}
probability <- predict(object=mod, newdata=data, type='response')
prediction <- ifelse(probability <= cutoff, 0, 1)
confusion <- table(Prediction=prediction, Real=y)
TP <- confusion[1, 1]; FP <- confusion[1, 2]
FN <- confusion[2, 1]; TN <- confusion[2, 2]
list(confusion=confusion,
sensitivity=TP / (TP + FN),
specificity=TN /( TN + FP),
accuracy=(TP + TN) / sum(confusion),
false_discovery_rate=FP / (FP + TP))
}
con <- confusion(mod = fit, data = train, cutoff = 0.5)
con

## $confusion
##           Real
## Prediction   0   1
##          0 369  85
##          1  59 197
## 
## $sensitivity
## [1] 0.8621495
## 
## $specificity
## [1] 0.6985816
## 
## $accuracy
## [1] 0.7971831
## 
## $false_discovery_rate
## [1] 0.1872247

• De acuerdo al accuracy del modelo, se obtiene que el modelo clasifica de manera correcta un 80.7% de las observaciones, aproximadamente.

• De acuerdo con la especificididad, el modelo predice valores negativos un 70% de las veces.

• Según la sensibilidad, el modelo tiene una probabilidad del 87.27% de predecir correctamente un individuo clasificado como sobreviviente.

3 Evaluación de supuestos en regresión logística

3.1 Linealidad

Uno de los supuestos más importantes en regresión logística es que la relación entre el \(logit\) o log-odds de la variable respuesta y cada variable predictora o variable independiente es lineal; este supuesto se verifica únicamente para las variables numéricas continuas que se tengan en el modelo.

El \(logit\) es el logaritmo del \(odds \ ratio\), en donde \(p=\) Probabilidad de éxito, que para el conjunto de datos, se refiere a la probabilidad de sobrevivencia del pasajero al hundimiento del barco.

\[logit(p)=log\left(\frac{p}{1-p}\right)\]

3.1.1 Test de Box Tidwell

Para verificar este supuesto, es posible usar el Test Box-Tidwell, el cual se realiza por medio de la función boxTidwell del paquete car. El test se plantea tomando como hipótesis nula el cumplimineto del supuesto de linealidad. Según esto, se puede apreciar el siguiente resultado:

logodds <- fit$linear.predictors
boxTidwell(logodds ~ train$Edad)

##  MLE of lambda Score Statistic (z) Pr(>|z|)
##        0.97816              0.0251     0.98
## 
## iterations =  3

La salida anterior presenta el valor p asociado al estadístico de prueba del test de Box Tidwell, a partir del cual, no se rechaza el supuesto de linealidad entre el \(logit\) y la variable Edad. Además, es posible apreciar el gráfico de la edad vs. el logaritmo de las probabilidades estimadas:

logEdad <- data.frame(logodds, Edad = train$Edad)
ggplot(data = logEdad, aes(x = Edad, y = logodds)) + 
  geom_point() +
  ggtitle("Edad vs Logaritmo del Odds") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

El diagrama de dispersión no parece tener un patrón no lineal fuerte, por lo que se puede concluir que no se observa una violación al supuesto de linealidad.

Nota: Para la comprobación de este supuesto en Python, no hay una función integrada de la misma manera que en R, por lo que es posible verificar este supuesto agregando los términos de interacción entre las variables continuas y su correspondiente transformación logarítmica dentro del modelo, o en su defecto, por medio del diagrama de dispersión entre las variables continuas y el logaritmo de las probabilidades estimadas.

3.2 Independencia entre las observaciones

La independencia entre observaciones es uno de los supuestos más importantes. Suele ser el resultado del proceso de recolección de los datos, por lo que usualmente no es fácil detectarlo mediante la evaluación de residuales del modelo. La mejor forma de evaluar este supuesto usualmente es comprender el proceso mediante el cual se recopilaron los datos. Sin embargo, si los datos no fueron recopilados a lo largo del tiempo, este supuesto puede ser verificado por medio del gráfico de los residuos vs. los residuos rezagados un período de tiempo. En un escenario ideal, no se debe observar ningún tipo de patrón aparente.

En el caso de datos espaciales, la independencia puede probarse graficando los residuos vs. variables explicativas espaciales, como latitud o longitud. De igual forma, para la verificación del supuesto de independencia, no se debe observar ningún tipo de patrón aparente en estos gráficos.

Esta suposición de independencia se cumple automáticamente para el conjunto de datos de Titanic. Suposición que sería más preocupante en datos de series de tiempo, donde la autocorrelación puede ser un problema. No obstante, es pñosible verificar la independencia de las observaciones para estos datos gráficamente. Aquí, la ‘variable de tiempo’ es el orden de las observaciones, es decir, los números índice.

En particular, podemos crear la gráfica de los residuales de desviación del modelo logit contra los números índice de las observaciones.

Indice <- seq(1,710,1)
Residuales <- fit$residuals
residuales <- data.frame(Indice, Residuales)
ggplot(data = residuales, aes(x = Indice, y = Residuales)) +
  geom_point() +
  ggtitle("Residuales del modelo ajustado") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

El gráfico de dispersión permite observar más residuos positivos que negativos, sin embargo, los puntos se distribuyen de manera homogénea y no se observa ningún patrón aparente que pudiera indicar un grado de dependencia entre las observaciones, por lo que se puede concluir que se cumple el supuesto de independencia entre observaciones para el conjunto de datos de Titanic, en particular.

3.3 Ausencia de multicolinealidad

Se dice que hay presencia de multicolinealidad cuando dos o más variables dentro del conjunto de datos mantienen una relación lineal. Existe una colinealidad perfecta cuando alguna covariable puede ser representada como combinación lineal de una o más covariables presentes en los datos. Cuando esto sucede, debe eliminarse aquella covariable que genera un grado de dependencia con las demás.

Usualmente es posible trabajar con un grado de correlación moderado, ya que cuando la correlación entre variables es muy alta, se genera un incremento en los errores estándar, lo cual hace que los coeficientes estimados no sean confiables, y en consecuencia, las estimaciones sean poco creíbles.

La metodología a seguir para probar este supuesto consiste en evaluar inicialmente los coeficientes de correlación entre las variables explicativas. Luego, es necesario verificar colinealidad por medio de los valores de tolerancia o los factores infladores de varianza (VIF). Finalmente, si a partir de estas medidas hay sospecha de multicolinealidad, esto puede ser confirmado por medio del índice de condición, los valores propios y las proporciones de varianza.

3.3.1 Correlación

Inicialmente se presenta un mapa de calor con la correlación existente entre las dos variables numéricas que fueron consideradas en el ajuste del modelo.

cormat <- round(cor(titanic[,c(4,5)]),2)
library(reshape2)
melted_cormat <- melt(cormat)
library(ggplot2)
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile()

El anterior mapa de calor permite concluir que existe una correlación negativa entre las variables Edad y Hermanos para los datos que se estan analizando, dicha correlación esta alrededor de \(-0.25\).

Para verificar este valor de correlación se realiza el siguiente gráfico, donde se presenta el histograma para cada variable, valor de correlación y gráfico de dispersión.

correlaciones <- cor(titanic[,c("Edad", "Tarifa")])

my_data <- titanic[, c(4,5)]
chart.Correlation(my_data, histogram=TRUE, pch=19)

Se observa que la variable edad no posee una distribución simétrica, ya que es sesgada a la derecha. Además,la variable hermanos es una variable numérica discreta. El gráfico de dispersión no muestra patrones claros de correlación, sin embargo, la correlación de estas dos variables es de \(-0.3\), lo cual no es preocupante en términos de violación del supuesto.

3.3.2 Factores infladores de varianza o VIF

Los factores infladores de varianza o VIF determinan el grado de relación entre variables independientes. Además, vienen determinados por el ajuste de una regresión entre ambas variables. La tolerancia está definida como el recíproco de los valores VIF, por lo que se puede verificar el supuesto con cualquiera de estas dos medidas.

A continuación se presentan los valores para los estadísticos VIF para las variables numéricas que se tienen en el modelo.

kable(vif(fit), digits = 2)

	GVIF	Df	GVIF^(1/(2*Df))
Clase	1.46	2	1.10
Sexo	1.13	1	1.06
Edad	1.52	1	1.23
Hermanos	1.20	1	1.10

Sería problemático obtener valores VIF mayores a 5 o 10, ya que esto generaría un sesgo en el modelo, por lo que sería viable descartar alguna de las variables predictoras. Para el conjunto de datos considerado, los valores en la tabla parecen indicar que no hay un problema evidente de multicolinealidad.

4 Conclusiones

El conjunto de datos seleccionado contiene información acerca de las personas que viajaban en el barco al momento de la tragedia, después de ajustado el modelo y validado los tres supuestos principales para este tipo de modelos de clasificación se determina que el conjunto de datos cumple con todos los supuestos, por lo tanto, es un modelo adecuado para predecir si un pasajero del Titanic sobrevive o no al hundimiento del barco.

5 Referencias

• Box, G. E. P. and Tidwell, P. W. (1962) Transformation of the independent variables. Technometrics 4, 531-550. Fox, J. (1997) Applied Regression, Linear Models, and Related Methods. Sage.

• Dunn & Smyth, Generalized Linear Models With Examples in R, Springer (2018).

• Hernandez, Freddy, and Mazo Mauricio. 2020. Modelos de regresión con R. https://fhernanb.github.io/libro_regresion/.

• Jaén, M. E. F. (n.d.). Feir 45: Regresión Logística. Retrieved May 22, 2022, from https://gauss.inf.um.es/feir/45/#26_supuestos_del_modelo.

• Senaviratna, N. A. M. R., & Cooray, T. M. J. A. (2019). Diagnosing multicollinearity of logistic regression model. Asian Journal of Probability and Statistics, 5(2), 1-9.

• Thilakshasilva. (2017, December 14). Predicting titanic survival using five algorithms. Kaggle. Retrieved May 22, 2022, from https://www.kaggle.com/code/thilakshasilva/predicting-titanic-survival-using-five-algorithms#introduction.