Statistical methods

The distribution of covariates was assessed for the entire cohort and across waves. Chi-square test was used to assess independence. Mortality rates and person-time of observation were calculated for the main exposure groups and all covariates of interest. The strength of the association was quantified using incidence rate ratios (IRR), and the statistical significance using 95%CIs and p values. Survival across the different waves was explored using time-to-event analysis and log-rank to test the significance of the difference between the survival curves.

Rates were modelled using Cox regression for the multivariable analysis. Proportional hazard assumption was supported (graphically and by testing for a zero slope in Schoenfeld residuals).

A causal model was built using a stepwise backward approach where (non-forced) pre-defined covariates were retained in the model unless there were problems with multicollinearity. Age and gender were considered a priori confounders (forced variables). Age was fitted using restricted cubic splines, with knots positioned so numbers of events between knots were approximately equally distributed. The full model included age, gender and all variables found to confound the crude association between wave and mortality (non-forced variables). A change in the magnitude ≥ 5.5% was considered an indication of confounding. ICU admission was included in the model as a non-forced variable regardless of the degree of confounding of the main association. Problems with multicollinearity on the main effect in the full model, were resolved using RMSE (Root Mean Square Error) reduction for backward deletion of non-forced variable, with RMSE for the full model used as reference for each step [21].

Following the same methodology, we did a sub-analysis among those requiring ICU admission. Data management and statistical analysis were carried out using R (R Core Team version 3.6.3, Vienna, Austria).