Regresión logística condicional

Autor/a

Tamara Ricardo

Introducción

En algunos estudios de casos y controles, se utiliza un diseño apareado para controlar la influencia de factores de confusión como la edad, el sexo o el lugar de residencia. La regresión logística condicional (RLC) es una extensión de la regresión logística clásica que incorpora el apareamiento de los datos, estratificando el análisis por grupos definidos por las variables de emparejamiento (estratos).

Este modelo permite estimar el efecto de las covariables sobre una variable dependiente binaria mientras controla de manera efectiva los efectos de las variables de emparejamiento o agrupamiento. Es especialmente útil en estudios de casos y controles apareados, investigaciones multicéntricas y estudios observacionales que abarcan múltiples ubicaciones geográficas, proporcionando un enfoque robusto para manejar la heterogeneidad dentro de los datos.

Modelo logístico condicional

El modelo logístico condicional utiliza una función de verosimilitud condicional que elimina la necesidad de estimar los efectos de los factores de emparejamiento como interceptos adicionales. Este enfoque resulta especialmente útil cuando hay un gran número de grupos o estratos, ya que permite concentrarse exclusivamente en los coeficientes de las covariables de interés.

Se modela la probabilidad de que un evento ocurra en función de las covariables dentro de cada estrato, definido por las variables de emparejamiento. La función de verosimilitud condicional se puede expresar matemáticamente como:

\[ logit(p) = \alpha_1 + \alpha_2 z_2 + ... + \alpha_s z_s + \beta_1 x_1 + ... + \beta_p x_p \]

donde:

  • \(\alpha_1 + \alpha_2 z_2 + ... + \alpha_s z_s\): Coeficientes asociados a las variables de estratificación \(z_s\), que identifican los estratos.

  • \(z_1, z_2, ..., z_s\): Variables de emparejamiento o estratos.

  • \(\beta_1, ..., \beta_p\): Efectos de las variables independientes ajustadas por los estratos.

  • \(x_1, ... x_p\): Variables independientes incluídas en el modelo.

El uso de la verosimilitud condicional simplifica el modelo al evitar la estimación de un gran número de interceptos y reduce el sesgo en situaciones con pocos datos por estrato.

Supuestos de la regresión logística condicional

  • Independencia dentro de los estratos: Las observaciones son independientes entre los diferentes estratos o grupos definidos por el emparejamiento, pero pueden estar correlacionadas dentro de los mismos.

  • Relación logit-lineal: Existe una relación lineal entre el logit de la probabilidad del evento y las covariables explicativas.

  • No existencia de confusión residual: Las variables de emparejamiento deben ser suficientes para controlar por los factores de confusión.

  • Tamaño adecuado del estrato: Cada estrato debe tener al menos un caso y un control.

Ejemplo práctico en R

En R podemos ajustar modelos de regresión logística condicional usando la función clogit() del paquete survival (Therneau 2024).

Cargamos los paquetes necesarios:

# Análisis estadístico
library(gtsummary)
library(survival)

# Chequeo de supuestos
library(easystats)

# Apareamiento de datos
library(MatchIt)

# Manejo y exploración de datos
library(skimr)
library(janitor)
library(tidyverse)

Usaremos la base “higado_graso.txt” que contiene datos de un estudio de casos y controles apareados por sexo y edad para muertes por hígado graso no alcohólico:

# Carga datos
datos <- read_delim("datos/higado_graso.txt")|> 
  
  # Cambia nivel de referencia de fib4
  mutate(fib4_cat = fct_rev(fib4_cat))

# Explora datos
glimpse(datos)
Rows: 2,467
Columns: 11
$ id           <dbl> 19, 25, 26, 32, 56, 57, 68, 77, 78, 99, 100, 120, 126, 13…
$ tipo         <dbl> 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, …
$ edad         <dbl> 35, 39, 55, 45, 51, 81, 51, 55, 32, 33, 50, 82, 56, 69, 7…
$ sexo         <dbl> 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, …
$ bmi          <dbl> 30.49939, 24.26164, 30.59225, 22.11900, 36.42498, 26.3815…
$ hdl          <dbl> 50.09091, 52.50000, 49.13043, 70.55556, 45.10526, 67.7500…
$ p_sistolica  <dbl> NA, NA, 134.2500, NA, 164.0000, 148.0833, 124.5000, NA, N…
$ p_diastolica <dbl> NA, NA, 134.2500, NA, 164.0000, 148.0833, 124.5000, NA, N…
$ fib4_cat     <fct> bajo riesgo, bajo riesgo, bajo riesgo, bajo riesgo, bajo …
$ fuma         <dbl> 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, …
$ status       <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, …

El dataset tiene 2467 observaciones y 11 columnas. Las variables de interés son las siguientes:

  • id: identificador único del paciente.

  • tipo: Indica si el paciente tiene hígado graso no alcohólico (caso: 1) o no presenta la condición (control: 0).

  • edad: Edad del paciente al momento de su ingreso al hospital.

  • sexo: Sexo biológico del paciente (1: masculino, 0: femenino).

  • bmi: Índice de masa corporal del paciente al momento de la hospitalización.

  • hdl: Nivel de lipoproteínas de alta densidad (HDL) en sangre.

  • p_sistolica: Presión sistólica en milímetros de mercurio (mmHg).

  • p_diastolica: Presión diastólica en milímetros de mercurio (mmHg).

  • fib4_cat: Índice de fibrosis hepática categorizado (bajo riesgo, alto riesgo).

  • fuma: Tabaquismo (0: no, 1: sí).

  • status: Variable dependiente que toma el valor de 0 si la persona sobrevivió y 1 si falleció.

Exploremos más en profundidad los datos usando la función tabyl() del paquete janitor y la función skim() del paquete skimr:

# Número de casos y controles
tabyl(datos$tipo) |> 
  adorn_pct_formatting()
 datos$tipo    n percent
          0 1685   68.3%
          1  782   31.7%
# Niveles fib4
tabyl(datos$fib4_cat)
 datos$fib4_cat    n   percent
    bajo riesgo 2109 0.8548845
    alto riesgo  358 0.1451155
# Explora valores ausentes
skim(datos)
Data summary
Name datos
Number of rows 2467
Number of columns 11
_______________________
Column type frequency:
factor 1
numeric 10
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
fib4_cat 0 1 FALSE 2 baj: 2109, alt: 358

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
id 0 1.00 8644.27 4998.09 19.00 4486.00 8507.00 12864.00 17563.0 ▇▇▇▇▇
tipo 0 1.00 0.32 0.47 0.00 0.00 0.00 1.00 1.0 ▇▁▁▁▃
edad 0 1.00 53.60 14.73 18.00 44.00 53.00 64.00 98.0 ▂▆▇▃▁
sexo 0 1.00 0.50 0.50 0.00 0.00 0.00 1.00 1.0 ▇▁▁▁▇
bmi 0 1.00 31.22 7.51 9.21 26.02 30.14 35.04 84.4 ▂▇▁▁▁
hdl 0 1.00 50.92 14.94 19.83 40.00 48.53 58.50 134.0 ▅▇▂▁▁
p_sistolica 928 0.62 131.41 17.11 40.00 121.53 132.50 142.67 211.0 ▁▁▇▃▁
p_diastolica 928 0.62 131.41 17.11 40.00 121.53 132.50 142.67 211.0 ▁▁▇▃▁
fuma 0 1.00 0.64 0.48 0.00 0.00 1.00 1.00 1.0 ▅▁▁▁▇
status 0 1.00 0.12 0.32 0.00 0.00 0.00 0.00 1.0 ▇▁▁▁▁

Las variables p_sistolica y p_diastolica tienen gran cantidad de valores ausentes, por lo que no serán consideradas para el análisis.

Como no tenemos definido a que caso corresponde cada control, vamos a aparear casos y controles usando la función matchit() del paquete MatchIt (Ho et al. 2011), especificando que queremos dos controles por caso con el argumento ratio = 2:

# Aparea casos y controles por sexo y edad
match <- matchit(tipo ~ edad + sexo, 
                 data = datos, 
                 exact = "sexo",
                 ratio = 2)

# Resumen del apareamiento
summary(match)

Call:
matchit(formula = tipo ~ edad + sexo, data = datos, exact = "sexo", 
    ratio = 2)

Summary of Balance for All Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.3268        0.3124          0.2552     1.0589    0.0507
edad           51.0243       54.7923         -0.2550     1.0302    0.0486
sexo            0.4962        0.5003         -0.0083          .    0.0041
         eCDF Max
distance   0.1284
edad       0.1274
sexo       0.0041

Summary of Balance for Matched Data:
         Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance        0.3268        0.3181          0.1549     1.1495    0.0301
edad           51.0243       53.2270         -0.1490     1.1373    0.0288
sexo            0.4962        0.4962          0.0000          .    0.0000
         eCDF Max Std. Pair Dist.
distance   0.0882          0.1578
edad       0.0838          0.1528
sexo       0.0000          0.0000

Sample Sizes:
          Control Treated
All          1685     782
Matched      1564     782
Unmatched     121       0
Discarded       0       0

Obtenemos los datos apareados:

datos <- match.data(match, drop.unmatched = T) 

names(datos)
 [1] "id"           "tipo"         "edad"         "sexo"         "bmi"         
 [6] "hdl"          "p_sistolica"  "p_diastolica" "fib4_cat"     "fuma"        
[11] "status"       "distance"     "weights"      "subclass"

Ajustamos el modelo de regresión logística condicional, usando la variable subclass como variable de estratificación:

mod_cond <- clogit(status ~ bmi + hdl + fuma + fib4_cat +
                     strata(subclass),
                   data = datos)

Revisemos la salida del modelo:

summary(mod_cond)
Call:
coxph(formula = Surv(rep(1, 2346L), status) ~ bmi + hdl + fuma + 
    fib4_cat + strata(subclass), data = datos, method = "exact")

  n= 2346, number of events= 251 

                         coef exp(coef)  se(coef)      z Pr(>|z|)    
bmi                  0.005507  1.005522  0.012814  0.430    0.667    
hdl                 -0.006502  0.993519  0.006980 -0.931    0.352    
fuma                -0.237017  0.788978  0.186629 -1.270    0.204    
fib4_catalto riesgo  1.072163  2.921692  0.224099  4.784 1.72e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
bmi                    1.0055     0.9945    0.9806     1.031
hdl                    0.9935     1.0065    0.9800     1.007
fuma                   0.7890     1.2675    0.5473     1.137
fib4_catalto riesgo    2.9217     0.3423    1.8831     4.533

Concordance= 0.625  (se = 0.039 )
Likelihood ratio test= 27.94  on 4 df,   p=1e-05
Wald test            = 25.92  on 4 df,   p=3e-05
Score (logrank) test = 28.23  on 4 df,   p=1e-05

Probemos eliminar una de las variables explicativas no significativas a la vez:

# (-) Tabaquismo
mod_cond1 <- clogit(status ~ bmi + hdl + fib4_cat +
                     strata(subclass),
                   data = datos)

# (-) IMC
mod_cond2 <- clogit(status ~ hdl + fuma + fib4_cat +
                     strata(sexo, edad),
                   data = datos)

# (-) HDL
mod_cond3 <- clogit(status ~ bmi + fuma +  fib4_cat +
                     strata(sexo, edad),
                   data = datos)

Comparamos los modelos

# Compara modelos
compare_performance(mod_cond, mod_cond1, mod_cond2, mod_cond3,
                    metrics = "common")
# Comparison of Model Performance Indices

Name      |  Model | AIC (weights) | BIC (weights) | Nagelkerke's R2 |  RMSE
----------------------------------------------------------------------------
mod_cond  | clogit | 393.1 (0.452) | 416.2 (0.044) |           0.073 | 0.230
mod_cond1 | clogit | 392.8 (0.548) | 410.0 (0.956) |           0.069 | 0.230
mod_cond2 | clogit | 921.1 (<.001) | 938.3 (<.001) |           0.055 | 0.276
mod_cond3 | clogit | 923.5 (<.001) | 940.8 (<.001) |           0.052 | 0.276
# Compara modelos y ordena según performance
compare_performance(mod_cond, mod_cond1, mod_cond2, mod_cond3,
                    metrics = "common", rank = T)
# Comparison of Model Performance Indices

Name      |  Model | Nagelkerke's R2 |  RMSE | AIC weights | BIC weights | Performance-Score
--------------------------------------------------------------------------------------------
mod_cond1 | clogit |           0.069 | 0.230 |       0.548 |       0.956 |            95.02%
mod_cond  | clogit |           0.073 | 0.230 |       0.452 |       0.044 |            71.65%
mod_cond2 | clogit |           0.055 | 0.276 |   1.04e-115 |   1.81e-115 |             3.55%
mod_cond3 | clogit |           0.052 | 0.276 |   3.12e-116 |   5.45e-116 |             0.40%

El modelo sin fuma tiene mejor performance que el modelo saturado:

# Salida modelo condicional 1
summary(mod_cond1)
Call:
coxph(formula = Surv(rep(1, 2346L), status) ~ bmi + hdl + fib4_cat + 
    strata(subclass), data = datos, method = "exact")

  n= 2346, number of events= 251 

                         coef exp(coef)  se(coef)      z Pr(>|z|)    
bmi                  0.004445  1.004455  0.012740  0.349    0.727    
hdl                 -0.006922  0.993102  0.006941 -0.997    0.319    
fib4_catalto riesgo  1.064462  2.899278  0.223595  4.761 1.93e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
bmi                    1.0045     0.9956    0.9797     1.030
hdl                    0.9931     1.0069    0.9797     1.007
fib4_catalto riesgo    2.8993     0.3449    1.8705     4.494

Concordance= 0.641  (se = 0.038 )
Likelihood ratio test= 26.33  on 3 df,   p=8e-06
Wald test            = 24.56  on 3 df,   p=2e-05
Score (logrank) test = 26.63  on 3 df,   p=7e-06
# (-) BMI
mod_cond1a <- clogit(status ~ hdl + fib4_cat +
                     strata(subclass),
                   data = datos)

# (-) HDL
mod_cond1b <- clogit(status ~ bmi + fib4_cat +
                     strata(subclass),
                   data = datos)

# Compara modelos y ordena según performance
compare_performance(mod_cond1, mod_cond1a, mod_cond1b, 
                    metrics = "common", rank = T)
# Comparison of Model Performance Indices

Name       |  Model | Nagelkerke's R2 |  RMSE | AIC weights | BIC weights | Performance-Score
---------------------------------------------------------------------------------------------
mod_cond1a | clogit |           0.069 | 0.230 |       0.492 |       0.601 |            72.00%
mod_cond1  | clogit |           0.069 | 0.230 |       0.192 |       0.013 |            50.00%
mod_cond1b | clogit |           0.066 | 0.230 |       0.316 |       0.386 |            32.90%
# Salida modelo condicional 1b
summary(mod_cond1b)
Call:
coxph(formula = Surv(rep(1, 2346L), status) ~ bmi + fib4_cat + 
    strata(subclass), data = datos, method = "exact")

  n= 2346, number of events= 251 

                        coef exp(coef) se(coef)     z Pr(>|z|)    
bmi                 0.007979  1.008011 0.012127 0.658    0.511    
fib4_catalto riesgo 1.067510  2.908129 0.223684 4.772 1.82e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
bmi                     1.008     0.9921    0.9843     1.032
fib4_catalto riesgo     2.908     0.3439    1.8759     4.508

Concordance= 0.601  (se = 0.04 )
Likelihood ratio test= 25.32  on 2 df,   p=3e-06
Wald test            = 23.74  on 2 df,   p=7e-06
Score (logrank) test = 25.68  on 2 df,   p=3e-06
# (-) bmi
mod_cond1c <- clogit(status ~ fib4_cat +
                     strata(subclass),
                   data = datos)

# Compara modelos
compare_performance(mod_cond1b, mod_cond1c,
                    metrics = "common")
# Comparison of Model Performance Indices

Name       |  Model | AIC (weights) | BIC (weights) | Nagelkerke's R2 |  RMSE
-----------------------------------------------------------------------------
mod_cond1b | clogit | 391.8 (0.313) | 403.3 (0.025) |           0.066 | 0.230
mod_cond1c | clogit | 390.2 (0.687) | 396.0 (0.975) |           0.065 | 0.230

Nos quedamos con el modelo que tiene solo fib4_cat como variable explicativa:

tbl_regression(mod_cond1c, exponentiate = T)
Characteristic OR 95% CI p-value
fib4_cat


    bajo riesgo — —
    alto riesgo 2.94 1.90, 4.55 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

En el modelo final, las personas con alto riesgo de fibrosis hepática presentan 3,7 veces más probabilidades de morir (\(95% IC: 2,45-5,57\)) que aquellas sin la condición.

Volver arriba

Referencias

Agresti, Alan. 2015. Foundations of Linear and Generalized Linear Models. Wiley Series En Probability y Statistics. Hoboken, New Jersey: John Wiley & Sons, Inc.
Ho, Daniel E., Kosuke Imai, Gary King, y Elizabeth A. Stuart. 2011. «MatchIt: Nonparametric Preprocessing for Parametric Causal Inference» 42. https://doi.org/10.18637/jss.v042.i08.
Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Brenton M. Wiernik, Etienne Bacher, Rémi Thériault, y Dominique Makowski. 2022. «easystats: Framework for Easy Statistical Modeling, Visualization, and Reporting». https://doi.org/10.32614/CRAN.package.easystats.
Therneau, Terry M. 2024. «A Package for Survival Analysis in R». https://CRAN.R-project.org/package=survival.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. «Welcome to the tidyverse» 4: 1686. https://doi.org/10.21105/joss.01686.

Reutilización