Multistate Modeling of COVID-19 Patients Using a Large Multicentric Prospective Cohort of Critically Ill Patients

The mortality of COVID-19 patients in the intensive care unit (ICU) is influenced by their state at admission. We aimed to model COVID-19 acute respiratory distress syndrome state transitions from ICU admission to day 60 outcome and to evaluate possible prognostic factors. We analyzed a prospective French database that includes critically ill COVID-19 patients. A six-state multistate model was built and 17 transitions were analyzed either using a non-parametric approach or a Cox proportional hazard model. Corticosteroids and IL-antagonists (tocilizumab and anakinra) effects were evaluated using G-computation. We included 382 patients in the analysis: 243 patients were admitted to the ICU with non-invasive ventilation, 116 with invasive mechanical ventilation, and 23 with extracorporeal membrane oxygenation. The predicted 60-day mortality was 25.9% (95% CI: 21.8%–30.0%), 44.7% (95% CI: 48.8%–50.6%), and 59.2% (95% CI: 49.4%–69.0%) for a patient admitted in these three states, respectively. Corticosteroids decreased the risk of being invasively ventilated (hazard ratio (HR) 0.59, 95% CI: 0.39–0.90) and IL-antagonists increased the probability of being successfully extubated (HR 1.8, 95% CI: 1.02–3.17). Antiviral drugs did not impact any transition. In conclusion, we observed that the day-60 outcome in COVID-19 patients is highly dependent on the first ventilation state upon ICU admission. Moreover, we illustrated that corticosteroid and IL-antagonists may influence the intubation duration.


Supplementary Material of "Multistate modeling of COVID-19
patients using a large multicentric prospective cohort of critically ill patients"

Statistical methods
Multistate models allow for extending the standard survival model to more than two states and, consequently, more than one transition [1][2]. A multistate model describes the individual path across states in continuous time. We considered 6 states: 1. Discharge alive from hospital; 2. Discharge alive from ICU; 3. ICU non-invasive, defined as: in ICU without invasive mechanical ventilation (i.e., with Optiflow or Continuous Positive Airway Pressure (CPAP)); 4. ICU invasive, defined as: in ICU with invasive mechanical ventilation (barometric and positive endexpiratory pressure PEEP ≤ 10 or volumetric and PEEP > 10); 5. ECMO: in ICU with extracorporeal membrane oxygenation; 6. Death.
Patients can start in state 3, 4 or 5. State 1 and 6 are called absorbing state since once the patient has entered one of them, s/he won't move anymore. Moreover, ICU patients are continuously-observed, that is, the state is known at each day; specifically, each patient is associated with the worst state s/he encountered during the day.
Denote by ( ) the state occupied at time and by ( ) the transition intensity that expresses the instantaneous risk of a transition from state into state ℎ at time . It can be expressed as This definition implies that the probability of moving to a future state depends only on the present state and not on the history, therefore multistate model is called Markovian. The cumulative transition hazard from state into state ℎ is computed as ( ) = ( ) . Let ( ) be a 6 x 6 matrix, with elements ( ), Equation (1) can be estimated using the Nelson-Aalen estimator for the cumulative intensities ( ), that is .
Cumulative transition hazards and transition probabilities can be also estimated conditioning on covariates .
In this case, ( ; ) is now estimated for given value of the covariate vector = using a Cox proportional hazard model with Breslow method for handling ties [3].
We first used a non-parametrical approach, using the Nelson-Aalen estimator for the cumulative intensities.
We computed the prediction of state occupancy, Regarding the covariate analysis, a semi-parametric modeling approach via Cox proportional hazard regression, Breslow method for handling ties and robust variance was used. For each ordinal variable, each category was coded as a dummy variable that takes the value of 1 if the variable is at least higher to the lower cut-off of the category, and 0 otherwise. This coding was chosen to better interpret the results after variable selection technique, since more categories can collapse into a single one.
Univariable analysis was first performed and covariates associated with a p-value lower or equal to 0.2 were retained for the multivariate analysis. Due to the small sample size, no interaction was tested. Therefore, we assumed an additive effect for the drugs. Covariates were added only on transitions with more than 10 events and when all covariate categories can be represented. Then, the final model was achieved using a stepwise backward-forward selection for the multivariable analysis using the BIC criterion. The proportional hazards assumption was tested using the scaled Schoenfeld residuals. Due to the possible computation approximation instabilities, the estimated cumulative hazard function conditional to specific a covariate set = were linearly interpolated in order to have values in a denser time space before using Eq. (1).
The effect of corticosteroids and Tocilizumab/Anakinra in the ICU population was tested using a Gcomputation approach [5]. Let's denote the first two components of the vector as the variables and that represent the presence or absence of corticosteroids and IL-antagonists, respectively. For sake of simplicity, we do not write the rest of the vector and we focus only on these two patient's covariates that will be changed in the G-computation approach. The average ICU population effect can be estimated as where , and , denotes the presence ( , = 1, , = 1) or absence ( , = 0, , = 0) of corticosteroids and Tocilizumab/Anakinra, respectively, at admission.
In order to compute confidence intervals for probabilities of state occupancy , (0, ), a probabilistic sensitivity analysis (PSA) was performed via an asymptotic Monte Carlo approximation [6]. Maximum likelihood estimates, by the Cox proportional hazard model, were sample from an asymptotic multivariate normal distribution, with mean equals to the estimated parameters and variance-covariance matrix given by the estimation process. Hundred Monte Carlo runs were performed, and confidence intervals were obtained using 0.025 and 0.975 percentiles (we also checked the impact on the results of increasing the Monte Carlo runs up to 500 and we obtain negligible differences, therefore we opted for 100 runs).
The mean sojourn time at each state is computed as the integral of the probability of being in the state, , (0, ), between zero, the starting point, and the selected final day. Confidence interval for the mean sojourn are then computed using the PSA results runs. To approximate the integral, we used the approximation implemented in the "ELOS" function of mstate R package.

Missing values
Two missing patient states were imputed ad hoc looking at the patient's multistate path. When possible, missing values of the number of days from first symptom to ICU was assumed equal to the number of days from first diagnosis to ICU. Other missing covariate variables were imputed to the median of the ICU population, except for BMI where the median according to the sex was used.

Univariable analyses details
The following variables could not be tested because all their modalities were not represented in the transition: Age >50, and SOFA >8 could not be tested in transition from ICU non-invasive to death; sex and age >70 could not be tested in transition from ICU invasive to ECMO; age > 70, hydroxychloroquine could not be tested from ECMO to ICU invasive; age >70, SOFA >3 could not be tested from ECMO to death . Missing data are shown at the end of each variable. On the right, the discretized variable version used for the semi-parametric model and stepwise variable selection. Null missing data were not listed. BMI: Body mass index; SAPS: Simplified Acute Physiology Score; ICU: intensive care unit; SOFA: sequential organ failure assessment score; CRP: C-reactive protein serum. Figure S1 Stacked plot of predicted probabilities of state occupancy resulting from Gcomputation.
Legend. On the top left, the results when corticosteroids and tocilizumab/anakinra were administered at the admission; on the top right, when corticosteroids were administered to patients without tocilizumab/anakinra; on the bottom left when tocilizumab/anakinra were administered to patients without corticosteroid; on the bottom right when none of these treatments was administered to the ICU population. ICU=intensive care unit; ICU non-invasive: in ICU without mechanical ventilation; ICU invasive: in ICU with mechanical ventilation; ECMO=extracorporeal membrane oxygenation. Figure S2: Probability of state occupancy plots