The funnel plot methodology has been applied to five case studies, corresponding to different stages of the COVID-19 pandemic, chosen in view of their relevance with respect to the spread of VoC's or flaws of the diagnostic infrastructure. Two case studies refer to England (initial spread of Omicron and large failure of a diagnostic lab), and the other three to Italy (initial spread of Omicron), India (first emergence of the Delta variant), and South Africa (first emergence of the Omicron variant). In addition, the nine English regions are monitored over a 18-month period from December 2020 to June 2022.
In all cases, we focus on four dates of interest. The first analyzed date corresponds to a situation of statistical homogeneity: when variants are uniformly spread in the country and contact rates do not vary much across regions, differences between estimated \({R}_{t}\)'s are due to natural variability alone and the regional \({R}_{t}\)’s are expected to lie within the funnel, centered around the national \({R}_{t}\) (see Methods). The second and the third dates refer to the disruption of the natural variability: when a new VoC starts spreading, it first colonizes a few territories, whose behavior becomes abnormal with respect to the national one. This is highlighted by the fact that the corresponding \({R}_{t}\)’s cross the funnel limits. The last date corresponds to a new homogeneity, typically established around a higher \({R}_{t}\): the VoC is now uniformly spread in the country and the condition of natural variability has been restored. Finally, in order to have a snapshot of the whole period of analysis, the standardized \({R}_{t}\)'s with ± 3.09 sigma are plotted on a Bonferroni control chart (see Methods for details). Due to its statistical background, the scope of the new control method is not restricted to VoC monitoring but can detect other kinds of anomalies, e.g., those related to testing availability or malfunctioning buffer factories, an example being discussed in the Immensa case study.
Omicron spread in Italy.
We first apply the funnel plot methodology to the Italian regional data in the period from 4 December to 3 January 2022, based on epidemiological indicators released daily by the Civil Protection Department, which provides 21 regional time series (for 19 regions and the 2 autonomous provinces of Trento and Bolzano). The Delta variant was dominant in Italy until December 2021, when the Omicron variant started to spread across the country.
The results are summarized in Fig. 1. In the Panels a-d, the estimates of Italian regional \({R}_{t}\)'s are plotted against the infectious cases on four selected dates. On 7 December 2021 (see panel a), differences between estimated \({R}_{t}\)'s were due to natural variability alone and the 21 points lay within the funnel limits. On 22 December 2021, Lombardy (dark red) crossed the alarm limit (see panel b) and on 24 December 2021 (see panel c) it was definitely outside the upper alarm limit. In fact, as confirmed by a survey by the Italian National Institute of Health published on 31 December 2021 [26], Lombardy was the first Italian region to be colonized by the Omicron variant. As other regions became increasingly colonized by the Omicron variant, their \({R}_{t}\)'s rose as well and, by 2 January 2022, Lombardy was absorbed again within a funnel, now with a higher mean with respect to the early December’s one (see panel d).
The trend can be monitored by plotting the standardized \({R}_{t}\)'s on a Bonferroni control chart with ± 3.09 sigma limits (see panel e), where the arrival of the Omicron variant in Lombardy in mid-December is clearly detectable.
England statistical monitoring for a year and a half.
The figure displays the Bonferroni control chart of normalized \({R}_{t}\)'s of the nine English regions during 18 months, from the end of November 2020 to the beginning of June 2022. Irrespective of the current national Rt, under natural variability conditions all the normalized curves are expected to lie within the limits. Points outside the limits highlight a disruption of the statistical homogeneity across regions that is worth being investigated in order to retrieve the root cause of the anomaly. In Fig. 2, seven major events are visible, labelled from A to G. In the figure caption, plausible explanations are conjectured: they refer to the emergence or the arrival of the VoCs (Alpha [46, 47], Delta [48], Omicron [31, 45] and Omicron sub-variants [42]), the malfunction of swab factories (also depicted in Fig. 4) [14, 15], some incidents of violation of lockdown restrictions [38, 39, 40], and changes in the testing policies [41].
Delta emergence in India.
We applied our methodology to epidemic data from India in the period 13 February − 5 March 2021, when the Delta variant emerged and started spreading from the state of Maharashtra. In Fig. 3, panels a-d, funnel plots on four selected days are shown, with colour-coded circles representing the \({R}_{t}\)’s of the 36 Indian states. While on February 13 all circles fell within the funnel, on February 16 the state of Maharashtra (dark red) exceeded the alert threshold (in correspondence of the initial spread of the Delta variant), further departing from the mean on February 22. Lastly, on 4 March 2021, the \({R}_{t}\)’s of all regions but Kerala (orange) shaped a new funnel with a higher mean, which again incorporated Maharashtra. The peculiar dropping of Kerala’s \({R}_{t}\) below the lower alert threshold, despite the very high number of infectious cases, might be explained by the overlapping of Alpha and Delta variants during the same period, resulting in a lower \({R}_{t}\) than in the areas predominantly hit by the Delta variant.
In the Bonferroni control chart (see panel e), the rise of the Delta variant in Maharashtra is clearly visible since mid-February 2021. One month later, on March 17, it was made public that a 10-lab research consortium had alerted the Union Health Ministry about a new variant spreading in Maharashtra [27] and a week later a press release was released on the new VoC [28]. This case study suggests that the use of statistical control methods would have anticipated the detection of the variant.
Omicron emergence in South Africa.
From 7 November to 4 December 2021, the Omicron variant colonized South Africa, starting with the province of Gauteng. In Fig. 3 (see panels f-i) four funnel plots are shown, with colour-coded circles representing the \({R}_{t}\)’s of the South African provinces. Until the very beginning of November 2021, the Delta variant prevailed and the differences in\({ R}_{t}\) levels across provinces were merely a result of natural fluctuations (see panel f). By mid-November the Gauteng province crossed the upper alert threshold (see panel g) and then further diverged (see panel h). This is precisely the timing when the Omicron variant was first identified, as declared by the WHO [29], and became a threat [30]. By 3 December 2021, Gauteng was reabsorbed within the funnel, now with a definitely higher mean, following the spread of Omicron in the other provinces and the consequent rise of their \({R}_{t}\)’s (see panel i). On the Bonferroni control chart with ± 3.09 sigma limits (see panel j), the out-of-control trajectory of the Gauteng province (red) is plainly evident.
Omicron spread in England.
From 4 December 2021 to 1 January 2022, the Omicron variant massively spread in England. The four funnel plots show colour-coded circles corresponding to the \({R}_{t}\)’s of the English regions, see panels k-n of Fig. 3. On December 4, all the regions were within the alarm limits (panel k). By 10 December 2021, the London region was out of the funnel (panel l), further diverging from the upper limit on 15 December (panel m). This suggests that Omicron was more prevalent in London than in the rest of England and, indeed, on 13 December 2021 it was reported that 20% of the cases in England and over 44% in London were due to Omicron [31]. As the other regions were colonized, the distribution of their \({R}_{t}\)’s moved upward and, on 23 December 2021, the London region was again inside the funnel (panel n). An early detection would have been allowed by the Bonferroni control chart, where London first crossed the alarm limit in early December (panel o).
Immensa scandal in England.
Our last case study concerns England in the period from 27 August to 25 September 2021. In Panels a-d of Fig. 4, four funnel plots in selected dates are displayed, with colour-coded circles corresponding to the \({R}_{t}\)’s of the English regions. On 5 September 2021, all English regions were within the funnel (panel a). By 9 September 2021, the South West (red) had moved below the lower alarm limit (panel b) and remained below the lower limit for about two weeks (panel e). The timing of this swing coincides with the period during which the Immensa lab in Wolverhampton gave some 43,000 incorrect negative tests relative to South West and West Midlands territories [14, 15]. While the suspension of lab operations came in mid-October, the Bonferroni control chart indicated an out-of-control condition already in early September and would have allowed a much earlier detection of the anomaly.
Trajectory plots and funnel movies
In order to gain a better understanding of the time evolution of regional \({R}_{t}\)’s, it is useful to complement funnel plots with the time dimension. This can be done in two ways: the first approach displays the trajectories of the regional \({R}_{t}\)’s and an example relative to South Africa is provided in Fig. 5. Alternatively, the funnel plot can be animated yielding a “funnel movie”, where both the trajectories and the shapes of the funnels are iteratively updated (see Supplementary {2} for the description and {3},{4},{5},{6} for the movies). As observed in [51], when \({R}_{t}\) is plotted against the infective subjects, the trajectories exhibit a peculiar clock-wise spiral-shaped pattern. An analogous behavior was observed also in [50, Fig. 2], where trajectories were plotted in the plane of infected subjects against “cooperators”, a variable connected with \({R}_{t}\).
Monitoring regional \({R}_{t}\)’s
In addition to the Bonferroni control charts, a further monitoring and visualization tool is obtained by plotting the regional \({R}_{t}\)’s together with the national one and control limits defining the in-control band for the specific region. A distinct plot is needed for each region because the width of the in-control band depends on the number of the infectious subjects in that region. On the one hand, this means that, differently from the Bonferroni control chart previously introduced, as many plots are needed as the number of regions. On the other hand, this visualization provides a direct display of the regional \({R}_{t}\) and may therefore be more easily understandable. Examples of these plots for Italy, England, India and South Africa are given in the Supplementary {1}.