6 Regresión logística

Este capítulo introduce la regresión logística como el método más sencillo para crear modelos predictivos en problemas de clasificación que es el principal objetivo de la asignatura. Se cubrirán los siguientes temas:

• Conocer la función logística
• Cómo interpretar los coeficientes de los modelos
• Cómo evaluar la capacidad predictiva de un modelo
• Cómo interpretar variables
• Ilustrar un ejemplo de análisis completo
• Aprender a hacer nomogramas fijos y dinámicos

Hasta ahora, nuestra variable de resultado era continua. Pero si la variable de resultado es binaria (0/1, “No”/“Sí”), entonces nos enfrentamos a un problema de clasificación. El objetivo de la clasificación es crear un modelo capaz de clasificar el resultado — y, cuando se usa el modelo para la predicción, nuevas observaciones— en una de dos categorías. La regresión logística se introduce en el contexto de la epidemiología como un modelo de regresión que extiende el modelo lineal cuando nuestra variable respuesta es binaria, pero tambié es, probablemente, el método estadístico más utilizado para la clasificación y el más sencillo. Una de las grandes ventajas de estos modelos respecto a otros que veremo más adelante es que este método produce un modelo de probabilidad para nuestra variable resultado. En otras palabras, los valores ajustados en un modelo logístico o logit no son binarios sino que son probabilidades que representan la probabilidad de que el resultado pertenezca a una de nuestras dos categorías.

Desafortunadamente, debemos afrontar nuevas complicaciones cuando trabajamos con regresión logística, lo que hace que estos modelos sean inherentemente más difíciles de interpretar que los modelos lineales. Las complicaciones surgen del hecho de que con la regresión logística modelamos la probabilidad de que \(y\) = 1, y la probabilidad siempre se escala entre 0 y 1. Pero el predictor lineal, \(X_i \beta\), oscila entre \(\pm \infty\) (donde \(X\) representa todos los predictores del modelo). Esta diferencia de escala requiere transformar la variable de resultado, lo cual se logra con la función logit:

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

La función logit asigna el rango del resultado (0,1) al rango del predictor lineal \((-\infty, +\infty)\). El resultado transformado, \(\text{logit} (x)\), se expresa en logaritmos de probabilidades (\(\frac{x}{1-x}\) se conoce como probabilidades del resultado - razón de odds en inglés - momios en castellano). Así que el modelo también se puede escribir como:

\[\text{Pr}(y_i = 1) = p_i\]

\[\text{logit}(p_i) = X_i\beta\]

Las probabilidades logarítmicas (e.g. el log-odds) no tienen interpretación (que no sea el signo y la magnitud) y deben transformarse nuevamente en cantidades interpretables, ya sea en probabilidades, usando el logit inverso, o en razones de probabilidades, mediante el uso de la función exponencial. Discutimos ambas transformaciones a continuación.

6.1 La función logit inversa

El modelo logístico se puede escribir, alternativamente, usando el logit inverso:

\[ \operatorname{Pr}(y_i = 1 | X_i) = \operatorname{logit}^{-1}(X_i \beta) \]

donde \(y_i\) es la respuesta binaria, \(\operatorname{logit}^{- 1}\) es la función logit inversa y \(X_i \beta\) es el predictor lineal. Podemos interpretar esta formulación diciendo que la probabilidad de que \(y = 1\) es igual al logit inverso del predictor lineal \((X_i, \ beta)\). Por lo tanto, podemos expresar los valores ajustados del modelo logístico y los coeficientes como probabilidades utilizando la transformación logit inversa. Pero, ¿qué es exactamente el logit inverso? Pues es:

\[\operatorname{logit}^{-1}(x) = \frac{e^{x}}{1 + e^{x}}\]

donde \(e\) es la función exponencial.

Podemos tener una idea de cómo la función logit inversa transforma el predictor lineal mediante una gráfica. Aquí usamos un rango arbitrario de valores de x en (-6, 6) para demostrar la transformación.

x <- seq(-6, 6, .01)
y <- exp(x)/(1 + exp(x))
ggplot(data.frame(x = x, y = y), aes(x, y)) +
  geom_line() + 
  ylab(expression(paste(logit^-1,"(x)"))) +
  ggtitle(expression(paste("y = ", logit^-1,"(x)")))

Los valores \(x\), que van de -6 a 6, son comprimidos por la función logit inversa en el rango 0-1. El logit inverso es curvo, por lo que la diferencia esperada en \(y\) correspondiente a una diferencia fija en \(x\) no es constante. A valores bajos y valores altos de \(x\), un cambio de unidad corresponde a un cambio muy pequeño en \(y\), mientras que en la mitad de la curva un pequeño cambio en \(x\) corresponde a un cambio relativamente grande en \(y\). En la regresión lineal, la diferencia esperada en \(y\) correspondiente a una diferencia fija en \(x\) es, por el contrario, constante. Por lo tanto, cuando interpretamos los resultados logísticos debemos elegir en qué parte de la curva queremos evaluar la probabilidad del resultado, dado el modelo.

6.2 Ejemplo de regresión logística

Ilustremos estos conceptos utilizando el conjunto de datos “Default” del ISLR. Este conjunto de datos simulado contiene una variable binaria que representa el incumplimiento en los pagos de la tarjeta de crédito (variable “default”), que modelaremos como una función de la variable “balance” (la cantidad de deuda que tiene la tarjeta) y la variable “income” (ingresos). Primero visualizaremos cómo es esta relación.

library(ISLR)
data(Default)
str(Default)

'data.frame':   10000 obs. of  4 variables:
 $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
 $ balance: num  730 817 1074 529 786 ...
 $ income : num  44362 12106 31767 35704 38463 ...

ggplot(Default, aes(x = balance, y = income, col = default)) +
  geom_point(alpha = .4) +
  ggtitle("Balance vs. Income by Default")

grid.arrange(
ggplot(Default, aes(default, balance)) + 
  geom_boxplot() +
  ggtitle("Balance by Default") + 
  ylab("balance"),
ggplot(Default, aes(default, income)) + 
  geom_boxplot() +
  ggtitle("Income by Default") + 
  ylab("income"),
ncol = 2)

Claramente, los valores altos de saldo están asociados con el incumplimiento en todos los niveles de ingresos. Esto sugiere que los ingresos en realidad no son un fuerte predictor de incumplimiento, en comparación con el saldo, que es exactamente lo que vemos en los diagramas de cajas.

Exploremos estas relaciones usando la regresión logística. En R ajustamos un modelo logístico usando la función glm () con family = binomial.⁶ Centraremos y escalaremos las variables para facilitar la interpretación.

glm(default ~ balance + income, 
      data = Default, 
      family = binomial) %>%
    standardize %>%
    display

glm(formula = default ~ z.balance + z.income, family = binomial, 
    data = Default)
            coef.est coef.se
(Intercept) -6.13     0.19  
z.balance    5.46     0.22  
z.income     0.56     0.13  
---
  n = 10000, k = 3
  residual deviance = 1579.0, null deviance = 2920.6 (difference = 1341.7)

NOTA: Se puede apreciar la ventaja del uso de tidyverse (pipe) - no se necesita crear las variables estandarizadas, ni guardar el resultado para luego hacer un print

Interpretamos esta salida exactamente como lo haríamos para un modelo lineal con un predictor centrado y escalado:

• Intercept: -6,13 representa las probabilidades logarítmicas (log odds) de incumplimiento cuando el saldo es promedio (835.37) y el ingreso es promedio (3.351698^{4} ). (Promedio porque las variables se han centrado).
• z.balance: 5.46 representa el cambio predicho en las probabilidades logarítmicas de incumplimiento asociado con un aumento de 1 unidad en el saldo (z.balance), manteniendo constante el ingreso (z.income). Un aumento de 1 unidad en el saldo (z.balance) es equivalente a un aumento de 2 desviaciones estándar en el saldo (balance) (967.43). Este coeficiente es estadísticamente significativo ya que 5.46 - 2 x .22 > 0 (el IC del 95% que no contiene 0 indica significación estadística).
• z.income: .56 representa el cambio predicho en las probabilidades logarítmicas (log odds) de incumplimiento asociadas con un aumento de 1 unidad en el ingreso (z.income), manteniendo constante el balance (z.balance). Un aumento de 1 unidad en el ingreso (z.income) es equivalente a un aumento de 2 desviaciones estándar en el ingreso (income) (2.667328^{4}). Este coeficiente también es estadísticamente significativo ya que .56 - 2 x .13 > 0.

¿Qué significa que las probabilidades logarítmicas de incumplimiento aumenten en 5.46 o .56? En términos precisos, ¿quién sabe? Para que estas cantidades tengan una mejor interpretación, necesitamos transformarlas, ya sea en probabilidades (odds) o en razones de probabilidades (razón de odds -> odds ratio). Sin embargo, debemos señalar que el signo y la magnitud de los coeficientes si son informativas: la relación con el incumplimiento del pago es positiva en ambos casos y, como ya se había visto de forma gráfica en los diagramas de cajas, el efecto del saldo (balance) es mucho mayor que el del ingreso (income).

6.3 Coeficientes de regresión logística como probabilidades

Podemos dar una interpretación más específica de la regresión logística más allá del efecto y magnitud. Para ello, podemos usar la función logit inversa para convertir las probabilidades logarítmicas (log-odds) de incumplimiento de pago en las tarjetas (cuando el saldo y los ingresos son promedio) en una probabilidad:

invlogit <- function(x) exp(x)/(1 + exp(x))
 invlogit(-6.13 + 5.46 * 0 + .56 * 0)

[1] 0.002171854

La probabilidad de incumplimiento para aquellos con un saldo promedio de tarjeta de crédito de (835.37) y un ingreso promedio de (3.351698^{4}) es de hecho bastante bajo: solo 0.002. Asimismo, podemos calcular el cambio en la probabilidad de incumplimiento en el pago asociado con un aumento de 1 unidad en el saldo, manteniendo el ingreso constante en el promedio (z.ingreso=0). Esto equivaldría a aumentar el saldo en casi 1000$, de 835.37 a 1802.8.

invlogit(-6.13 + 5.46 * 1) - invlogit(-6.13 + 5.46 * 0)

[1] 0.336325

6.4 Coeficientes de regresión logística como razones de odds

También podemos interpretar los coeficientes de regresión logística como razones de odds (OR).⁷ Si dos resultados tienen probabilidades \((p, 1-p)\), entonces \(\frac {p} {1-p}\) se conoce como odds (probabilidades o momio) del resultado. Las odds son simplemente diferentes formas de representar la misma información: \(\text{odds} = \frac{p}{1-p}\) y \(p = \frac{\text{odds}} {1+ \text{odds}}\). Por ejemplo, una odds de 1 es equivalente a una odds de .5 — es decir, resultados igualmente probables para \(p\) y \(1-p\): \(\text{odds(p = .5)} = \frac{.5}{1-.5} = 1\) y \(p(\text{oods} = 1) = \frac{\text{1}}{1 + 1} = .5.\)

La razón de dos odds es una OR:

\[ \frac{\frac{p_2}{1-p_2}}{\frac{p_1}{1-p_1}} \]

Una razón de odds se puede interpretar como un cambio en la probabilidad. Por ejemplo, un cambio en la probabilidad de \(p_1 = .33\) a \(p_2 = .5\) da como resultado un OR de 2, de la siguiente manera:

\[ \frac{\frac{.5}{.5}}{\frac{.33}{.66}} = \frac{1}{.5} = 2 \]

También podemos interpretar el OR como el aumento porcentual de las probabilidades de un evento. Aquí, un OR de 2 equivale a aumentar las probabilidades en un 100%, de 0,5 a 1.

Recordemos que representamos el modelo logit de esta manera:

\[ \text{log} \left(\frac{p}{1-p}\right) = \alpha + \beta x \]

La parte izquierda de la ecuación, expresado como logaritmos de probabilidades, está en la misma escala que la derecha derecho: \(\pm \infty\). Por lo tanto, no hay no linealidad en esta relación, y aumentar \(x\) en 1 unidad tiene el mismo efecto que en la regresión lineal: cambia el resultado predicho en \(\beta\). Entonces, pasar de \(x\) a \(x + 1\) equivale a sumar \(\beta\) a ambos lados de la ecuación anterior. Centrándonos solo en el lado izquierdo, tenemos, después de exponenciar, las probabilidades originales multiplicadas por \(e^\beta\):

\[ e^{\text{log} \left(\frac{p}{1-p}\right) + \beta} = \frac{p}{1-p} *e^ \beta \]

(ya que \(e^{a+b} = e^a*e^b\)).

Podemos pensar en \(e^\beta\) como el cambio de la odds del resultado cuando \(x\) cambia en 1 unidad, que se puede representar, utilizando la formulación anterior, como una OR:

\[ \frac{\frac{p_2}{1-p_2}}{\frac{p_1}{1-p_1}} = \frac{\frac{p_1}{1-p_1} * e^\beta }{\frac{p_1}{1-p_1}} = e^\beta. \]

Por lo tanto, \(e^\beta\) se puede interpretar como el cambio en las probabilidades asociadas con un aumento de 1 unidad en \(x\), expresado en términos porcentuales. En el caso de OR = \(\frac{1}{. 5} = 2\), el porcentaje de aumento en las probabilidades del resultado es del 100%. Cuando la OR es \(>2\) se suele expresar como x-veces más (OR=3.5, hay 3.5 veces más probabilidad de observar el evento que no observarlo cuando se cambia \(x\) en 1 unidad), y cuando la OR es \(<1\) se suele hablar de protección a no tener el evento y el porcentaje se calcula como 1-OR.

Apliquemos esta información a nuestro modelo anterior de aplicando la exponencial a los coeficientes de saldo e ingresos:

exp(5.46)

[1] 235.0974

exp(.56)

[1] 1.750673

Podemos interpretar estos ORs como el porcentaje o cambio multiplicativo en las probabilidades asociadas con un aumento de 1 unidad en el predictor (mientras se mantienen constantes los demás), de 1 a 235 (un aumento de 23,400%) en el caso de balance, y de 1 a 1,75 (un aumento del 75%) para los ingresos. Por ejemplo, podemos decir que la probabilidad de incumplimiento es un 75% mayor cuando los ingresos aumentan en 1 unidad.

Cuando las variables predictoras son categóricas (como en biomedicina: sexo, estadío tumoral, fumar, beber, …) la interpretación se hace más sencilla porque el cambio de 1 unidad en estas variables, supone el cambio de una categoría respecto a la basal (ya que se usan dummy variables). Así, por ejemplo, si nuestro outcomes tener cáncer de pulmón o no, y nuestro predictor es ser fumanor o no, si nuestros análisis nos dan una OR de 6 asociado a ser fumador, la interpretación sería: “La odds (probabilidad, abusando de lenguaje - también riesgo si el outcome es poco frecuente) de sufrir cáncer de pulmón es 6 veces mayor en los fumadores que en los no fumadores.

6.5 Capacidad predictiva de un modelo de clasificación

Podemos evaluar el rendimiento (es decir, la capacidad predictiva) del modelo logístico utilizando el AIC, así como mediante el uso de otras medidas como: la desviación (deviance) residual, la precisión, la sensibilidad, la especificidad y el área bajo la curva (AUC).

Al igual que AIC, la deviance es una medida de error, por lo que una deviance más baja significa un mejor ajuste a los datos. Esperamos que la desviación disminuya en 1 para cada predictor, por lo que con un predictor informativo (e.g. variable imporante para el modleo), la deviance disminuirá en más de 1. Deviance = \(-2ln(L)\), donde \(ln\) es el logaritmo natural y \(L\) es la función de verosimilitud . Veámoslo con nuestro ejemplo:

logistic_model1 <- 
  glm(default ~ balance, 
      data = Default, 
      family = binomial)
logistic_model1$deviance

[1] 1596.452

logistic_model2 <- 
  glm(default ~ balance + income, 
      data = Default, 
      family = binomial)
logistic_model2$deviance

[1] 1578.966

En este caso, la deviance se redujo en más de 1, lo que indica que los ingresos mejoran el ajuste del modelo. Podemos hacer una prueba formal de la diferencia usando, como para los modelos lineales, la prueba de razón de verosimilitud:

lrtest(logistic_model1, logistic_model2)

Model 1: default ~ balance
Model 2: default ~ balance + income

  L.R. Chisq         d.f.            P 
1.748541e+01 1.000000e+00 2.895205e-05 

También podemos traducir las probabilidades de un modelo logístico para el incumplimiento del pago en predicciones de clase asignando “Sí” (predeterminado) a probabilidades mayores o iguales a .5 y “No” (sin valor predeterminado) a probabilidades menores que .5, y luego contar el número de veces que el modelo asigna la clase predeterminada correcta. Si dividimos este número por el total de observaciones, habremos calculado la “precisión”. La precisión se utiliza a menudo como medida del rendimiento del modelo.

preds <- ifelse(fitted(logistic_model2) >= .5, "Yes", "No") 
(length(which(preds == Default$default)) / nrow(Default))*100

[1] 97.37

El modelo es 97.37% preciso. Valores superirores al 50% mostrarían una mejora en la predicción ya que por azar, se espera que el modelo tenga una precisión del 50%.

EJERCICIO (Entrega Moodle P2-Precisión)

Utiliza una simulación sencilla para demostrar que is asignamos por azar que una persona va a incumplir o no con los pagos, el valor esperado de la precisión de la variable “Default$default” es del 50%. NOTA: realiza 1000 simulaciones.

Una forma sencilla de obtener un buen modelo de clasificación sería asignar a todos la categoría más frecuente. Por ejemplo, en nuestros datos, la clase mayoritaria es “No” para la variable incimpliminto por un amplio margen (9667 a 333). La mayoría de las personas no incumplen. ¿Cuál es nuestra precisión, entonces, si simplemente predecimos “No” para todos los casos? La proporción de “No” en los datos es 96.67%, por lo que si siempre predijimos “No” esa sería nuestra precisión (9667 / (9667 + 333) = .9667). El modelo logístico, sorprendentemente, ofrece solo una ligera mejora.

Sin embargo, al evaluar el rendimiento del clasificador, debemos reconocer que no todos los errores son iguales y que la precisión tiene limitaciones como métrica de rendimiento. En algunos casos, el modelo puede haber predicho el incumplimiento cuando no lo había. Esto se conoce como “falso positivo”. En otros casos, el modelo puede haber predicho que no hubo incumplimiento cuando hubo incumplimiento. Esto se conoce como “falso negativo”. En la clasificación, utilizamos lo que se conoce como matriz de confusión para resumir estos diferentes tipos de errores del modelo, denominados así porque la matriz resume cómo se confunde el modelo. Usaremos la función confusionMatrix () de la librería caret para calcular rápidamente estos valores:

confusionMatrix(as.factor(preds), Default$default, positive = "Yes")

Confusion Matrix and Statistics

          Reference
Prediction   No  Yes
       No  9629  225
       Yes   38  108
                                          
               Accuracy : 0.9737          
                 95% CI : (0.9704, 0.9767)
    No Information Rate : 0.9667          
    P-Value [Acc > NIR] : 3.067e-05       
                                          
                  Kappa : 0.4396          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       
                                          
            Sensitivity : 0.3243          
            Specificity : 0.9961          
         Pos Pred Value : 0.7397          
         Neg Pred Value : 0.9772          
             Prevalence : 0.0333          
         Detection Rate : 0.0108          
   Detection Prevalence : 0.0146          
      Balanced Accuracy : 0.6602          
                                          
       'Positive' Class : Yes             
                                          

Esta función produce una gran cantidad de resultados. Podemos ver que informa la misma precisión que calculamos anteriormente: .97. La matriz de confusión 2 x 2 está en la parte superior. Podemos caracterizar estos 4 valores en la matriz de la siguiente manera:

1. 9629 negativos verdaderos (TN): cuando el modelo predice correctamente “No”
2. 108 verdaderos positivos (TP): cuando el modelo predice correctamente “Sí”
3. 225 falsos negativos (FN): cuando el modelo predice incorrectamente “No”
4. 38 falsos positivos (FP): cuando el modelo predice incorrectamente “Sí”

La siguiente tabla resume estas posibilidades:

	Reference
Predicted	No	Yes
No	TN	FN
Yes	FP	TP

Hay dos medidas clave, además de la precisión, para caracterizar el rendimiento del modelo. Mientras que la precisión mide el error general, la sensibilidad y la especificidad miden errores específicos de clase.

• Precisión = 1 - (FP + FN) / Total:

1 - (38 + 225) / 10000

[1] 0.9737

• Sensibilidad (o la tasa de verdaderos positivos): TP / (TP + FN). En este caso, la sensibilidad mide la proporción de incumplimientos que se clasificaron correctamente como tales.

108 / (108 + 225) 

[1] 0.3243243

• Especificidad (o la tasa de verdaderos negativos): TN / (TN + FP). En este caso, la especificidad mide la proporción de no incumplimientos que se clasificaron correctamente como tales.

9629 / (9629 + 38) 

[1] 0.9960691

¿Por qué deberíamos considerar estos errores específicos de clase? Todos los modelos tienen errores, pero no todos los errores del modelo son igualmente importantes. Por ejemplo, un falso negativo — prediciendo incorrectamente que un prestatario no incurrirá en incumplimiento — puede ser un error costoso para un banco, si el incumplimiento se hubiera podido prevenir mediante la intervención. Pero, por otro lado, un falso positivo, que predice incorrectamente que un prestatario incurrirá en incumplimiento, puede desencadenar una advertencia innecesaria que irrita a los clientes. Los errores que comete un modelo se pueden controlar ajustando el umbral de decisión utilizado para asignar probabilidades predichas a las clases. Usamos un umbral de probabilidad de .5 para clasificar los incumplimientos en los pagos. Si el umbral se establece en .1, por el contrario, la precisión general disminuiría, pero también lo haría el número de falsos negativos, lo que podría ser deseable. El modelo luego atraparía a más morosos, lo que ahorraría dinero al banco, pero eso tendría un costo: más falsos positivos (clientes potencialmente irritados).

preds <- as.factor(ifelse(fitted(logistic_model2) >= .1, "Yes", "No"))
confusionMatrix(preds, Default$default)$table

          Reference
Prediction   No  Yes
       No  9105   90
       Yes  562  243

La pregunta de cómo establecer el umbral de decisión —. 5, .1 o algo más — debe responderse con referencia al contexto empresarial.

Una curva de característica operativa del receptor (ROC por sus siglas en ingles) visualiza las compensaciones entre los tipos de errores trazando la especificidad frente a la sensibilidad. El cálculo del área bajo la curva ROC (AUC) nos permite, además, resumir el rendimiento del modelo y comparar modelos. La curva en sí muestra los tipos de errores que cometería el modelo en diferentes umbrales de decisión. Para crear una curva ROC usamos la función roc () del paquete pROC, y mostramos los valores de sensibilidad / especificidad asociados con los umbrales de decisión de .1 y .5:

library(pROC)
library(plotROC)
invisible(plot(roc(factor(ifelse(Default$default == "Yes", 1, 0)), 
                   fitted(logistic_model2)), print.thres = c(.1, .5), 
               col = "red", print.auc = T))

Un modelo con una precisión del 50%, es decir, un clasificador aleatorio, tendría una curva ROC que siguiera la línea de referencia diagonal. Un modelo con una precisión del 100%, un clasificador perfecto, tendría una curva ROC siguiendo los márgenes del triángulo superior. Cada punto de la curva ROC representa un par de sensibilidad / especificidad correspondiente a un umbral de decisión particular. Cuando establecimos el umbral de decisión en .1, la sensibilidad fue .73 (243 / (243 + 90)) y la especificidad fue .94 (9105 / (9105 + 562)). Ese punto se muestra en la curva. Del mismo modo, cuando establecimos el umbral de decisión en .5, la sensibilidad fue .32 y la especificidad fue .996. Ese punto también está en la curva. ¿Qué umbral de decisión es óptimo? Nuevamente, depende del problema (pensad en cáncer o en este ejemplo de dinero). Las curvas ROC nos permiten elegir los errores específicos de clase que podemos cometer.

El AUC es el resumen de cómo funciona el modelo en cada umbral de decisión. En general, los modelos con AUC más altos son mejores. Esta medida nos servirá para comparar métodos de aprendizaje automático que iremos aprendiendo durante el curso.

6.6 Ejemplo de regresión logística: modelización de riesgo diabetes

Para este ejemplo usaremos el conjunto de datos Pima, incluido en la librería MASS que contienen esta información:

Una población de mujeres que tenían al menos 21 años, de ascendencia indígena Pima y que vivían cerca de Phoenix, Arizona, se sometieron a pruebas de diabetes de acuerdo con los criterios de la Organización Mundial de la Salud. Los datos fueron recopilados por el Instituto Nacional de Diabetes y Enfermedades Digestivas y Renales de EE. UU. Usaremos información para 532 mujeres con datos completos después de eliminar los datos (principalmente faltantes) sobre la insulina sérica.

El conjunto de datos incluye las siguientes variables:

1. npreg: número de embarazos
2. glu: concentración de glucosa en plasma a 2 horas en una prueba de tolerancia oral a la glucosa
3. bp: presión arterial diastólica (mm Hg)
4. piel: grosor del pliegue cutáneo del tríceps (mm)
5. bmi: índice de masa corporal (peso en kg / (altura en m) ^ 2)
6. ped: función del pedigrí de la diabetes
7. age: Edad (años)
8. type: Sí o No (diabetes)

La variable resultado es “type”, que indica si una persona tiene diabetes. Los datos están divididos en un conjunto de entrenamiento y otro test que combinaremos para ilustrar este ejemplo.

library(MASS)
data("Pima.tr")
data("Pima.te")
d <- rbind(Pima.te, Pima.tr)
str(d)

'data.frame':   532 obs. of  8 variables:
 $ npreg: int  6 1 1 3 2 5 0 1 3 9 ...
 $ glu  : int  148 85 89 78 197 166 118 103 126 119 ...
 $ bp   : int  72 66 66 50 70 72 84 30 88 80 ...
 $ skin : int  35 29 23 32 45 19 47 38 41 35 ...
 $ bmi  : num  33.6 26.6 28.1 31 30.5 25.8 45.8 43.3 39.3 29 ...
 $ ped  : num  0.627 0.351 0.167 0.248 0.158 0.587 0.551 0.183 0.704 0.263 ...
 $ age  : int  50 31 21 26 53 51 31 33 27 29 ...
 $ type : Factor w/ 2 levels "No","Yes": 2 1 1 2 2 2 2 1 1 2 ...

Todos los predictores son enteros o numéricos. Nuestro objetivo es construir un modelo logístico de diabetes para ilustrar cómo interpretar los coeficientes del modelo.

6.6.1 Modelo simple

Comencemos con un modelo simple.

bin_model1 <- 
  glm(type ~ bmi + age, 
      data = d, 
      family = binomial)
display(bin_model1)

glm(formula = type ~ bmi + age, family = binomial, data = d)
            coef.est coef.se
(Intercept) -6.26     0.67  
bmi          0.10     0.02  
age          0.06     0.01  
---
  n = 532, k = 3
  residual deviance = 577.2, null deviance = 676.8 (difference = 99.6)

• Intercept: -6.26 es el logaritmo de la probabilidad de tener diabetes cuando bmi = 0 y edad = 0. Dado que ni la edad ni el bmi pueden ser iguales a 0, el intercept no es interpretable en este modelo. Por tanto, tendría sentido centrar las variables para facilitar la interpretación.
• bmi: .1 es el cambio previsto en el log-odds de diabetes asociado con un aumento de 1 unidad en el bmi, manteniendo la edad constante. Este coeficiente es estadísticamente significativo ya que .1 - 2 x .02 > 0. (Un IC del 95% que no contiene 0 indica significancia estadística) y también porque su p-valor asociacio mediante el test de score es \(<0.05\) (usar la función summary () . Podemos traducir este coeficiente en un OR mediante la exponencial: \(e^.1\) = 1.11. Un aumento de 1 unidad en el IMC, manteniendo la edad constante, se asocia con un aumento del 11% en la odds (o, más coloquialmente, la probabilidad) de diabetes.
• edad: .06 es el cambio predicho en el log-oods de diabetes asociado con un aumento de 1 unidad en la edad, manteniendo constante el bmi. Este coeficiente también es estadísticamente significativo ya que .06 - 2 x .01> 0. El OR para la edad es \(e^.06\) = 1.06 lo que indica un aumento del 6% en la probabilidad de sufrir diabetes asociada con un aumento de 1 unidad en la edad.

6.6.2 Agregar predictores y evaluar el ajuste

Ahora ajustaremos un modelo completo (excluyendo “skin”, ya que parece medir casi lo mismo que bmi). ¿Mejora el ajuste?

bin_model2 <- 
  glm(type ~ bmi + age + ped + glu + npreg + bp , 
      data = d, 
      family = binomial) 
display(bin_model2)

glm(formula = type ~ bmi + age + ped + glu + npreg + bp, family = binomial, 
    data = d)
            coef.est coef.se
(Intercept) -9.59     0.99  
bmi          0.09     0.02  
age          0.03     0.01  
ped          1.31     0.36  
glu          0.04     0.00  
npreg        0.12     0.04  
bp          -0.01     0.01  
---
  n = 532, k = 7
  residual deviance = 466.5, null deviance = 676.8 (difference = 210.3)

La deviance disminuye de 577 en el modelo anterior a 466.5 en este modelo, muy por encima de los 4 puntos que debería bajar simplemente al incluir 4 predictores adicionales. Además, el LRT nos indica que estas diferencias son estadísticamente significativas:

lrtest(bin_model1, bin_model2)

Model 1: type ~ bmi + age
Model 2: type ~ bmi + age + ped + glu + npreg + bp

L.R. Chisq       d.f.          P 
  110.6664     4.0000     0.0000 

El segundo modelo es mucho mejor que el primero, lo que también es evidente cuando observamos las matrices de confusión (con el umbral de decisión establecido en .5)

preds <- as.factor(ifelse(fitted(bin_model1) > .5, "Yes", "No"))
confusionMatrix(preds, d$type)$table

          Reference
Prediction  No Yes
       No  312 112
       Yes  43  65

preds <- as.factor(ifelse(fitted(bin_model2) > .5, "Yes", "No"))
confusionMatrix(preds, d$type)$table

          Reference
Prediction  No Yes
       No  318  73
       Yes  37 104

Como era de esperar, el modelo completo comete menos errores. Podemos formalizar esta impresión calculando y comparando la precisión del modelo: 1 - (112 + 43) / (112 + 43 + 312 + 65) = 0.71 para el primer modelo, en comparación con 1 - (73 + 37) / (73 + 37 + 318 + 104) =r round (1 - (73 + 37) / (73 + 37 + 318 + 104 ), 2) para el segundo modelo más grande. Los números de verdaderos negativos son cercanos, pero el modelo más grande aumenta sustancialmente el número de verdaderos positivos y reduce el número de falsos negativos, prediciendo incorrectamente que una persona no tiene diabetes (aunque esta sigue siendo la clase de error más grande). Deberíamos comprobar para ver que estos modelos mejoran la precisión sobre un modelo que siempre predice la clase mayoritaria. En este conjunto de datos, “No” es la categoría mayoritaria con 66,7%. Entonces, si siempre predijimos “No”, estaríamos en lo correcto el 66.7% de las veces, que es una precisión menor que cualquiera de los modelos. Las curvas ROC proporcionan una descripción más sistemática del rendimiento del modelo en términos de errores de clasificación errónea.

invisible(plot(roc(d$type,
                   fitted(bin_model1)),
               col = "red", 
               main = "ROC curves: logistic model 1 (red) vs. logistic model 2 (blue)"))
invisible(plot(roc(d$type,
                   fitted(bin_model2)),
               print.auc = T, 
               col = "blue", 
               add = T))

El modelo más grande es claramente mejor: el AUC es más alto.

6.6.3 Análisis de interacciones

Agreguemos una interacción entre dos predictores continuos, ped y glu. Centraremos y escalaremos las entradas para que el coeficiente de la interacción sea interpretable.

bin_model3 <- 
  standardize(glm(type ~ bmi  + ped * glu +  age + npreg + bp , 
      data = d, 
      family = binomial)) 
display(bin_model3)

glm(formula = type ~ z.bmi + z.ped * z.glu + z.age + z.npreg + 
    z.bp, family = binomial, data = d)
            coef.est coef.se
(Intercept) -1.02     0.13  
z.bmi        1.29     0.26  
z.ped        1.02     0.24  
z.glu        2.31     0.27  
z.age        0.53     0.30  
z.npreg      0.87     0.29  
z.bp        -0.19     0.26  
z.ped:z.glu -1.15     0.41  
---
  n = 532, k = 8
  residual deviance = 460.0, null deviance = 676.8 (difference = 216.7)

La interacción mejora el modelo pero no cambia la imagen general del ajuste del modelo en el gráfico de residuos agrupados (no incluido).

• Intercept : -1.02 es el log-odds de diabetes cuando todos los predictores son promedio (ya que hemos centrado y escalado las entradas). La probabilidad de tener diabetes para la mujer promedio en este conjunto de datos es logit\(^{- 1}\) (- 1.02) = 0.27.
• z.bmi: 1.29 es el log-odds de diabetes asociado con un aumento de 1 unidad en z.bmi, manteniendo constantes los demás predictores. \(e^{1.29}\) = 3.63 por lo que un aumento de 1 unidad en z.bmi, manteniendo constantes los otros predictores, se asocia con un aumento del 263% en la probabilidad de sufrir diabetes.
• z.ped: 1.02 es el log-odds de diabetes asociado con un aumento de 1 unidad en z.ped, cuando z.glu = 0 y manteniendo constantes los otros predictores. \(e^{1.02}\) = 2.77 por lo que un aumento de 1 unidad en z.ped, cuando z.glu = 0 y manteniendo los otros predictores constantes, se asocia con un aumento del 177% en la probabilidad de sufrir diabetes.
• z.glu: 2.31 es el log-odds de diabetes asociado con un aumento de 1 unidad en z.glu, cuando z.ped = 0 y manteniendo constantes los demás predictores. \(e^{2.31}\) = 10.07 por lo que un aumento de 1 unidad en z.glu, cuando z.ped = 0 y manteniendo los otros predictores constantes, se asocia con un aumento del 907% en la probabilidad de sufrir diabetes.
• …. el resto de predictores igual
• z.ped:z.glu : se añade -1.15 al log-odds de diabetes de z.ped cuando z.glu aumenta en 1 unidad, manteniendo constantes los otros predictores. O, alternativamente, se añade -1.15 al log-odds de diabetes de z.glu por cada unidad adicional de z.ped. Calculamos el OR, como en los otros casos, exponenciando: \(e^ {- 1.15}\) = 0.32. El OR para z.ped disminuye en un 68% (1 - .32 = .68) cuando z.glu aumenta en 1 unidad, manteniendo constantes los demás predictores. O, alternativamente, el OR para z.glu disminuye en un 68% con cada unidad adicional de z.ped.

6.6.4 Gráfico de la interacción

Como siempre, debemos visualizar la interacción para comprenderla. Esto es especialmente necesario cuando las relaciones se expresan en términos de probabilidades logarítmicas y razones de probabilidades. Como hemos hecho anteriormente para fines de visualización, dicotomizaremos los predictores en la interacción y, en este caso, para facilitar la interpretación, presentaremos las relaciones en términos de probabilidades. El propósito de los gráficos es la comprensión y la ilustración, por lo que no nos preocupa demasiado la precisión estadística. Resumiremos las relaciones usando una curva loess (estimación no paramétrica de la regresión) para capturar la no linealidad del efecto del predictor (\(\pm \infty\)) al rango del resultado binario (0, 1).

d$ped_bin <- ifelse(d$ped > mean(d$ped), "above avg","below avg")
d$glu_bin <- ifelse(d$glu > mean(d$glu), "above avg","below avg")
d$prob <- fitted(bin_model3)
d$type_num <- ifelse(d$type == "Yes", 1, 0)
ggplot(d, aes(glu, type_num)) +
  geom_point() +
  stat_smooth(aes(glu, prob, col = ped_bin), se = F) +
  labs(y = "Pr(diabetes)", title = "Diabetes ~ glu varying by ped_bin")

La relación entre glu y diabetes depende claramente de los niveles de ped. Como en el caso lineal, las líneas no paralelas indican una interacción. El coeficiente de log-odds negativo para la interacción del modelo indica que a medida que ped aumenta, la fuerza de la relación entre glu y type (diabetes) disminuye. Esto es exactamente lo que vemos en este gráfico:

ggplot(d, aes(ped, type_num)) +
  geom_point() +
  stat_smooth(aes(ped, prob, col = glu_bin), se = F) +
  labs(y = "Pr(diabetes)", title = "Diabetes ~ ped varying by glu_bin")

La interacción es más difícil de ver aquí porque la escala de ped está comprimida, con la mayoría de las observaciones cercanas a 0. Sin embargo, podemos ver que a medida que glu aumenta, la fuerza de la relación entre ped y diabetes disminuye. Nuevamente, las líneas no paralelas indican la presencia de una interacción.

6.6.5 Uso del modelo para predecir probabilidades

El tamaño del efecto más grande en el modelo anterior con la interacción, con mucho, es z.glu. Por tanto, para comunicar los resultados de este modelo deberíamos concentrarnos en glu. Pero los coeficientes expresados como logaritmos de probabilidades son algo confusos y, lamentablemente, las razones de probabilidades no ayudan a aclarar mucho las cosas. Deberíamos ir al trabajo adicional de traducir los coeficientes del modelo en probabilidades, pero para hacerlo debemos identificar en qué parte de la curva de probabilidad nos gustaría evaluar glu. Escogeremos la región cercana al promedio de z.glu — 0 — y examinaremos el efecto de aumentar z.glu en 1 unidad (que es igual a 2 desviaciones estándar de glu) cuando los otros predictores son promedio. La forma más sencilla de hacer esto es crear una base de datos con los valores de predicción deseados para usar con la función predict ().

basal <- data.frame(z.bmi = 0, z.ped = 0, z.glu = 0, z.age = 0, z.npreg = 0, z.bp = 0)
glucosa <- data.frame(z.bmi = 0, z.ped = 0, z.glu = 1, z.age = 0, z.npreg = 0, z.bp = 0)
(lo <- as.numeric(invlogit(predict(bin_model3, newdata = basal))))

[1] 0.2652028

(hi <- as.numeric(invlogit(predict(bin_model3, newdata = ))))

  [1] 0.716146591 0.025911688 0.017983270 0.029013744 0.893018807 0.665405710 0.395033985 0.211066662
  [9] 0.431695033 0.225375587 0.032840692 0.443676401 0.022607120 0.418795847 0.194085484 0.411388311
 [17] 0.122928759 0.771984610 0.760725696 0.088470567 0.953485181 0.785645203 0.028409690 0.026232726
 [25] 0.055457646 0.038181372 0.757850071 0.007630220 0.110507313 0.066603557 0.213559028 0.017048355
 [33] 0.294579371 0.213686248 0.245378175 0.147939308 0.058529536 0.047092679 0.169073996 0.023421517
 [41] 0.434375748 0.123923484 0.796680824 0.257623355 0.236877492 0.032937565 0.075796893 0.459769771
 [49] 0.014969537 0.214416331 0.409898562 0.024597682 0.736084480 0.071421984 0.863089942 0.109813203
 [57] 0.415148971 0.192888308 0.671705524 0.707181753 0.222687051 0.074378276 0.047288218 0.250211485
 [65] 0.102873418 0.451885105 0.016412055 0.839500999 0.593725195 0.887724350 0.023636511 0.972427462
 [73] 0.421875342 0.240818110 0.193220347 0.275433738 0.029285367 0.858427986 0.792694816 0.428056807
 [81] 0.241389504 0.359146427 0.159848949 0.126229936 0.089106804 0.965646716 0.021565093 0.504101659
 [89] 0.255700597 0.303575110 0.316087819 0.695623695 0.059157940 0.023850831 0.703188586 0.768676985
 [97] 0.198583192 0.736445668 0.019832290 0.911158959 0.794437117 0.023391174 0.082335119 0.536540158
[105] 0.491636172 0.875915114 0.830550408 0.453101254 0.050024142 0.017749687 0.151803493 0.175667620
[113] 0.900329614 0.125898408 0.410586667 0.586540636 0.832688763 0.035561910 0.252327286 0.110864806
[121] 0.054988088 0.131662542 0.393097130 0.492496368 0.513182966 0.032671230 0.036078371 0.218212971
[129] 0.686096496 0.281975363 0.145884637 0.351411775 0.239547073 0.725330877 0.392954359 0.354585338
[137] 0.328110094 0.183872054 0.649535926 0.189835213 0.750660386 0.100154335 0.253081898 0.306902095
[145] 0.138876568 0.112315917 0.602945913 0.033225614 0.021023766 0.126008933 0.074805807 0.593242488
[153] 0.276438109 0.030025810 0.388107478 0.814468072 0.819988723 0.320406177 0.512019175 0.189866084
[161] 0.019102564 0.778045038 0.063532324 0.021660839 0.087123658 0.371697775 0.017831740 0.174573799
[169] 0.072586081 0.050795966 0.321374939 0.390300727 0.242283766 0.101281785 0.477098938 0.177832066
[177] 0.851704388 0.421578801 0.059594070 0.801660505 0.263016503 0.158261796 0.629712769 0.680147379
[185] 0.627289909 0.143229332 0.371021820 0.048405335 0.155912833 0.844592024 0.774738072 0.074975116
[193] 0.094313006 0.034029316 0.156403796 0.765185930 0.162409946 0.941491860 0.098203334 0.096782855
[201] 0.013696456 0.086137707 0.665930779 0.328426313 0.230933621 0.045364772 0.067599667 0.128835144
[209] 0.552086895 0.192049128 0.391808734 0.634609576 0.138796509 0.034537822 0.402157076 0.473285370
[217] 0.915578006 0.109882027 0.060813059 0.078347861 0.511021110 0.046236210 0.637413042 0.057352160
[225] 0.060578365 0.328120324 0.117022441 0.150080125 0.033383541 0.548421342 0.727231653 0.273012605
[233] 0.006665104 0.023959522 0.212622505 0.263913260 0.304803089 0.282446484 0.528587177 0.053303093
[241] 0.042787259 0.895978588 0.951451693 0.260631793 0.057999466 0.077448786 0.039445750 0.078320973
[249] 0.521759602 0.120793953 0.109749728 0.115212481 0.063160685 0.315382778 0.186348273 0.242141093
[257] 0.800887574 0.866564275 0.147824356 0.483279794 0.458255040 0.058388003 0.052876258 0.244883311
[265] 0.761115562 0.021394588 0.531937078 0.751371443 0.307274782 0.801133302 0.005240345 0.656336098
[273] 0.109766806 0.101435208 0.023548802 0.056804667 0.081720211 0.513341731 0.014456472 0.127221754
[281] 0.744369844 0.498650317 0.182496315 0.605638441 0.021059011 0.094691510 0.633177077 0.189201907
[289] 0.229657632 0.818999683 0.081933947 0.793780642 0.798845854 0.137588419 0.275762056 0.057395933
[297] 0.445742740 0.584089416 0.070998635 0.209775098 0.174014948 0.294145785 0.743857124 0.126042513
[305] 0.893349297 0.712229694 0.217225821 0.135596074 0.343239883 0.187767754 0.233598292 0.822021476
[313] 0.082764167 0.102210424 0.113711375 0.122461316 0.828791325 0.107099346 0.054386091 0.931938668
[321] 0.326454054 0.678618259 0.863704140 0.217542759 0.056249096 0.816489720 0.491806299 0.084796292
[329] 0.943486215 0.302404878 0.132747361 0.037305147 0.053261653 0.902996173 0.042527808 0.787771463
[337] 0.031525411 0.220902831 0.054775507 0.642266770 0.520264613 0.707148463 0.902568015 0.695527581
[345] 0.848650675 0.466017086 0.291730904 0.027601848 0.133315806 0.771592655 0.383831136 0.090423918
[353] 0.020481869 0.115239515 0.038125985 0.017056834 0.070184592 0.755542001 0.050653582 0.481007350
[361] 0.195086963 0.187977925 0.169394579 0.377212811 0.151859655 0.475032680 0.186324471 0.612248541
[369] 0.023804839 0.007054577 0.085213924 0.171343641 0.900710787 0.150595032 0.257766676 0.013926514
[377] 0.035750449 0.185105131 0.253500436 0.369424730 0.640522481 0.841160836 0.060643162 0.051633327
[385] 0.731126069 0.154040714 0.063202801 0.337255163 0.087712142 0.059196122 0.446522479 0.839581017
[393] 0.899946965 0.051162921 0.569809764 0.085461942 0.080315826 0.343678259 0.504516436 0.092545675
[401] 0.175360179 0.310167781 0.720652085 0.094860697 0.905656358 0.432122047 0.742494468 0.704660222
[409] 0.065240504 0.044936864 0.086384200 0.460776518 0.058306508 0.538740077 0.360529087 0.413937339
[417] 0.248787673 0.063168325 0.205549345 0.143823024 0.121451790 0.204643255 0.028047843 0.262047997
[425] 0.797845942 0.020192814 0.067782381 0.880235920 0.067820578 0.093376849 0.145058897 0.885115279
[433] 0.023833728 0.172919691 0.023186441 0.343149266 0.491733925 0.334336420 0.182785009 0.107602309
[441] 0.159470782 0.480057756 0.718226410 0.110307283 0.282105476 0.504892387 0.023457816 0.246497165
[449] 0.288819782 0.406440777 0.333499471 0.789228958 0.024672867 0.158285199 0.915532293 0.083660243
[457] 0.616972917 0.156244026 0.066104438 0.117576505 0.582609228 0.935907656 0.494078320 0.906231057
[465] 0.044752423 0.098120327 0.817658808 0.063067848 0.048103513 0.094446457 0.009656172 0.457728456
[473] 0.371287692 0.732871561 0.239173153 0.048538537 0.427026826 0.838680152 0.005916331 0.639747396
[481] 0.192937289 0.044098532 0.307301900 0.771646053 0.893235289 0.633989608 0.325705459 0.861507987
[489] 0.924405484 0.211878609 0.232857992 0.553168378 0.199630964 0.167297598 0.636698113 0.102745817
[497] 0.412260537 0.037681077 0.499347951 0.561847838 0.108622049 0.036902563 0.275679034 0.142289294
[505] 0.929821372 0.853111736 0.052948764 0.186298288 0.034369934 0.468239290 0.021331580 0.617803456
[513] 0.083865824 0.751272138 0.115022996 0.496356896 0.458485256 0.303064449 0.512265744 0.935001789
[521] 0.011659745 0.243551235 0.234201771 0.733582682 0.439827183 0.132566837 0.200183137 0.254534757
[529] 0.620459940 0.187218540 0.128750723 0.801532793

round(hi - lo, 2)

  [1]  0.45 -0.24 -0.25 -0.24  0.63  0.40  0.13 -0.05  0.17 -0.04 -0.23  0.18 -0.24  0.15 -0.07  0.15
 [17] -0.14  0.51  0.50 -0.18  0.69  0.52 -0.24 -0.24 -0.21 -0.23  0.49 -0.26 -0.15 -0.20 -0.05 -0.25
 [33]  0.03 -0.05 -0.02 -0.12 -0.21 -0.22 -0.10 -0.24  0.17 -0.14  0.53 -0.01 -0.03 -0.23 -0.19  0.19
 [49] -0.25 -0.05  0.14 -0.24  0.47 -0.19  0.60 -0.16  0.15 -0.07  0.41  0.44 -0.04 -0.19 -0.22 -0.01
 [65] -0.16  0.19 -0.25  0.57  0.33  0.62 -0.24  0.71  0.16 -0.02 -0.07  0.01 -0.24  0.59  0.53  0.16
 [81] -0.02  0.09 -0.11 -0.14 -0.18  0.70 -0.24  0.24 -0.01  0.04  0.05  0.43 -0.21 -0.24  0.44  0.50
 [97] -0.07  0.47 -0.25  0.65  0.53 -0.24 -0.18  0.27  0.23  0.61  0.57  0.19 -0.22 -0.25 -0.11 -0.09
[113]  0.64 -0.14  0.15  0.32  0.57 -0.23 -0.01 -0.15 -0.21 -0.13  0.13  0.23  0.25 -0.23 -0.23 -0.05
[129]  0.42  0.02 -0.12  0.09 -0.03  0.46  0.13  0.09  0.06 -0.08  0.38 -0.08  0.49 -0.17 -0.01  0.04
[145] -0.13 -0.15  0.34 -0.23 -0.24 -0.14 -0.19  0.33  0.01 -0.24  0.12  0.55  0.55  0.06  0.25 -0.08
[161] -0.25  0.51 -0.20 -0.24 -0.18  0.11 -0.25 -0.09 -0.19 -0.21  0.06  0.13 -0.02 -0.16  0.21 -0.09
[177]  0.59  0.16 -0.21  0.54  0.00 -0.11  0.36  0.41  0.36 -0.12  0.11 -0.22 -0.11  0.58  0.51 -0.19
[193] -0.17 -0.23 -0.11  0.50 -0.10  0.68 -0.17 -0.17 -0.25 -0.18  0.40  0.06 -0.03 -0.22 -0.20 -0.14
[209]  0.29 -0.07  0.13  0.37 -0.13 -0.23  0.14  0.21  0.65 -0.16 -0.20 -0.19  0.25 -0.22  0.37 -0.21
[225] -0.20  0.06 -0.15 -0.12 -0.23  0.28  0.46  0.01 -0.26 -0.24 -0.05  0.00  0.04  0.02  0.26 -0.21
[241] -0.22  0.63  0.69  0.00 -0.21 -0.19 -0.23 -0.19  0.26 -0.14 -0.16 -0.15 -0.20  0.05 -0.08 -0.02
[257]  0.54  0.60 -0.12  0.22  0.19 -0.21 -0.21 -0.02  0.50 -0.24  0.27  0.49  0.04  0.54 -0.26  0.39
[273] -0.16 -0.16 -0.24 -0.21 -0.18  0.25 -0.25 -0.14  0.48  0.23 -0.08  0.34 -0.24 -0.17  0.37 -0.08
[289] -0.04  0.55 -0.18  0.53  0.53 -0.13  0.01 -0.21  0.18  0.32 -0.19 -0.06 -0.09  0.03  0.48 -0.14
[305]  0.63  0.45 -0.05 -0.13  0.08 -0.08 -0.03  0.56 -0.18 -0.16 -0.15 -0.14  0.56 -0.16 -0.21  0.67
[321]  0.06  0.41  0.60 -0.05 -0.21  0.55  0.23 -0.18  0.68  0.04 -0.13 -0.23 -0.21  0.64 -0.22  0.52
[337] -0.23 -0.04 -0.21  0.38  0.26  0.44  0.64  0.43  0.58  0.20  0.03 -0.24 -0.13  0.51  0.12 -0.17
[353] -0.24 -0.15 -0.23 -0.25 -0.20  0.49 -0.21  0.22 -0.07 -0.08 -0.10  0.11 -0.11  0.21 -0.08  0.35
[369] -0.24 -0.26 -0.18 -0.09  0.64 -0.11 -0.01 -0.25 -0.23 -0.08 -0.01  0.10  0.38  0.58 -0.20 -0.21
[385]  0.47 -0.11 -0.20  0.07 -0.18 -0.21  0.18  0.57  0.63 -0.21  0.30 -0.18 -0.18  0.08  0.24 -0.17
[401] -0.09  0.04  0.46 -0.17  0.64  0.17  0.48  0.44 -0.20 -0.22 -0.18  0.20 -0.21  0.27  0.10  0.15
[417] -0.02 -0.20 -0.06 -0.12 -0.14 -0.06 -0.24  0.00  0.53 -0.25 -0.20  0.62 -0.20 -0.17 -0.12  0.62
[433] -0.24 -0.09 -0.24  0.08  0.23  0.07 -0.08 -0.16 -0.11  0.21  0.45 -0.15  0.02  0.24 -0.24 -0.02
[449]  0.02  0.14  0.07  0.52 -0.24 -0.11  0.65 -0.18  0.35 -0.11 -0.20 -0.15  0.32  0.67  0.23  0.64
[465] -0.22 -0.17  0.55 -0.20 -0.22 -0.17 -0.26  0.19  0.11  0.47 -0.03 -0.22  0.16  0.57 -0.26  0.37
[481] -0.07 -0.22  0.04  0.51  0.63  0.37  0.06  0.60  0.66 -0.05 -0.03  0.29 -0.07 -0.10  0.37 -0.16
[497]  0.15 -0.23  0.23  0.30 -0.16 -0.23  0.01 -0.12  0.66  0.59 -0.21 -0.08 -0.23  0.20 -0.24  0.35
[513] -0.18  0.49 -0.15  0.23  0.19  0.04  0.25  0.67 -0.25 -0.02 -0.03  0.47  0.17 -0.13 -0.07 -0.01
[529]  0.36 -0.08 -0.14  0.54

6.7 Creación de un modelo y validación

Todo lo explicado en la sección de selección de variables para un modelo lineal, aplica para el caso de la regresión logística. Las funciones reconocen que el objeto es un glm con familia binomial y realiza los cálculos requeridos para este tipo de regresión (progamación orientada a objetos). El tema de validación cruzada para evaluar un modelo también aplica. Veámoslo con un ejemplo.

Supongamos que queremos elegir el mejor modelo para predecir el riesgo de diabetes y queremos validarlo con valización cruzada. Todos los pasos y métodos que hemos aprendido en las lecciones anteriores, podemos realizarlos de la siguiente manera. Usaremos los datos train y test que hay por defecto (Pima.tr y Pima.te). Para la validación cruzada usaremos la librería caret que veremos en detalle más adelante.

Empecemos seleccionando el mejor modelo en los datos de entrenamiento con un stepwise hacia atrás

mod.all <- glm(type ~ ., data=Pima.tr, family="binomial")
mod <- stats::step(mod.all, trace=FALSE, direction="backward")
summary(mod)

Call:
glm(formula = type ~ npreg + glu + bmi + ped + age, family = "binomial", 
    data = Pima.tr)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0009  -0.6816  -0.3664   0.6467   2.2898  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -9.938059   1.541571  -6.447 1.14e-10 ***
npreg        0.103142   0.064517   1.599  0.10989    
glu          0.031809   0.006667   4.771 1.83e-06 ***
bmi          0.079672   0.032649   2.440  0.01468 *  
ped          1.811417   0.661048   2.740  0.00614 ** 
age          0.039286   0.020967   1.874  0.06097 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 256.41  on 199  degrees of freedom
Residual deviance: 178.47  on 194  degrees of freedom
AIC: 190.47

Number of Fisher Scoring iterations: 5

Podemos evaluar la capacidad predictiva en la muestra test mediante

preds <- predict(mod, newdata = Pima.te, type="response")
preds <- as.factor(ifelse(preds >= .5, "Yes", "No") )
confusionMatrix(preds, Pima.te$type)

Confusion Matrix and Statistics

          Reference
Prediction  No Yes
       No  199  42
       Yes  24  67
                                          
               Accuracy : 0.8012          
                 95% CI : (0.7542, 0.8428)
    No Information Rate : 0.6717          
    P-Value [Acc > NIR] : 1.116e-07       
                                          
                  Kappa : 0.5294          
                                          
 Mcnemar's Test P-Value : 0.03639         
                                          
            Sensitivity : 0.8924          
            Specificity : 0.6147          
         Pos Pred Value : 0.8257          
         Neg Pred Value : 0.7363          
             Prevalence : 0.6717          
         Detection Rate : 0.5994          
   Detection Prevalence : 0.7259          
      Balanced Accuracy : 0.7535          
                                          
       'Positive' Class : No              
                                          

y calcular la capacidad predictiva en la muestra train utilizando un método de 5-fold cross-validation con:

library(caret)
mod.test <- train(type ~ ., data=Pima.tr, 
                  trControl = trainControl(method="cv", number=5),
                  method = "glm", family="binomial")
mod.test

Generalized Linear Model 

200 samples
  7 predictor
  2 classes: 'No', 'Yes' 

No pre-processing
Resampling: Cross-Validated (5 fold) 
Summary of sample sizes: 160, 160, 159, 160, 161 
Resampling results:

  Accuracy   Kappa    
  0.7553283  0.4381235

6.8 Nomogramas

Una vez que ya tenemos creado y validad un modelo predictivo, nos puede interesar aplicarlo en otros individuos para poder tomar decisiones. Para ello, podemos usar nomogramas.

Un nomograma es una representación gráfica que permite realizar con rapidez cálculos numéricos aproximados. Dentro del campo de la medicina, es frecuente que este tipo de gráficos este asociado al calculo de probabilidades de ocurrencia de un evento o una característica asociada a una enfermedad. Es lo que se conoce como Medicina Translacional.

Aunque existen otro tipo de herramientas de cálculo vía web para estas probabilidades (Shiny), el uso de nomogramas esta muy extendido en el campo de la biomedicina como por ejemplo en el calculo de probabilidades de recurrencia en distintos tipos de cáncer, la probabilidad de supervivencia a un mes tras un infarto, o el pronóstico tras un diagnóstico de cáncer a cierto tiempo (en ese caso se usan modelos de supervivencia. Así pues, la regresión logística será una de las herramientas con una aplicación más directa y sencilla en el aprendizaje automático, donde el uso de los modelos predictivos suelen tener un aplicabilidad directa en la población.

Existen numerosas herramientas para crear nomogramas en R, empecemos con la creación de nomogramas sencillos. Para ello usaremos los datos del ejemplo de diabetes con el modelo que hemos obtenido y validado anteriormente. Para usar la librería rms necesitamos que el modelo esté estimado con la función lrm ()

library(rms)

t.data <- datadist(Pima.tr)
options(datadist = 't.data')

mod.lrm <- lrm(type ~ npreg + glu + bmi + ped + age, data=Pima.tr)

nom <- nomogram(mod.lrm, fun = plogis, funlabel="Risk of diabetes")
plot(nom)

Supongamos que llega a la consulta una persona con un bmi de 35. Eso sumaría 30 puntos (basta con subir hacia arriba y ver qué valor de ‘Points’ corresponde a un bmi de 35). Una edad de 50 años (~22 puntos), una función del pedigrí de la diabetes de 1.8 (~62 puntos), una glucosa de 120 (~50 puntos) y 0 embarazos que sumaría 0 puntos. En total, el paciente suma un total de 164 puntos. Si ahora vamos a la línea de ‘Total Points’ y proyectamos sobre el predictor lineal de aproximadamente ~1.9 que se asocia con un riesgo de diabetes ligeramente superior al 80% (proyectar sobre ‘Risk of diabetes’).

Obviamente estos cálculos se pueden hacer de forma más directa calculando la predicción sobre este individuo con el objeto de R

indiv <- data.frame(bmi=35, age=50, ped=1.8, glu=120, npreg=0)
predict(mod.lrm, newdata = indiv, type="lp")

       1 
1.892376 

predict(mod.lrm, newdata = indiv, type = "fitted")

        1 
0.8690262 

Estos cálculos se pueden programar en R y hace una función, o también una aplicación Shiny para aquellos médicos que no sepan usar R (los nomogramas se siguen utilizando). Otra opción es que hagamos uso de una librería para hacer nomogramas dinámicos (crea Shiny app) de forma sencilla con la librería DynNom

library(DynNom)
DynNom(mod, Pima.tr)

Con estas simples instrucciones obtendríamos esta aplicación Shiny (Figura abajo) donde cada intervalo de confianza corresponde a un cálculo obtenido variando alguna de las variables predictoras

Nomograma para el modelo de predicción para diabetes

EJERCICIO (Entrega en Moodle: P-Práctica regresión logística):

En esta página https://vincentarelbundock.github.io/Rdatasets/datasets.html tienes acceso a cientos de bases de datos que hay en R y que también puedes descargar como ficheros .csv.

Elige una base de datos donde la pregunta científica requiera el análisis de una variable respuesta binaria. Selecciona unas 8-15 variables independientes (intenta que hayan categóricas y continuas). Crea un pdf en el que muestres un análisis completo de los datos incluyendo:

1. Descripción de tu pregunta científica

2. Análisis descriptivo de las variables independientes en función de tu variable respuesta

3. Imputación de datos en caso de ser necesario

4. Creación de un modelo predictivo

5. Descripción de las ORs para el modelo final

6. Validación del modelo mediante CV

7. Descripción de la capacidad predictiva del modelo (Sensibilidad, Precisión, curva ROC, …)

8. Creación de un nomograma estático que incluya la ilustración del cálculo de riesgo para un individuo ficticio.

9. Apéndice con el código R utilizado

Se trata de crear un documento de unas 2-3 páginas (figuras y tablas aparte) explicando los principales resultados del estudio donde cada párrafo podría corresponder a cada una de las tareas anteriores. El documento debe incluir las tablas y/o figuras que creas conveniente que sostengan los resultados descritos en el documento. Se trata de intentar escribir el apartado de resultados de un artículo científico. Antes de escribir este apartado de resultados, el documento deberá empezar con un apartado breve de métodos describiendo los datos (brevemente - podéis usar la información que hay en la página de donde habéis descargado los datos) y los métodos estadísticos utilizados. Aquellos que quieran, lo pueden hacer con R Markdown. NOTA: El que no suba un pdf tendrá un 0 en la práctica.

---

6. GLM significa modelo lineal generalizado. Esta función se ajustará a una variedad de modelos no lineales. Por ejemplo, podríamos especificar una regresión de Poisson con family = poisson.↩︎

7. Para obtener más detalles, se puede consultar este texto [How to interpret odds ratio in logistic regression?] (Https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-interpret-odds-ratios-in-logistic-regression/).↩︎