1. Modelos generalizados

En la unidad anterior revisamos los conceptos del modelo de regresión lineal, donde la modelación se realiza bajo el supuesto de la distribución normal donde la media es una combinación lineal de variables explicativas. Ahora, en salud muchas veces trabajamos con variables categóricas donde la normalidad no aplica y se deben usar otro tipo de distribuciones para la respuesta condicionada a los predictores. Donde vamos a reemplazar la media condicionar de \(Y\) como una combinación lineal de los predictores por una función de esa combinación lineal.

Los modelos generalizados se conforman de las siguientes tres partes:

  • La distribución de la varible respuesta condicionada a los predictores, que debe pertenecer a la familia exponencial.
  • Los predictores deben ser una combinación lineal.
  • Debemos usar una función de enlace \(g()\), que define como la combinación lineal de predictores \((\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_{p - 1} x_{p - 1})\) está nrelacionada con la respuesta condicionada a los predictores, \(\text{E}[Y \mid {\bf X} = {\bf x}]\).

Dependiendo del tipo de variable respuesta que nos interese modelar se utilizan ciertas funciones de enlace, en la siguiente tabla se muestran unas de las más comunes:

Regresión lineal Regresión Poisson Regresión logística
\(Y \mid {\bf X} = {\bf x}\) \(N(\mu({\bf x}), \sigma^2)\) \(\text{Pois}(\lambda({\bf x}))\) \(\text{Bern}(p({\bf x}))\)
Distribución Normal Poisson Bernoulli (Binomial)
\(\text{E}[Y \mid {\bf X} = {\bf x}]\) \(\mu({\bf x})\) \(\lambda({\bf x})\) \(p({\bf x})\)
Soporte Real: \((-\infty, \infty)\) Entero: \(0, 1, 2, \ldots\) Entero: \(0, 1\)
Uso Dstos numéricos Conteos Binario
Nombre enlace Identidad Log Logit
Función de enlace \(\eta({\bf x}) = \mu({\bf x})\) \(\eta({\bf x}) = \log(\lambda({\bf x}))\) \(\eta({\bf x}) = \log \left(\frac{p({\bf x})}{1 - p({\bf x})} \right)\)
media de la función \(\mu({\bf x}) = \eta({\bf x})\) \(\lambda({\bf x}) = e^{\eta({\bf x})}\) \(p({\bf x}) = \frac{e^{\eta({\bf x})}}{1 + e^{\eta({\bf x})}} = \frac{1}{1 + e^{-\eta({\bf x})}}\)

Fuente: Weisberg (2005) – Applied Linear Regression

En los modelos generalizados se relaciona el promedio de la respuesta con una combinación lineal de los predictores a través de una función de enlace. De esta forma,

\[ \eta({\bf x}) = g\left(\text{E}[Y \mid {\bf X} = {\bf x}]\right). \] donde el promedio es:

\[ \text{E}[Y \mid {\bf X} = {\bf x}] = g^{-1}(\eta({\bf x})). \]

2. Modelo de Regresión Logístico

2.1 Definición

Cuando tenemos una variable binaria (bernoulli) o dicotómica, donde la respuesta es si o no (1 o 0), usamos el modelo logístico. Este se define de la siguiente forma:

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = \beta_0 + \beta_1 x_1 + \ldots + \beta_{p - 1} x_{p - 1} \] Notemos que la parte derecha de la ecuación es la misma que trabajamos en regresión lineal, el lado izquierdo es el logaritmo de los odds (oportunidad), recordemos que el odds es la razón de dos probabilidades, donde en el numerador tenemos la probabilidad de que ocurra el evento sobre su complemento (no ocurra el evento). Luego el logaritmo del odds es la tranformación logit aplicada a \(p({\bf x})\)

\[ \text{logit}(\xi) = \log\left(\frac{\xi}{1 - \xi}\right) \] Aqui debemos definir la función inversa del logit, conocida como la función “logística” o sigmoide, donde los valores están entre 0 y 1.

\[ \text{logit}^{-1}(\xi) = \frac{e^\xi}{1 + e^{\xi}} = \frac{1}{1 + e^{-\xi}} \]

Ahora el modelo logístico queda de la forma:

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = \beta_0 + \beta_1 x_1 + \ldots + \beta_{q} x_{q} \]

Aplicando la transformación logit inversa podemos obtener una expresión para \(p({\bf x})\)

\[ p({\bf x}) = P[Y = 1 \mid {\bf X} = {\bf x}] = \frac{e^{\beta_0 + \beta_1 x_{1} + \cdots + \beta_{p-1} x_{(p-1)}}}{1 + e^{\beta_0 + \beta_1 x_{1} + \cdots + \beta_{p-1} x_{(p-1)}}} \]

2.2. Ajuste del modelo logístico

Con n observaciones, la ecuación del modelo logístico indexado en \(i\) es de la siguiente forma:

\[ \log\left(\frac{p({\bf x_i})}{1 - p({\bf x_i})}\right) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)} \] Para obtener la probabilidad de cada observación, se puede aplicar el inverso de la transformación logit:

\[ p({\bf x_i}) = P[Y_i = 1 \mid {\bf X_i} = {\bf x_i}] = \frac{e^{\beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)}}}{1 + e^{\beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)}}} \] \[ 1 - p({\bf x_i}) = P[Y_i = 0 \mid {\bf X} = {\bf x_i}] = \frac{1}{1 + e^{\beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)}}} \] Para ajustar este modelo, es decir, obtener las estimaciones de los diferentes beta´s, se utiliza el método de máxima verosimilitud. El cual busca el valor para beta que maximice la probabilidad de los datos observados.

\[ \boldsymbol{{\beta}} = [\beta_0, \beta_1, \beta_2, \beta_3, \ldots, \beta_{p - 1}] \]

Primero, se escribe la verosimilitud dado los datos observados:

\[ L(\boldsymbol{{\beta}}) = \prod_{i = 1}^{n} P[Y_i = y_i \mid {\bf X_i} = {\bf x_i}] \] Al reemplazar las probabilidades descritas arriba en el modelo, que son una función de los parámetros, se obtiene lo siguiente:

\[ L(\boldsymbol{{\beta}}) = \prod_{i = 1}^{n} p({\bf x_i})^{y_i} (1 - p({\bf x_i}))^{(1 - y_i)} \] \[ L(\boldsymbol{{\beta}}) = \prod_{i : y_i = 1}^{n} p({\bf x_i}) \prod_{j : y_j = 0}^{n} (1 - p({\bf x_j})) \] \[ L(\boldsymbol{{\beta}}) = \prod_{i : y_i = 1}^{} \frac{e^{\beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)}}}{1 + e^{\beta_0 + \beta_1 x_{i1} + \cdots + \beta_{p-1} x_{i(p-1)}}} \prod_{j : y_j = 0}^{} \frac{1}{1 + e^{\beta_0 + \beta_1 x_{j1} + \cdots + \beta_{p-1} x_{j(p-1)}}} \]

Contrario al modelo de regresión lineal, para este proceso de maximización no hay una solución analítica, por tanto, requiere una aproximación numérica, que R hace usando un algoritmo de iteración.

Pruebas sobre los parámetros

Como en la regresión lineal, en la regresión logística se pueden probar hipótesis para uno o múltiples parámetros:

Test de Wald

Las hipótesis a plantear en el modelo son:

\[ H_0: \beta_j = 0 \quad \text{vs} \quad H_1: \beta_j \neq 0 \]

Recordemos que en este caso tenemos el modelo:

\[ \log\left(\frac{p({\bf x})}{1 - p({\bf x})}\right) = \beta_0 + \beta_1 x_1 + \ldots + \beta_{p - 1} x_{p - 1} \] Con las hipótesis planteadas anteriormente, el estadístico de prueba toma la forma de una prueba z, asumiendo aproximación a la normal, teniendo una muestra lo suficientemente grande.

\[ z = \frac{\hat{\beta}_j - \beta_j}{\text{SE}[\hat{\beta}_j]} \overset{\text{approx}}{\sim} N(0, 1) \]

Test de razón de verosimilitud (LR test)

Cuando queremos comparar dos modelos, existen diferentes pruebas o estadísticos que ayudan a tomar una decisión, sóbre cuál sería el mejor. El LR test, se aplica cuando comparamos modelos anidados, es decir, cuando un modelo contiene el otro, o simplemente cuando queremos saber si al agregar más variables en un modelo, este tendría un mejor comportamiento. A continuación se explica cómo funciona esta prueba:

Considerando el siguiente modelo full o saturado con \(p\) parámetros. Se denotará la función de verosimilitud (MLE) de estos \(\beta\)-parametros como \(\hat{\beta}_{\text{Full}}\).

\[ \log\left(\frac{p({\bf x_i})}{1 - p({\bf x_i})}\right) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_{(p-1)} x_{i(p-1)} \] Ahora considerando el modelo nulo o reducido

\[ \log\left(\frac{p({\bf x_i})}{1 - p({\bf x_i})}\right) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}+ \cdots + \beta_{(q-1)} x_{i(q-1)} \]

donde \(q < p\). Este modelo tiene \(q - 1\) predictores, para un total de \(q\) \(\beta\)-parametros. Se denotará la función MLE de estos \(\beta\)-parametros como \(\hat{\beta}_{\text{Null}}\).

La diferencia entre estos dos modelos, se puede ver a través del siguiente planteamiento de hipótesis:

\[ H_0: \beta_q = \beta_{q+1} = \cdots = \beta_{p - 1} = 0. \] Lo anterior implica que, el modelo reducido está anidado en el modelo full.

Se define el siguiente estadístico \(D\), para la prueba de hipótesis anterior.

\[ D = -2 \log \left( \frac{L(\boldsymbol{\hat{\beta}_{\text{Null}}})} {L(\boldsymbol{\hat{\beta}_{\text{Full}}})} \right) = 2 \log \left( \frac{L(\boldsymbol{\hat{\beta}_{\text{Full}}})} {L(\boldsymbol{\hat{\beta}_{\text{Null}}})} \right) = 2 \left( \ell(\hat{\beta}_{\text{Full}}) - \ell(\hat{\beta}_{\text{Null}})\right) \]

donde \(L\) denota la verosimilitud y \(\ell\) el logaritmo de la verosimilitud. Para una muestra grande, este estadístico tiene una distribución aproximada chi2.

\[ D \overset{\text{approx}}{\sim} \chi^2_{k} \] donde \(k = p - q\), es la diferencia en el número de parámetros entre los dos modelos.

Este test es denominado Razón de verosimilitud o en inglés Likelihood-Ratio Test, sería el análogo al test \(F\) que se tiene en regresión lineal múltiple.

2.3. Evaluación del Modelo logístico

Para determinar el ajuste del modelo logístico se puede usar el test de Hosmer-Lemeshow, el cual plantea como hipótesis nula que el modelo se ajusta bien a los datos, lo que se busca es no recharzarla.

Para la evaluación del modelo logístico, con un objetivo predictivo, se pueden usar varias métricas cómo son la exactitud (accuracy), precisión, tasa de verdaderos positivos (recall), F1 score y el area bajo la curva ROC. Los cuales se calculan de la siguiente forma:

  • Accuracy: proporción de los correctamente clasificados.

\[ Accuracy= \frac{VerdaderosPositivos+VerdaderosNegativos}{Total} \]

  • Precision: se enfoca en la exactitud de las predicciones positivas

\[ Precision= \frac{VerdaderosPositivos}{VerdaderosPositivos+FalsosPositivos} \]

  • Tasa de verdaderos positivos (Recall): mide la proporción de los casos positivos correctamente predichos entre todos los posibles positivos.

\[ Recall= \frac{VerdaderosPositivos}{VerdaderosPositivos+FalsosNegativos} \]

  • F1 score: es la media armónica de la precisión y el recall.

\[ F1Score= 2*\frac{Precision*Recall}{Precision+Recall} \]

Fuente: https://www.geeksforgeeks.org/auc-roc-curve/

  • Verdaderos positivos: casos positivos predichos como positivos.
  • Verdaderos negativos: casos negativos predichos como negativos.
  • Falsos positivos: casos negativos predichos como positivos.
  • Falsos negativos: casos positivos predichos como negativos.

2.2.1. Area bajo la curva ROC

El area bajo la curva (AUC) característica de operaciones (ROC) representa el desempeño de la clasificación binaria del modelo. En el eje X muestra la proporción de falsos positivos mientras en el eje Y, la proporción de verdaderos positivos identificados (sensibilidad), ver matriz de confusión. El área está entre cero y uno, entre más cercana a uno sea, el rendimiento es mejor. El objetivo es maximizar el area, en orden de tener la sensibilidad más alta y el porcentaje de falsos positivos más bajo. EL AUC mide la probabilidad de que el modelo asigne a un caso positivo elegido al azar una probabilidad predicha más alta en comparación con un caso negativo elegido al azar.

El AUC Representa la probabilidad del modelo para distinguir entre las dos clases (variable respuesta).

Fuente: https://www.youtube.com/watch?app=desktop&v=vrMu2XwxTok

La relación entre la proporción de verdaderos positivos y falsos positivos es directa, un aumento en los verdaderos positivos lleva a un aumento en los falsos positivos.

2.4. Aplicación en R

Para el modelo de regresión logístico donde la respuesta es una variable dicotómica se usará una base de datos de la librería mlbench con las variables: diabetes (positivo, negativo), glucose (concentración de glucosa mg/dL), pressure (presión diastólica mm Hg), mass (IMC) y age (edad).

Para cargar los datos:

#install.packages("mlbench")
library("mlbench")
data("PimaIndiansDiabetes2", package = "mlbench")

Para renombrar la base de datos:

data_diabetes<-PimaIndiansDiabetes2
names(data_diabetes)
## [1] "pregnant" "glucose"  "pressure" "triceps"  "insulin"  "mass"     "pedigree"
## [8] "age"      "diabetes"
#omitimos los datos perdidos
data_diabetes = na.omit(data_diabetes)

La estimación del modelo logístico teniendo como variable respuesta la diabetes (Si/No) se realiza a través de la función \(glm\), teniendo en cuenta la familia de la distribución binomial.

modelo2 <- glm( diabetes ~ glucose + pressure + mass + age, data = data_diabetes, family = binomial)
summary(modelo2)
## 
## Call:
## glm(formula = diabetes ~ glucose + pressure + mass + age, family = binomial, 
##     data = data_diabetes)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.582963   1.136690  -8.431  < 2e-16 ***
## glucose      0.036352   0.004931   7.373 1.67e-13 ***
## pressure    -0.002339   0.011422  -0.205 0.837752    
## mass         0.078953   0.020789   3.798 0.000146 ***
## age          0.054867   0.013810   3.973 7.10e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 498.10  on 391  degrees of freedom
## Residual deviance: 354.32  on 387  degrees of freedom
## AIC: 364.32
## 
## Number of Fisher Scoring iterations: 5

Teniendo en cuenta la definición del modelo, recordemos que la función logit aplica logaritmo al odds, por tanto, las estimaciones de los beta´s se encuentran en logaritmo, para calcular los OR´s debemos exponenciar de la siguiente forma.

exp(coef(modelo2)) # exponenciar los coeficientes
##  (Intercept)      glucose     pressure         mass          age 
## 6.889251e-05 1.037021e+00 9.976638e-01 1.082153e+00 1.056400e+00
exp(confint(modelo2)) # exponenciar los límites del intervalo de confianza
## Waiting for profiling to be done...
##                    2.5 %       97.5 %
## (Intercept) 6.621589e-06 0.0005775797
## glucose     1.027395e+00 1.0475076526
## pressure    9.755741e-01 1.0204584438
## mass        1.039763e+00 1.1283792669
## age         1.028770e+00 1.0861971947

Ahora vamos a interpretar los OR´s:

  • Para la glucosa, por cada unidad de aumento en mg/dL la oportunidad de presentar diabetes aumenta 3.7%, independiente del efecto de la presión diastólica, IMC y edad.
  • Para la presión diastólica, por cada unidad de aumento en mmHg, la oportunidad de presentar diabetes disminuye en 0.3%, independiente del efecto de la presión diastólica, IMC y edad.
  • Por cada unidad de aumento en el IMC, la oportunidad de presentar diabetes aumenta 8.2%, independiente del efecto de la presión diastólica, IMC y edad.
  • Por cada año de aumento en la edad, la oportunidad de presentar diabetes aumenta 5.6%, independiente del efecto de la presión diastólica, IMC y edad.

A través de los intervalos de confianza podemos ver que en el caso de la presión diastólica, este va de 0.975 a 1.020, lo cual muestra que toma la unidad y valores muy cercanos, si nos devolvemos a la tabla de estimación podemos observa que es la única variable que no resulta significativa dado que su valor p es grande (mayor a 0.05, usual)

Para el análisis de los residuales de devianza, se puede obtener la tabla con la siguiente función:

anova(modelo2, test="Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: diabetes
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       391     498.10              
## glucose   1  111.432       390     386.67 < 2.2e-16 ***
## pressure  1    3.855       389     382.81 0.0496090 *  
## mass      1   11.595       388     371.22 0.0006612 ***
## age       1   16.893       387     354.32 3.956e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En la tabla anterior se observa como la devianza de los residuales del modelo nulo (sin covariables) es mayor respecto a los modelos con las covariables que se agregaron al modelo, lo cual es lo que se espera, que al agregar covariables se aumente la explicación de la oportunidad de presenatr diabetes.

Para analizar el test de **test de Hosmer-Lemeshow*, se utlizará la función creada en la librería ResourceSelection, de la siguiente forma:

#install.packages("ResourceSelection")
library(ResourceSelection)

hoslem.test(data_diabetes$diabetes, fitted(modelo2), g=10)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  data_diabetes$diabetes, fitted(modelo2)
## X-squared = 392, df = 8, p-value < 2.2e-16

El resultado anterior muestra que, se rechaza la hipótesis nula, por tanto el modelo no se ajusta bien a los datos.

Otra forma de evaluar la bondad de ajuste del modelo se puede hacer uso del pseudo \(R^2\) de McFadden.

#install.packages("pscl")
library(pscl)
pR2(modelo2)
## fitting null model for pseudo-r2
##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -177.1618509 -249.0489027  143.7741034    0.2886463    0.3070315    0.4268161

Este muestra que este modelo explica el 28.8% de variabilidad de la oportunidad de diabetes, el cual es bajo.

Evaluación del modelo

Para generar la matriz de confusión se utilizan las siguientes funciones:

#se estima de nuevo el modelo excluyendo los datos perdidos
modelo2 <- glm( diabetes ~ glucose + pressure + mass + age, data = data_diabetes, family = binomial)
data_diabetes$diabetes <- ifelse(data_diabetes$diabetes == "neg", 0, 1)
predictions <- predict(modelo2, type = "response")

#para convertir las probabilidades a predicciones binarias
predicted_clases <- ifelse(predictions > 0.5, 1, 0)

# matrix de confusion
conf_matrix <- table(data_diabetes$diabetes, predicted_clases)
library(caret)
confusionMatrix(conf_matrix)
## Confusion Matrix and Statistics
## 
##    predicted_clases
##       0   1
##   0 234  28
##   1  52  78
##                                           
##                Accuracy : 0.7959          
##                  95% CI : (0.7526, 0.8347)
##     No Information Rate : 0.7296          
##     P-Value [Acc > NIR] : 0.001468        
##                                           
##                   Kappa : 0.5172          
##                                           
##  Mcnemar's Test P-Value : 0.010127        
##                                           
##             Sensitivity : 0.8182          
##             Specificity : 0.7358          
##          Pos Pred Value : 0.8931          
##          Neg Pred Value : 0.6000          
##              Prevalence : 0.7296          
##          Detection Rate : 0.5969          
##    Detection Prevalence : 0.6684          
##       Balanced Accuracy : 0.7770          
##                                           
##        'Positive' Class : 0               
##

El modelo clasifica de manera correcta personas con diabetes y no, con una excatitud del 79.6% IC95%(75.3%; 83.5%). La sensibilidad es de 81.8%, es decir que, la capacidad del modelo de identificar correctamente los casos de diabetes es de 82 casos aproximandamente por cada 100. Mientras que, la especificidad indica que el modelo descarta diabetes en 74 casos aproximadamente por cada 100.

  • Accuracy: la exactitud del modelo en la predicción correcta de diabetes basado en las variables del modelo es del 79.6%.
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(paste("Accuracy:", accuracy))
## [1] "Accuracy: 0.795918367346939"
  • Precisión: en cuanto a la proporción de predicciones positivas correctas en relación con el total de predicciones positivas (VP+FP), el modelo identifica el 73.6%.
precision <- conf_matrix[2,2] / sum(predicted_clases)
print(paste("Precision:", precision))
## [1] "Precision: 0.735849056603774"
  • Recall: capacidad del modelo de identificar correctamente los casos positivos (VP) respecto al total de casos positivos reales (VP+FN), que en este caso corresponde a 60%.
recall <- conf_matrix[2,2] / sum(data_diabetes$diabetes)
print(paste("Recall:", recall))
## [1] "Recall: 0.6"
  • F1 Score: corresponde con la media armónica entre la precisión y el Recall, entre más cercano a uno sea, es mejor. Para este modelo fue del 66.1%.
f1_score <- 2 * (precision * recall) / (precision + recall)
print(paste("F1 Score:", f1_score))
## [1] "F1 Score: 0.661016949152542"
  • Curva ROC y AUC: para gráficar la curva ROC se utiliza la librería pROC.
#install.packages("pROC")
library(pROC)
roc_curve <- roc(data_diabetes$diabetes, predictions)
plot(roc_curve, main = "Curve ROC")

auc<- auc(roc_curve)
print(auc)
## Area under the curve: 0.8488

El area bajo la curva corresponde a 84.8% indicando un buen ajuste del modelo.

Finalmente, si quisieramos comparar este modelo frente a otro que no contenga una o varias variables, como la edad, podemos hacerlo a través del LR test, de la siguiente forma:

modelo_completo <- glm( diabetes ~ glucose + pressure + mass + age, data = data_diabetes, family = binomial)

modelo_reducido <- glm( diabetes ~ glucose + pressure + mass , data = data_diabetes, family = binomial)

# Likelihood Ratio Test (LR test)
anova(modelo_reducido, modelo_completo, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: diabetes ~ glucose + pressure + mass
## Model 2: diabetes ~ glucose + pressure + mass + age
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       388     371.22                          
## 2       387     354.32  1   16.892 3.956e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Lo cual nos dice que, el agregar la edad (variable adicional en el modelo completo), mejora significativamente el ajuste (valor p=0.0000395).

Evaluación de un término de interacción en el modelo

En los estudios epidemiológicos siempre es importante evaluar las variables que puedan sesgar, confundir o modificar los efectos medidos. Para la confusión, el ajuste de variables que se realiza en los modelos múltiples, es parte del proceso que permite evaluar su efecto ajustado. Pero, para la modificiación del efecto, es necesario incluir términos de interacción y medir su significancia para establecer si existe o no este efecto.

A continuación se simulará un ejemplo para mostrar, cómo se ingresa la interacción a un modelo de regresión logística, además cómo se interpreta.

Primero se simulará una base de datos de 200 observaciones, con tres variables: enfermedad cardiovascular, fumar y sexo.

set.seed(123)

# Simular datos
n <- 200
datos <- data.frame(
  enfermedad = rbinom(n, 1, 0.3),                  # 30% con enfermedad
  fuma = rbinom(n, 1, 0.5),                        # 50% fuma
  sexo = rbinom(n, 1, 0.5)                         # 0 = mujer, 1 = hombre
)

# asignar etiquetas
datos$sexo <- factor(datos$sexo, labels = c("Mujer", "Hombre"))
datos$fuma <- factor(datos$fuma, labels = c("No fuma", "Fuma"))

head(datos)
##   enfermedad    fuma   sexo
## 1          0 No fuma Hombre
## 2          1    Fuma  Mujer
## 3          0    Fuma Hombre
## 4          1    Fuma Hombre
## 5          1 No fuma  Mujer
## 6          0    Fuma  Mujer

Se ajusta el modelo de regresión logística, queriendo evaluar si existe un efecto modificar del sexo sobre la evaluación de la asociación entre fumar y enfermedad cardiovascular, con el término de interacción (fuma * sexo):

# Modelo con interacción
modelo <- glm(enfermedad ~ fuma * sexo, family = binomial, data = datos)

# Resumen
summary(modelo)
## 
## Call:
## glm(formula = enfermedad ~ fuma * sexo, family = binomial, data = datos)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)          -0.6325     0.3001  -2.108   0.0351 *
## fumaFuma             -0.7538     0.4638  -1.625   0.1041  
## sexoHombre           -0.1907     0.4123  -0.462   0.6438  
## fumaFuma:sexoHombre   0.4138     0.6529   0.634   0.5262  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 235.27  on 199  degrees of freedom
## Residual deviance: 231.99  on 196  degrees of freedom
## AIC: 239.99
## 
## Number of Fisher Scoring iterations: 4

Al ingresar un término de interacción, el modelo, estima los coeficiente individuales y del término multiplicativo. Ahora interpretemos los OR´s:

exp(coef(modelo))
##         (Intercept)            fumaFuma          sexoHombre fumaFuma:sexoHombre 
##           0.5312500           0.4705882           0.8263989           1.5125868

Del resultado anterior, podemos decir que:

  • La oportunidad de enfermedad cardiovascular en personas que fuman es 53% menor, independiente del sexo o en mujeres.
  • La oportunidad de enfermedad cardiovascular en hombres es 17% menor, independiente de fumar o en personas que no fuman.
  • La oportunidad de enfermedad cardiovascular en hombres que fuman es 51% mayor que en mujeres.O tambien se puede decir que, el efecto de fumar sobre la oportunidad de enfermedad cardiovascular es 51% mayor es hombres comparado con mujeres. Sin embargo, este término no fue significativo (valor p=0.5262), por tanto no hay un efecto modificador.

3. Modelo de Regresión Poisson

Siguiendo con el interés de modelar variables respuesta con distribuciones diferentes a la normal, en epidemiología muchas veces las variables que medimos hacen referencia a conteos o tasas, es aqui donde el modelo de regresión Poisson cobra relevancia.

Una variable Poisson es un conteo por unidad de tiempo o espacio, donde el promedio y la varianza es el parámetro \(\lambda\). Si lo comparamos con la regresión lineal, el parámetro de interés es la respuesta promedio, \(\mu_{i}\) y este es modelado como una función lineal de las variables explicativas. Para el modelo Poisson, también buscamos modelar el parámetro \(\lambda\) como una función lineal de otras variables, sin embargo, un modelo de la forma \(\lambda=\beta_0 + \beta_1 x_i\) no funciona bien para datos de tipo Poisson, dado que una linea podría tomar valores negativos pero \(\lambda\) no, sólo entre 0 y \(\infty\). Una forma de evitar estos problemas es modelar el \(log(\lambda)\), quedando el modelo de la forma:

\[ log(\lambda_i)=\beta_0 + \beta_1 x_i \]

Notemos que la expresión anterior del modelo poisson no contiene el término del error como en el modelo de regresión lineal, dado que, \(\lambda\) determina la media y varianza (propiedad de equidispersión de la Poisson).

Supuestos del modelo Poisson

  • La variable respuesta es un conteo por unidad de tiempo o espacio.
  • Las observaciones deben ser independientes.
  • La media es igual a la varianza.
  • El log de la media debe ser \(log(\lambda)\) debe ser una función lineal de \(x\).

3.1. Aplicación

Para este modelo se utilizará una base de datos que se encuentra en la librería ISwR:

#install.packages("ISwR")
library(ISwR)

Los datos corresponden a casos de cáncer en diferentes ciudades (datos agregados), con las siguientes variables: ciudad, edad, población y casos.

Para cargar la base de datos:

data(eba1977)
datos_cancer = eba1977
head(datos_cancer)
##         city   age  pop cases
## 1 Fredericia 40-54 3059    11
## 2    Horsens 40-54 2879    13
## 3    Kolding 40-54 3142     4
## 4      Vejle 40-54 2520     5
## 5 Fredericia 55-59  800    11
## 6    Horsens 55-59 1083     6

Para modelar casos con una distribución Poisson, se debe notar que, la base de datos tiene una forma particular, está organizada de forma agregada, por tanto variables como la ciudad se repite en la filas, para mostrar la cantidad de casos por cada grupo de edad, que es otra covariable, si se tuviera además sexo, por ejemplo, entonces tambien se debe repetir la ciudad para mostrar las combinaciones grupo de edad y sexo.

Para introducir la población como una variable offset en el modelo -la cual se incluye como predictora fija con coeficiente igual a uno y con el objetivo de ajustar las escalas de las tasas-, y calcular la razón de tasas, se tranforma a logaritmo.

logpop = log(datos_cancer[ ,3])   #se crea la variable logaritmo de la población
datos_cancer2 = cbind(datos_cancer, logpop)#se pega la nueva variable a la base de datos y se crea una nueva base.
datos_cancer2 
##          city   age  pop cases   logpop
## 1  Fredericia 40-54 3059    11 8.025843
## 2     Horsens 40-54 2879    13 7.965198
## 3     Kolding 40-54 3142     4 8.052615
## 4       Vejle 40-54 2520     5 7.832014
## 5  Fredericia 55-59  800    11 6.684612
## 6     Horsens 55-59 1083     6 6.987490
## 7     Kolding 55-59 1050     8 6.956545
## 8       Vejle 55-59  878     7 6.777647
## 9  Fredericia 60-64  710    11 6.565265
## 10    Horsens 60-64  923    15 6.827629
## 11    Kolding 60-64  895     7 6.796824
## 12      Vejle 60-64  839    10 6.732211
## 13 Fredericia 65-69  581    10 6.364751
## 14    Horsens 65-69  834    10 6.726233
## 15    Kolding 65-69  702    11 6.553933
## 16      Vejle 65-69  631    14 6.447306
## 17 Fredericia 70-74  509    11 6.232448
## 18    Horsens 70-74  634    12 6.452049
## 19    Kolding 70-74  535     9 6.282267
## 20      Vejle 70-74  539     8 6.289716
## 21 Fredericia   75+  605    10 6.405228
## 22    Horsens   75+  782     2 6.661855
## 23    Kolding   75+  659    12 6.490724
## 24      Vejle   75+  619     7 6.428105

La estimación del modelo se realiza de la sigueinte forma:

modelo3=glm(cases ~ age+ offset(logpop), family = poisson(link = "log"), data = datos_cancer2)
summary(modelo3)
## 
## Call:
## glm(formula = cases ~ age + offset(logpop), family = poisson(link = "log"), 
##     data = datos_cancer2)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.8623     0.1741 -33.676  < 2e-16 ***
## age55-59      1.0823     0.2481   4.363 1.29e-05 ***
## age60-64      1.5017     0.2314   6.489 8.66e-11 ***
## age65-69      1.7503     0.2292   7.637 2.22e-14 ***
## age70-74      1.8472     0.2352   7.855 4.00e-15 ***
## age75+        1.4083     0.2501   5.630 1.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 129.908  on 23  degrees of freedom
## Residual deviance:  28.307  on 18  degrees of freedom
## AIC: 136.69
## 
## Number of Fisher Scoring iterations: 5

Así como en el modelo logístico obteniamos los coeficientes en logaritmo, igual con el modelo Poisson, para calcular los riesgos relativos (RR) o razones de tasas, dado que se tiene la población expuesta, debemos exponenciar los beta´s.

Para obtener los Riesgos Relativos y sus intervalos de confianza:

exp(coef(modelo3)) # exponenciar los coeficientes
## (Intercept)    age55-59    age60-64    age65-69    age70-74      age75+ 
## 0.002844828 2.951583534 4.489204489 5.756252481 6.342176843 4.088919211
exp(confint(modelo3)) # exponenciar los límites del intervalo de confianza
## Waiting for profiling to be done...
##                   2.5 %       97.5 %
## (Intercept) 0.001981286  0.003928784
## age55-59    1.809509477  4.808601155
## age60-64    2.859538726  7.111934189
## age65-69    3.685432697  9.084775159
## age70-74    4.005487299 10.110915401
## age75+      2.495074580  6.684025923

Ahora vamos a interpretar los RR del resultado anterior, teniendo en cuenta que la edad es una variable categórica, se generan \(p-1\) coeficientes, siendo \(p\) el número de categorías de la variable. La interpretación se realiza de la siguiente forma, teniendo en cuenta que se realizó la estimación de un modelo simple, sólo con la edad.

  • La tasa de cancer en personas de 55-59 años es 1.80 veces la tasa del grupo entre 40-54 años.
  • La tasa de cancer en personas de 60-64 años es 2.86 veces la tasa del grupo entre 40-54 años.
  • La tasa de cancer en personas de 65-69 años es 3.68 veces la tasa del grupo entre 40-54 años.

4. Modelo de Regresión de Cox

El análisis de sobrevida es muy usado en estudios de cohorte o ensayos clínicos, donde se miden tiempos hasta que se presenta un evento (muerte, recaída, enfermedad, etc). Aqui el interés está en medir el riesgo dependiente del tiempo, sin embargo, se debe tener en cuenta que en estudios longitudinales generalmente no se tienen seguimientos completos, se presentan observaciones censuradas, además los tiempos generalmente tienen distribuciones asimétricas, no normales.

De esta forma el análisis de sobrevida y el modelo cox permite la comparación completa de sobrevida surante el seguimiento para cauqlueir variable con desenlace dicotómico. Revisar la conceptualización del análisis de sobrevida y regresión de cox en los siguientes textos:

https://www.fisterra.com/formacion/metodologia-investigacion/analisis-supervivencia/ https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8651375/pdf/OMCL2021-1302811.pdf

Para este análisis de sobrevida se utilizará una base de datos sobre un ensayo clínico de dos tratamientos para cancer de pulmón en veteranos de guerra de la librería survival.

library(survival)
library(dplyr)
library(ggplot2)
library(ggfortify)

Las variables contenidas en la base de datos son:

  • trt: 1=estándar 2=prueba
  • celltype: 1=escamoso, 2=pequeña, 3=adeno, 4=grande
  • time: tiempo de sobrevida en días
  • status: censura (0=censura, 1=muerte)
  • karno: puntaje Karnofsky (100=bueno)
  • diagtime: meses desde el diagnóstico hasta la aleatorización
  • age: edad en años
  • prior: terapia previa 0=no, 10=si.

Para ver el encabezado de la base de datos:

data(veteran)
head(veteran)
##   trt celltype time status karno diagtime age prior
## 1   1 squamous   72      1    60        7  69     0
## 2   1 squamous  411      1    70        5  64    10
## 3   1 squamous  228      1    60        3  38     0
## 4   1 squamous  126      1    60        9  63    10
## 5   1 squamous  118      1    70       11  65    10
## 6   1 squamous   10      1    20        5  49     0

Inicialmente se analizará la tabla de sobrevida de Kaplan-Meier. Esta permite calcular la sobrevida para cada punto en el tiempo en que se presenta un evento, además de la función de sobrevida acumulada.

km_fit <- survfit(Surv(time, status) ~ 1, data=veteran)  
#se debe incluir la variable tiempo y el evento/censura

summary(km_fit, times = c(1,30,60,90*(1:10)))
## Call: survfit(formula = Surv(time, status) ~ 1, data = veteran)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    137       2    0.985  0.0102      0.96552       1.0000
##    30     97      39    0.700  0.0392      0.62774       0.7816
##    60     73      22    0.538  0.0427      0.46070       0.6288
##    90     62      10    0.464  0.0428      0.38731       0.5560
##   180     27      30    0.222  0.0369      0.16066       0.3079
##   270     16       9    0.144  0.0319      0.09338       0.2223
##   360     10       6    0.090  0.0265      0.05061       0.1602
##   450      5       5    0.045  0.0194      0.01931       0.1049
##   540      4       1    0.036  0.0175      0.01389       0.0934
##   630      2       2    0.018  0.0126      0.00459       0.0707
##   720      2       0    0.018  0.0126      0.00459       0.0707
##   810      2       0    0.018  0.0126      0.00459       0.0707
##   900      2       0    0.018  0.0126      0.00459       0.0707

Al final del seguimiento la sobrevida es del 7%, y la probabilidad del evento (incidencia) es del 93%. Esta función se puede observar graficamente:

Gráfico de Kaplan-Meier:

autoplot(km_fit)

Interesaría observar la función de sobrevida diferenciada por tratamiento.

Gráfico de Kaplan-Meier diferenciado por tratamiento

km_trt_fit <- survfit(Surv(time, status) ~ trt, data=veteran)
autoplot(km_trt_fit)

En el gráfico anterior podemos observar que las curvas de sobrevida (con intervalos de confianza) diferenciadas por tratamiento se cruzan, por tanto, muy posiblemente no se encuentren diferencias entre tratamientos. Al aplicar la prueba para la comparación de curvas (log rank test), se observa lo siguiente:

surv_obj <- Surv(time = veteran$time, event = veteran$status == 1)
log_rank <- survdiff(surv_obj ~ trt, data = veteran)
log_rank
## Call:
## survdiff(formula = surv_obj ~ trt, data = veteran)
## 
##        N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=1 69       64     64.5   0.00388   0.00823
## trt=2 68       64     63.5   0.00394   0.00823
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9

Con un valor p de 0.9, mayor que el nivel de significancia de 0.05, no existen diferencias en las curvas de sobrevida entre tratamientos.

Luego, para la estimación del modelo cox se utiliza la función \(coxph\):

cox <- coxph(Surv(time, status) ~ trt + celltype + karno + diagtime + age + prior , data = veteran)
summary(cox)
## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno + 
##     diagtime + age + prior, data = veteran)
## 
##   n= 137, number of events= 128 
## 
##                         coef  exp(coef)   se(coef)      z Pr(>|z|)    
## trt                2.946e-01  1.343e+00  2.075e-01  1.419  0.15577    
## celltypesmallcell  8.616e-01  2.367e+00  2.753e-01  3.130  0.00175 ** 
## celltypeadeno      1.196e+00  3.307e+00  3.009e-01  3.975 7.05e-05 ***
## celltypelarge      4.013e-01  1.494e+00  2.827e-01  1.420  0.15574    
## karno             -3.282e-02  9.677e-01  5.508e-03 -5.958 2.55e-09 ***
## diagtime           8.132e-05  1.000e+00  9.136e-03  0.009  0.99290    
## age               -8.706e-03  9.913e-01  9.300e-03 -0.936  0.34920    
## prior              7.159e-03  1.007e+00  2.323e-02  0.308  0.75794    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## trt                  1.3426     0.7448    0.8939    2.0166
## celltypesmallcell    2.3669     0.4225    1.3799    4.0597
## celltypeadeno        3.3071     0.3024    1.8336    5.9647
## celltypelarge        1.4938     0.6695    0.8583    2.5996
## karno                0.9677     1.0334    0.9573    0.9782
## diagtime             1.0001     0.9999    0.9823    1.0182
## age                  0.9913     1.0087    0.9734    1.0096
## prior                1.0072     0.9929    0.9624    1.0541
## 
## Concordance= 0.736  (se = 0.021 )
## Likelihood ratio test= 62.1  on 8 df,   p=2e-10
## Wald test            = 62.37  on 8 df,   p=2e-10
## Score (logrank) test = 66.74  on 8 df,   p=2e-11

Interpretando la exponencial de los coeficientes para quitar el logaritmo, se estiman los *Hazard ratios, su interpretación es: - La probabilidad de morir en el grupo con el tratamiento de prueba es 1.34 veces la probabilidad del tratamiento estándar para todo momento del seguimiento** independiente del efecto de las demás covariables presentes en el modelo (valor p>0.05)

En el caso de la edad que es una variable cuantitativa: - Por cada año de aumento en la edad, la probabilidad de morir disminuye en 1% para todo momento del seguimiento independiente del efecto de las demás covariables presentes en el modelo (valor p>0.05).

La función de sobrevida estimada por el modelo:

cox_fit <- survfit(cox)
autoplot(cox_fit)

Para el modelo de cox se debe validar el supuesto de proporcionalidad, este evalúa que el riesgo en la exposición de interés multiplica el riesgo de base por un factor constante en cualquier punto durante el seguimiento. Para la validación del supuesto de proporcionalidad se utiliza la siguiente función:

test.ph <- cox.zph(cox)
test.ph
##            chisq df       p
## trt       0.2644  1 0.60712
## celltype 15.2274  3 0.00163
## karno    12.9352  1 0.00032
## diagtime  0.0129  1 0.90961
## age       1.8288  1 0.17627
## prior     2.1656  1 0.14113
## GLOBAL   34.5525  8 3.2e-05

Al interpretar el valor global, la hipótesis nula plantea la proporcionalidad, teniendo en cuenta el valor de p, para este modelo no se cumple ese supuesto, por tanto se debe reconsiderar el modelo.

Por último, la validación a través de los residuales de Schoenfeld se pueden obtener de la siguiente forma:

#install.packages("survminer")
library("survminer")
ggcoxzph(test.ph)

El gráfico anterior se obtiene por cada variable explicativa del modelo, y el comportamiento que se espera no observar una tendencia y que la linea no tenga una pendiente, es decir, que no dependa del tiempo.

5. Referencias

Para la preparación de este material se utilizaron las siguientes referencias: