Un modelo de riesgos proporcionales ajustado a los datos del COVID-19 en México

Introducción

El coronavirus SARS-Cov-2 es un virus cuya primera aparición fue en Wuhan, China en diciembre de 2019. Este virus provoca la enfermedad COVID-19. En México, el virus fue detectado por primera vez a finales de febrero de 2020. Para el 15 de mayo del mismo año, se registraron un total de 45,032 casos confirmados, 4,767 defunciones y 37,325 recuperados de dicha enfermedad ¹Gobierno de México. Coronavirus.. A nivel mundial, para el 15 de mayo de 2020 el total de casos confirmados era de 4,347,935 mientras que las muertes registradas ascendían a 297,236²WHO. Who coronavirus disease $($covid-19) dashboard.. El presente artículo resume el proceso y resultados obtenidos al realizar un análisis de supervivencia basado en la información de los pacientes de coronavirus recabada por el gobierno de México. Nuestros principales objetivos fueron identificar las principales condiciones que aumentan el riesgo de fallecimiento por causas relacionadas directamente con el coronavirus y establecer las relaciones existentes entre las características individuales de los infectados y el desenlace del padecimiento.
El trabajo en cuestión fue elaborado con los https://datos.gob.mx/busca/dataset/informacion-referente-a-casos-covid-19-en-mexico datos liberados por el gobierno de México asociados a COVID-19 hasta el 15 de mayo de 2020 por la Dirección General de Epidemiología de la Secretaría de Salud³Gobierno de México. Coronavirus..

Variables y escala de medición

Al 15 de mayo del 2020, la base de datos contaba con 163,691 observaciones de 34 variables. Las siguientes variables son las que elegimos para llevar a cabo el análisis de supervivencia.

• Sexo: identifíca el sexo del paciente. Valores:
  ·0 = Hombre                              ·1 = Mujer
• Días internado: identifíca el número de días que lleva internado el paciente; esta variable la construimos nosotros.
• Edad: identifica la edad del paciente medida en años.
• Delta: identifica el tipo de observación; esta variable la construimos nosotros. Valores:
  ·0 = censurada                              ·1 = exacta .
  Las variables presentadas a continuación son binarias y han sido codificadas para tomar los valores:
  ·0 = No                              ·1 = Sí
• Neumonía: identifica si al paciente se le diagnosticó neumonía.
• Diabetes: identifica si el paciente tiene un diagnóstico de diabetes.
• EPOC: identifica si el paciente tiene un diagnóstico de EPOC.
• Asma: identifica si el paciente tiene un diagnóstico de asma.
• Inmunosupresión: identifica si el paciente presenta inmunosupresión.
• Hipertensión: identifica si el paciente tiene un diagnóstico de hipertensión.
• Cardiovascular: identifica si el paciente tiene un diagnóstico de enfermedades cardiovasculares.
• Obesidad: identifica si el paciente tiene un diagnóstico de obesidad.
• Renal crónica: identifica si el paciente tiene un diagnóstico de insuficiencia renal crónica.
• Tabaquismo: identifica si el paciente tiene hábito de tabaquismo.

Aunque la base de datos contaba con información acerca de todos los pacientes a los cuales se les hubiese aplicado una prueba para coronavirus al igual que la fecha en la que se presentaron síntomas por primera vez, se optó por descartar todas las observaciones que no fueran positivas para el virus al igual que la fecha de primeros síntomas. Esto último bajo la lógica que la respuesta que da alguien cuando se le pregunta cuándo mostró síntomas por primera vez es altamente subjetiva, pues no todos padecen de los mismos síntomas en la misma gravedad o cronología. Por lo tanto, se optó por tomar como fecha de entrada al estudio la fecha en la que el paciente fue ingresado con una prueba positiva a COVID-19, considerando que distintos pacientes presentarían distintos niveles de progreso de la enfermedad. De esta forma, la muestra obtenida presenta censura generalizada del tipo I por la derecha, pues ante distintas fechas de entrada al estudio, la fecha de censura $($la fecha de actualización de la base de datos$)$ hace que la cantidad de días que llevan en el estudio los pacientes sea distinta. Cabe recalcar que el planteamiento del trabajo no considera truncamiento, pues los objetivos del estudio son solamente aquellos pacientes que ya tengan un diagnostico positivo. Por la misma razón, las conclusiones del estudio aplicarán solamente a individuos que sean COVID-19 positivos.

Análisis exploratorio

Si una enfermedad no incide considerablemente en el riesgo de fallecer para un paciente dado, entonces las proporciones de los pacientes positivos a la condición deberían ser similares entre aquellos pacientes que aun no han fallecido y aquellos que si. Bajo esta lógica, resultó interesante conocer el porcentaje de individuos COVID-19 positivos que tenían una enfermedad subyacente y el porcentaje entre los difuntos COVID-19 positivo que sufrían de estas condiciones con el objetivo de comparar la prevalencia de estas condiciones entre las dos poblaciones.
Al tomar la razón de incidencia entre los pacientes COVID-19 positivos que fallecieron teniendo el virus y aquellos que son COVID-19 positivos, es claro que un número mayor a 1 indica que la incidencia es mayor entre lo fallecidos, mientras que un valor cercano a 1 indica que no hay gran cambio en las proporciones. Definiendo a $p_{C_i}$ como los pacientes de coronavirus que sufrían la condición $i$, y a $m_{C_i}$ como los fallecidos que padecían la condición $i$, representamos a continuación gráficamente el cociente $\frac{m_{C_i}}{p_{C_i}}$ para las distintas condiciones.

La Figura 1 indica que la población de pacientes con problemas renales crónicos que murió fue 3 veces mayor que aquella que seguía viva al 15 de mayo. Para aquellos con EPOC esta razón fue de 2.7, y con neumonía 2.5 veces. De forma preliminar, esto parece identificar a estas 3 condiciones como candidatas a ser factores de riesgo para aquellos que padecen de coronavirus. Por otro lado, podemos notar que la obesidad y tabaquismo, al tener valores cercanos a 1, no cambian su presencia en la población enferma y los difuntos. El asma fue la única condición subyacente en los pacientes de COVID-19 que presentó una menor cantidad de población muerta que viva.

Otro aspecto que resulta interesante analizar es la distribución de muertes por edad.
Observamos que la mayor cantidad de defunciones por coronavirus suceden en pacientes entre los 50 y 70 años con una media de 59.2 años. Para los individuos entre 0 y 25 años y los mayores a 80 no se observan tantas muertes. Sin embargo, esto puede ser por la estructura de la pirámide de población en el país y la forma en la que interactúan miembros de estos distintos estratos generacionales. También cabe recalcar que la distribución de muertes por edad en México es distinta a la de otros países, como Italia, donde la media es de 85 años⁴WHO. Who coronavirus disease $($covid-19) dashboard..
Teniendo en mente las limitaciones que conlleva contar a la cantidad de fallecidos por edad, consideramos la tasa de mortalidad observada cuyos resultados se presentan en la Figura 2. Se observa cómo de los 0 a 10 años, la tasa de mortalidad es mayor que de los 10 a los 40 años, mientras que a partir de los 30 años, la tasa de mortalidad comienza a aumentar y a partir de los 60 años el valor es superior a 0.2, alcanzando el máximo valor de dicha tasa en el intervalo de 90 a 100 años, lo cual sugiere que la edad puede ser un factor de riesgo para aquellos que padecen de coronavirus.

Análisis inferencial

A partir de la información proporcionada por la Figura 1, decidimos analizar las funciones de supervivencia y de riesgo acumulado para pacientes de diferente sexo y posteriormente para las condiciones de fallas renales crónicas y tabaquismo. Para esto, ajustamos la función de supervivencia utilizando el estimador Kaplan-Meier para modelar el tiempo de muerte de un paciente con coronavirus y el estimador Nelson-Aalen para modelar la función de riesgo acumulada. En la Figura 3 se puede observar que la función de supervivencia ajustada para los hombres toma valores más pequeños que para las mujeres, mientras que de manera esperada figura 4, muestra que la función de riesgo ajustada es mayor para los hombres que para las mujeres.

La condición de fallos renales crónicos puesto que 1 es la que tiene una mayor razón de cambio entre los difuntos y pacientes por coronavirus. La diferencia en los valores de la función de supervivencia entre pacientes con fallas renales crónicas y los que no tienen, es de casi 30 puntos porcentuales a partir del día 50. El valor mínimo que alcanza la función de supervivencia para individuos que no sufren fallas renales crónicas es 0.87, mientras que para los que sí la sufren es de 0.59. En cuanto a la función de riesgo acumulado, para pacientes sin fallos renales crónicas, los valores son inferiores a 0.2 y, en cambio, para los pacientes que sí tienen esa condición, los valores son superiores a 0.50 a partir del día 40.

Nos pareció importante ahondar en el padecimiento de tabaquismo debido a que no hay una diferencia tan fuerte entre la prevalencia de esta entre pacientes vivos y difuntos contagiados de coronavirus. Resulta ser que las curvas de supervivencia para individuos COVID-19 positivos que fuman y que no fuman son muy parecidas pues la máxima diferencia entre sus respectivas funciones es de 0.02. Para individuos que fuman la probabilidad de sobrevivir para el día 25 es del 86% mientras que para los que no fuman es del 87.5%.

Mediante la función survdiff en R se aplicó una prueba de log-rank para probar la hipótesis nula de que las curvas de cada grupo son iguales, para cada variable. La función regresó, no sorprendentemente, significancias observadas muy por debajo de la significancia de la prueba del 5%, por lo que en cada instancia se rechaza la hipótesis nula de que los estimadores de Kaplan-Meier, al estratificar, son iguales.

Un modelo de riesgo

Para investigar la asociación que hay entre características propias a individuos $($sexo, edad y condiciones de salud subyacentes$)$ y los distintos niveles de riesgo a morir a causa de COVID-19, ajustamos el Modelo de Riesgos Proporcionales de Cox⁵David R Cox. Regression models and life-tables.Journal of the Royal Statistical Society:Series B$($Methodological), 34$($2):187–202, 1972. a los datos.
Supongamos que en una muestra de $n$ individuos se midieron $p$ distintos factores que podrían afectar la supervivencia $($conocidos como covariables$)$. Entonces denotamos a la muestra de un individuo como $\mathbf{X} = (X_{1}, X_{2}, \ldots, X_{p})$. Con esto, el Modelo de Riesgos Proporcionales es de la forma

\begin{align}
\label{modeloCox}
h(t | \mathbf{X}) = h_0(t) \exp{\sum_{i=1}^p \beta_i X_i },
\end{align}

en donde $t$ representa al tiempo, $h_0(t)$ es la función de riesgo base y los coeficientes $\beta_0, \ldots, \beta_p$ miden el impacto de las covariables.
El modelo recibe su nombre del hecho que el cociente de la función de riesgos para dos individuos con covariables respectivas $\mathbf{X},\mathbf{Z}$ es la constante

\begin{align}
\label{cocienteRiesgos}
\frac{h(t| \mathbf{X})}{h(t| \mathbf{Z})} = \exp{\sum_{i=1}^p \beta_i(X_i – Z_i) }.
\end{align}

En el caso de tener covariables que sean categóricas, el representar esta información por medio de variables indicadoras nos permite inferir la siguiente interpretación: $\textrm{exp}(\beta_{i}), i =1,\ldots, p $ $($conocidos como cocientes de riesgo y denotado HR por sus siglas en inglés$)$ corresponden al riesgo proporcional que conlleva que la covariable indicadora afirme que el individuo pertenece a esa categoría, teniendo todo lo demás constante.
Asi:

• $\beta_i < 0$ $(\textrm{exp}(\beta_{i}) < 1)$ indica que el individuo que pertenece a la categoría que indica $X_{i}$ está en menor riesgo a experimentar el evento que alguien que no pertenezca a la misma categoría.
• $\beta_{i} = 0$ $(\textrm{exp}(\beta_i)= 1)$, el riesgo no cambia.
• $\beta_{i} > 0$ $(\textrm{exp}(\beta_i)> 1)$, el riesgo es mayor.

Por lo tanto, los valores de HR nos muestran el riesgo que se tiene con respecto al caso base. Para efectos de este modelo, se considera como caso base a aquellas variables que se definen con un $0$ $($ser hombre y no padecer de ninguna de las afecciones$)$.

Regresión simple de Cox

Con la intención de observar como se traducen las conclusiones que se hicieron durante el análisis exploratorio al modelo de riesgos proporcionales ajustamos de manera preliminar un modelo de Cox a cada covariable sin ajustar para las demás. Los resultados de este proceso están resumidos en la tabla 1. De forma consistente a nuestras observaciones previas resaltan los cocientes de riesgo $($HR$)$ para variables explicativas como Sexo, Neumonía y Condición Renal Crónica. La primera por la disminución marcada de riesgo asociada a no ser del sexo masculino y las otras por el drástico aumento que esta asociado a presentar estas condiciones.

Variable	beta	HR	IC 95 % para HR	waldTest	pValue
SEXO NEUMONÍA DIABETES EPOC ASMA INMUSUPR HIPERTENSIÓN CARDIOVASCULAR OBESIDAD RENAL CRÓNICA TABAQUISMO	-0.4439 2.108 1.122 1.14 -0.3008 0.7236 1.062 0.8696 0.4141 1.34 0.1319	0.6415 8.228 3.071 3.127 0.7402 2.062 2.892 2.386 1.513 3.82 1.141	$($0.6033,0.6822$)$ $($7.697,8.796$)$ $($2.896,3.257$)$ $($2.775,3.524$)$ $($0.6137,0.8928$)$ $($1.755,2.422$)$ $($2.729,3.065$)$ $($2.111,2.696$)$ $($1.419,1.613$)$ $($3.432,4.251$)$ $($1.036,1.257$)$	200.4 3835 1398 350.1 9.89 77.49 1289 194.4 161.4 602.3 7.1	1.719e-45 0 4.823e-306 3.955e-78 0.00166 1.331e-18 3.537e-282 3.525e-44 5.494e-37 5.37e-133 0.007696

Dado que para muestras grandes las pruebas de Wald, Score y cociente de verosimilitud son prácticamente análogas⁶John P Klein and Melvin L Moeschberger. Statistics for Biology and Health. Springer, 2 edition, 2003., nos permitimos utilizar la prueba de Wald para llevar a cabo pruebas sobre la significancia de la $\beta$ asociada a cada modelo. Se puede observar en la tabla 1 cómo los coeficientes con menos significancia corresponden a las covariables cuyo múltiplo en la figura 1 era relativamente cercano a 1. Así mismo resulta curioso caer en cuenta de la similitud de los cocientes de riesgo y la figura 1. Sin embargo, esto no es más que un primer paso y resulta necesario construir un modelo en donde se ajuste para varias variables simultáneamente.

Regresión multivariada de Cox

Al ser un tipo de regresión, el problema de determinar cuáles covariables considerar como variables explicativas surge naturalmente. Optamos por un método iterativo retrospectivo basado en el Criterio de Información de Akaike conocido como el Criterio de Información Bayesiano $($BIC$)$, que penaliza aún más la inclusión de variables adicionales⁷Douglas C Montgomery, Elizabeth A Peck, and G Geoffrey Vining. Introduction to linear regression analysis. John Wiley & Sons, 5 edition, 2012. al modificar el AIC para que refleje el tamaño de la muestra. $($AIC$)$⁸Douglas C Montgomery, Elizabeth A Peck, and G Geoffrey Vining. Introduction to linear regression analysis. John Wiley & Sons, 5 edition, 2012..
El método iterativo de eliminación consiste en eliminar aquella covariable cuyo valor para AIC habiendo ajustado para las otras $p-1$ covariables sea el máximo, considerando que este valor sea menor al mínimo de la iteración pasada. El método termina cuando ya no se logra reducir el criterio más.
Para elegir los modelos de regresión, utilizamos la función step disponible en la documentación de R. Esta función lleva a cabo el procedimiento que acabámos de describir. El resultado se encuentra resumido en la tabla 2.

Variable	beta	HR	Error Estándar	z	p
SEXO HIPERTENSIÓN NEUMONÍA EDAD DIABETES EPOC INMUSUPR OBESIDAD RENAL CRÓNICA TABAQUISMO	-0.335880 0.166807 1.649977 0.035646 0.296953 0.184210 0.237748 0.315657 0.589242 -0.103752	0.714709 1.181526 5.206861 1.036289 1.345753 1.202269 1.268390 1.371160 1.802622 0.901449	0.032011 0.033231 0.035290 0.001075 0.032485 0.063127 0.083609 0.033409 0.056870 0.050501	-10.493 5.020 46.755 33.155 9.141 2.918 2.844 9.448 10.361 -2.054	2e-16 5.18e-07 2e-16 2e-16 2e-16 0.00352 0.00446 2e-16 2e-16 0.03993

En esta instancia, se utiliza una prueba local de Wald para probar la significancia de cada $\beta$. A comparación de los ajustes individuales de la tabla 1, este modelo muestra cocientes de riesgo relativamente cercanos a 1. Destaca el fuerte aumento del riesgo al presentar Neumonia y la disminución de riesgo que presenta el padecer Tabaquismo y ser del sexo femenino. Sin embargo, notamos que la significancia observada para el estadístico de Wald asociado a tabaquismo es relativamente alto con respecto a los otros valores.

Validación de supuestos del modelo

El supuesto clave detrás del ajuste a un modelo de riesgos proporcionales es, en efecto, que los riesgos sean proporcionales en cada covariable. Existen métodos gráficos y estadísticos para llevar a cabo esta validación. Ambas formas presentan sus ventajas y desventajas. Por su parte, los métodos estadísticos de ajuste son susceptibles al tamaño de la muestra⁹John P Klein and Melvin L Moeschberger. Statistics for Biology and Health. Springer, 2 edition, 2003., y para muestras grandes suelen rechazar la hipótesis nula de que el modelo se ajuste bien. Los métodos gráficos de residuos son una forma de apoyarse, pero fallan en el momento que se introducen interacciones entre las covariables o cuando se introducen interacciones. Dado que este modelo no incluye interacciones, consideramos apropiado llevar a cabo un diagnostico gráfico, como recomiendan Klein y Moeschberger¹⁰John P Klein and Melvin L Moeschberger. Statistics for Biology and Health. Springer, 2 edition, 2003..

Los residuos de Cox-Snell presentan una herramienta para verificar el supuesto de ajuste a un modelo de riesgos proporcionales, pues bajo esta hipótesis nula, estos residuos vienen de una distribución exponencial con media unitaria. Por lo tanto, si se grafíca el estimador de Nelson-Aalen para la función de riesgo acumulada utilizando a los residuos como pseudo-tiempos, el estimador debería ajustarse a la identidad a través del origen, con la excepción de la cola derecha, cuando la estimación pierde precisión drásticamente por su construcción¹¹John P Klein and Melvin L Moeschberger. Statistics for Biology and Health. Springer, 2 edition, 2003..
Incluimos el diagnóstico de los residuos Cox-Snell para el modelo 1 en la Figura 5, haciendo hincapié en el hecho de que el cuantil 0.63 de esta distribución corresponde aproximadamente a 1, por lo que el ajuste del modelo no es plenamente satisfactorio. Es decir, es difícil argumentar que la desviación del estimador de la identidad que se comienza a observar a partir de t = 1 es atribuible a estar cerca de la cola derecha.

Ahora bien, otro supuesto que va incluido tácitamente en el modelo 1 es que las covariables son funciones constantes del tiempo. Los residuos de martingalas proporcionan una herramienta valiosa para validar la mejor forma funcional de la variable continua Edad y los residuos de Schoenfeld sirven para verificar si los coeficientes $\beta$ son funciones constantes del tiempo¹²John P Klein and Melvin L Moeschberger. Statistics for Biology and Health. Springer, 2 edition, 2003. o si los cocientes de resigo cambian en algún momento del tiempo.

Se puede apreciar en la Figura 6 que el spline que ajusta a los residuos de martingala se aleja de cero conforme crece la edad, lo cual sugiere que una discretización de la covariable Edad es apropiada.

Refinación del modelo

Una forma de lidiar con covariables que parecen tener dependencia con el tiempo es estratificar¹³Patricia M Grambsch and Terry M Therneau. Modeling survival data: extending the cox model. Statistics for Biology and Health, 2000.. Estratificar consiste en considerar a grupos (estratos) distintos en la variable correspondiente. Una refinacion del modelo consiste en estratificar a la covariable de Edad como una variable categórica de 4 niveles:
$$[0,35], \qquad (35,45], \qquad (45,70] \qquad \text{y}\qquad (70,\infty).$$
Intuitivamente esto hace sentido al considerar las observaciones hechas anteriormente respecto a la tasa de mortalidad observada al estratificar a la población por décadas de vida. Posteriormente, se pueden considerar interacciones entre las nuevas variables y ajustar un nuevo modelo. Sin embargo esto se debe de hacer con cuidado, pues como mencionamos anteriormente, las pruebas gráficas no pueden ser utilizadas para validar supuestos, por lo que uno depende en pruebas estadísticas para llevar a cabo esta validación. Estas refinaciones se llevaron a cabo en el estudio original, cuya extensión y alcance es considerablemente mayor de lo que los autores desean transmitir en este artículo. Empero, las conclusiones apuntan en la misma dirección.

Conclusiones y limitaciones

Uno de los supuestos fundamentales que se hicieron para este trabajo fue asumir que las condiciones médicas con las que se cuenta con información son condiciones pre-existentes. Mientras que es posible que el paciente se haya presentado al hospital con neumonía, por ejemplo, es también posible que esta se haya desarrollado después de que haya sido internado. La posibilidad que un paciente desarrollara la condición después de haber sido internado no se tomó en cuenta.
Adicionalmente, se modeló al COVID-19 como una enfermedad crónica en el sentido de que cualquier persona que no hubiera muerto por complicaciones asociadas a COVID podría morir en algún tiempo futuro, sin contemplar la posibilidad de que un paciente se recuperara del virus. De forma similar a esto, se asumió que dentro de la población positiva a COVID no habían riesgos en competencia como lo podrían ser contraer alguna otra enfermedad al estar internado en el hospital, de la misma forma que no se incorporó una metodología que tomara en cuenta el efecto de la saturación de hospitales en el tiempo de vida de los pacientes.

Se reconoció desde un principio que el hecho de que un paciente estuviera intubado era un factor que artificialmente extiende los días de vida para un paciente que lo necesita. Al haber escasez de estos aparatos y verse la accesibilidad a ellos limitada, se podría dar el caso de que la supervivencia se recorriera a la derecha por este pequeño sesgo. Sin embargo, por la pequeña cantidad de pacientes intubados en relación a los hospitalizados y la gran cantidad de NA’s en este rubro $($arriba del 60%$)$, se decidió llevar a cabo el análisis tomando a estos pacientes en cuenta, sin utilizar al criterio como regresor.
El mismo supuesto de haber ajustado un modelo de riesgos proporcionales puede ser cuestionado. En principio, las pruebas estadísticas para determinar el ajuste del modelo suelen rechazar el supuesto de riesgos proporcionales para muestras suficientemente grandes. Si bien, pruebas gráficas como plots de Arjas, Andersen o procesos de Score proporcionan una herramienta poderosa para complementar este tipo de pruebas, la incorporación de interacciones al modelo para intentar entender mejor las relaciones que la muestra deja entrever se da a cuestas de pérdida de métodos para la validación de los supuestos básicos.
Ahora bien, en todos los modelos que fueron ajustados a los datos, se observó que las mujeres tienen menor riesgo que los hombres y que neumonía es la enfermedad que más se relaciona con el desenlace de muerte. También se le dió sustancia a las observaciones internacionales que muestran que el grupo de mayor en riesgo son las personas mayores.
Finalmente, hallamos crucial el recalcar que el modelo de riesgos proporcionales debe de ser utilizado para identificar a grandes rasgos cuáles factores contribuyen a aumentar el riesgo de morir para un individuo. Por la misma razón, el hecho de que los supuestos de proporcionalidad constantes a través del tiempo no se cumplan perfectamente no es un factor que invalide al modelo, dado que el objetivo en esta instancia es identificar los factores que pueden complicar un diagnostico de COVID-19, no llevar a cabo una predicción sobre tiempo de vida para un paciente diagnosticado. La importancia de la prevención, particularmente para la gente que cae en las categorías de pre-condiciones que han sido discutidas en este trabajo, es esencial para limitar la cantidad de muertes.