Likely community transmission of COVID-19 infections between neighboring, persistent hotspots in Ontario, Canada

Preparation of Ontario COVID-19 case data from ICES

Access to anonymized Ontario COVID-19 test results through the ICES Research and Analytic Environment (RAE) was obtained through Applied Health Research Questions (AHRQ) project #P0950.096 approved by the Ontario Ministry of Health and Long-Term Care (MOHLTC). ICES collects health care data for analysis and evaluation of the health system. Ontario-wide COVID-19 data tables maintained in ICES were pre-processed to ensure that only the chronologically first laboratory-confirmed viral RNA polymerase chain reaction (PCR) SARS-CoV-2 test of a specific individual was counted (while repeat infections occur, these are assumed to be significantly less frequent). PCs associated with each positive case were cross-referenced to their respective FSA (derived from the 2016 Canadian Census boundary file from Statistics Canada) to identify and filter out invalid entries. The PCs provided were sourced from the OLIS (Ontario Laboratories Information System) and the RPDB (Registered Persons) databases. The OLIS postal code was more likely to relate to a current address and was therefore prioritized. However, it is also the most error prone. Reasons for rejection include blank entries, invalid FSAs (e.g., out of province FSA), or intentional masking. The RPDB PC was used if the OLIS PC was rejected. Entries where both PCs were invalid were rare (<0.1% of all unique positive cases). Non-residents or those ineligible for the Ontario Health Insurance Plan (OHIP) are also not included. Cases were separately aggregated by FSAs and PCs, then converted into tables structured as rows (FSA or PC count) by columns (date) as input to ArcGIS Desktop. Cases were based on the date of sampling, rather than when test results were reported, to more accurately reflect the inception of infections.

As of June 2021, three separate and distinct waves of COVID-19 cases have been recognized in Ontario (and elsewhere). We delineated the inception and termination dates of each wave from the rates of change in infection rates (cases/day²) using the second derivative of daily case counts. Stochasticity in day-to-day variation due to testing and case reporting caused fluctuations in changes of daily counts, even when averaging cases over multiple consecutive days. Consecutive dates without fluctuations in counts (between -10 and +10 cases/day²) were rare. Therefore, dates were selected where the second derivative remained between -10 and +10 cases/day² immediately preceding rapid increases in cases that defined inception of a wave, or after infections ceased decreasing in counts, defining termination of a wave. The first wave (1) was found to span from March 26^th to June 16^th, 2020, second wave (2) began on September 25^th, 2020 and continued until February 14^th, 2021, and the third wave (3) covered March 17^th to June 23^rd, 2021.

Of the geostatistical methods tested, PC-level counts of infected individuals were the most geospatially resolved granular COVID-19 case data. Many PCs comprise of small geographic areas with low population densities. This made them more prone to the effects of day-to-day stochasticity in patient ascertainment and laboratory reporting, especially when case counts were low. To address this, PC-level case data was also aggregated over multi-day sliding windows, effectively smoothing out this source of variation. Although 3 and 7-day windows were initially analyzed, longer windows were not as useful for determining the inception of hotspots, a goal of this analysis, and were eventually discounted. 3-day windows have been used in other geostatistical studies of COVID-19.^13,14

Access to, pre-processing and geostatistical analysis of COVID-19 daily cases and associated metadata required the RAE computer system at ICES. To ensure patient privacy, ICES stipulated that case counts could not be exported from the ICES system unless FSAs and/or PCs exhibiting between 1 and 5 daily COVID-19 cases were masked. We modified counts of FSAs meeting these criteria and assessed their impact on the geostatistical analysis. Gi* analysis of ground truth case data was compared with the analysis of two modified FSA-level case datasets (where all instances of between 1 and 5 case counts were either converted to a single count or randomized between 1 and 5). More than one third of the locations that were deemed significant hotspots by Gi* from actual case counts (from March 1^st to September 24^th, 2020) were not flagged when locations were masked and assigned a count of 1; greater than one half were missed when case counts were randomized (Extended data,²⁹ Section 1 - Table S1). Discordant hotspots identified from actual vs. masked case counts were reduced in FSAs with high-case counts (≥10) but were nevertheless unacceptably frequent (14% or 28% discordant, respectively, by masking low case counts as 1 or by randomization). We concluded that errors introduced by masking count values precluded analysis of masked data and required all work to be carried out within the RAE computing environment. All geostatistical analyses presented in this study were performed against unmasked, true case counts.

FSA and PC boundary files

Shapefiles which store geospatial data [points, lines or polygons] with related attribute information and depict the 2016 census geographic boundaries of FSAs across Canada were obtained from Statistics Canada. PC boundary shapefiles created by DMTI (“DMTI_2019_CMPCS_LocalDeliveryUnitsRegion.shp”) were accessed through the Scholar's Geoportal at the University of Western Ontario. These files were used to compute latitude/longitude coordinates of the centroid for each FSA and PC (with the ArcGIS “Calculate Geometry” function), to validate PCs provided for each individual in the COVID-19 test dataset, and to convert spatial interpolation maps (such as those generated by kriging) into machine-readable text files (using ArcGIS “Intersection” function between the spatial map and with each boundary file). The ArcGIS “Split by Attribute” function was used to limit PC boundary files to the province of Ontario. This was necessary, both because data from other provinces were not available, and because of system limitations affecting processing and memory requirements for large geographic regions.

Minor discrepancies between the FSA and PC boundary files were identified which had little or no impact on our geostatistical analysis. Nearly 9% of all PCs in Ontario overlapped multiple FSAs. Most of these PCs were present in multiple locations and occurred in rural areas. Approximately 2.7% of PCs intersected with an incorrect FSA, which was likely the result of asynchronous creation of the corresponding shapefiles. These locations rarely contained many COVID-19 cases as only 11 PCs had >5 cases over any 3-day span in 2020, including: N9H 0E3 (Windsor, ON) which overlaps with the FSA N9G; N8N 0B3 (Tecumseh, ON) which overlaps the bordering FSA N8R, and M3H 5V9 (Toronto, ON) which intersected with neighboring FSA M3J. These discrepancies might impact the interpolation of cases in FSAs based on kriging analysis. Since the discrepant PCs were typically found on the border of FSAs, these conflicts had minimal effect on kriging accuracy.

Geostatistical analysis

COVID-19 cases from March 26^th, 2020 to June 25^th, 2021 at the FSA (Ontario-wide) and PC-level (for 6 populated areas: Hamilton, Kitchener/Waterloo, London Ottawa, Toronto, and Windsor, including Essex county) were evaluated using ArcGIS Desktop 10.7 [ESRI] with the Geostatistical Analyst extension. Daily cases compiled by FSA were evaluated with Cluster and Outlier analysis, Gi* ‘hotspot’ analysis, and with space-time integrated Gi* analysis. Cluster and Outlier analysis computes a Local Moran’s I index score for each feature. A positive score indicates that a feature is clustered with other similar high-case count regions, while a negative index indicates that a feature is an outlier. The z-score and p-value indicates the statistical significance of the index value. A false discovery rate (FDR) correction was applied to this and Gi* analysis, which changes the threshold of significance of p-value ≤0.05 to a reduced value within the 95% confidence interval [C.I.], thereby accounting for multiple testing and spatial dependency. Z-scores are also determined in Gi* analysis, with higher positive scores indicating greater clustering of cases in a region. Space-time Gi* analysis was based on a single day interval but incorporated temporal trends based on the days preceding and subsequent to each date. We set a minimum threshold of 10 confirmed cases for an FSA to be considered significant by these three approaches, as this threshold was found to provide reliable and conservative detection of known hotspots using the Gi* test. Over two thirds of FSAs with ≥10 cases between March 1^st and Sept. 24^th, 2020 were found to be significant by Gi* [99% C.I.] using both a fixed and inverse distance metric (Extended data,²⁹ Section 1 - Figure S1). Spatial relationships were evaluated to determine which distance metric best reflected the relationship between the test location and its surrounding FSAs by fixed distance (all FSAs within a given threshold distance [10km] are equally weighted with the test FSA), inverse distance (FSAs closer to the test FSA are weighted more highly), and by K Nearest Neighbors (KNN; comparing counts in the test FSA and its closest K neighbors).

PC-level case data (binned over 3-day sliding windows) for each municipality were evaluated by Spatial Autocorrelation (Global Moran’s I), Cluster and Outlier analysis, and EBK. Spatial Autocorrelation used FDR correction to account for multiple testing; it was not applied for PC-level Cluster and Outlier analysis due a limitation of the method in high-case count regions where all locations were deemed non-significant. The KNN distance metric was used to define the spatial relationship between features (with K = 3) for all PCs, including those with zero cases. PCs with ≥6 cases were reported if deemed significant by Cluster and Outlier analysis. This is consistent with case minimums used in other COVID-19 geostatistical studies.^13,14 Kriging generates a topographic contour map of the number of COVID-19 cases interpolated over the evaluated region, as indicated below.

COVID-19 testing data compiled by ICES contained metadata on tested individuals, such as gender, age, and relevant predicate health conditions. To identify geospatially significant cases in these subgroups, case data was stratified by gender (Male: 241,393; Female: 247,152; those without gender information were ignored), by age (≥65 years old: 63,746; <65 years: 424,854), and by presence or absence of at least one chronic health condition (187,242; includes chronic asthma, congestive heart failure, chronic obstructive pulmonary disease, dementia, hypertension and/or diabetes). 301,358 patients were without one of these conditions. Individuals with these conditions are grouped as an “at-risk population”, consistent with provincial terminology. These case subsets were then evaluated with Gi* (FSA-level), Spatial Autocorrelation, Cluster and Outlier analysis and kriging (PC-level).

Results from PC-level Cluster and Outlier analysis were further analyzed with a program script to identify significantly, highly clustered PCs over consecutive periods (i.e., termed high-case cluster “streaks”; Local_Morans_Analysis.Recurrent_Clustered_PC_Identifier.pl; see Data Availability section²⁹). We also identified which of these PC streaks occurred in close proximity (Local_Morans_Analysis.Clustered_Streak_Pairing_Program.pl) based on the neighbors of each postal code from PC centroid data (neighbors identified by Ontario_City_Closest_Postal_Code_Identification.pl; see FSA and PC boundary files subsection). Directional network graphs which visualize the interconnectivity of these paired high-case cluster streaks were created in MATLAB using built-in digraph and plot functions, where each edge (the connection between two PCs) represents a single pair of PCs with nearly concomitant or overlapping streaks.

Software automation was required to computationally evaluate daily and windowed case data over 16 months at the geographic resolution of FSAs and PCs across multiple cities with different geostatistical tests. Using the ArcPy package (ESRI), which enables Python programs to access ArcGIS’s geoprocessing tools and extensions, program scripts were developed to perform geostatistical analyses across Ontario. These scripts set conditional variables which define the spatial relationship used (fixed, inverse or KNN), the threshold distance (or for KNN, the number of neighbors of the target feature), the distance method (how distances were calculated; set to ‘Euclidean’), row standardization (set to ‘none’), and whether FDR correction was used (applied to FSA-level Gi* and Cluster and Outlier analysis, and PC-level Spatial Autocorrelation analysis). ArcPy scripts read in case data (in tabular format), iterated sequentially through each date skipping those without positive cases, applied the associated geostatistical tool, and converted the output into text files. Tabular output was subsequently filtered by removing PCs without cases to avoid producing outputs of prohibitive size.

These processes were then integrated and streamlined in a software package named the Geostatistical Epidemiology Toolbox (see Data Availability section²⁹ ), which is accessed through the ArcMap graphical user interface. This resource provides the same geostatistical operations and map imaging developed for this project in a simplified, user-friendly environment. This toolbox bundles manual pre-processing, analysis, and post-analytical steps into a set of linked, uninterrupted workflows. While the toolbox was designed to evaluate COVID-19 case data in Ontario, the software can be easily modified to accommodate other Canadian provinces and, with additional effort, other countries (see Geostatistical Epidemiology Toolbox user manual²⁹). To enable kriging functionality, the Advanced Geostatistical Analyst (ESRI) extension and the PC and FSA boundary shapefiles are required. These shapefiles are also required for PC-level Local Moran’s I and Global Moran’s I analyses. This package is not bundled with any of the program scripts developed to evaluate geostatistical output (besides data imaging), which are also provided in the accompanying Zenodo archive.²⁹

Empirical Bayesian Kriging analysis

Kriging estimates case counts in unsampled regions by interpolation between counts at known locations. These values are spatially correlated with degree of separation from centroids of PCs or FSA with known cases. A topographic contour map is generated which illustrates the interpolated case levels within the analyzed region. The spatial distribution of the disease outbreak was determined for Hamilton, Kitchener/Waterloo, London, Ottawa, Toronto, and Windsor/Essex county by kriging analysis of PC-level COVID-19 case data over consecutive, overlapping 3-day windows. Kriging was also performed for other municipalities bordering Toronto (Ajax, Brampton, Markham, Mississauga, Pickering, Richmond Hill and Vaughan; see results in Extended data,²⁹ Section 3). Only PCs with ≥1 case were included in this analysis, as locations without confirmed cases greatly outnumbered PCs with cases (≥95% of PCs on any given date evaluated) leading to the severe depression of kriging signals. Several types of kriging (Ordinary, Universal, Simple and EBK) were evaluated for accurate detection of high case clusters in municipalities with high and moderate population densities (Toronto and London, respectively). Contour maps generated by EBK best matched regions of known cases used for analysis. Ordinary and Universal kriging yielded similar results for locations with high case counts (≥50) but lacked the sensitivity of EBK in regions with moderate cases. Simple kriging failed to define most contours correctly, regardless of case count. Therefore, EBK was chosen for all kriging analyses presented in this study.

Semivariograms describe how the distances separating COVID-19 positive individuals and the semivariance of counts at these locations are correlated. The best model that fits this relationship using empirically defined distance data was used to predict COVID-19 count levels at unsampled locations between the centroids of FSAs or PCs. Semivariogram models assessed included Power, Linear, Thin Plate Spline (TPS), Whittle, Exponential, and K-Bessel. Kriging parameters for these models were considered, such as sector type, search radius, neighborhood type and the number of neighbors considered for each location. Through empirical testing, the Power semivariogram model (based on a parabolic relationship to distance between locations) was chosen based on interpolated values that best represented actual case counts in low/moderately dense populations (Extended data,²⁹ Section 1 - Figure S2A). Kriging accuracy improved by inferring interpolated values from 3 neighboring locations with known counts. However, the Power semivariogram to failed to detect some hotspots in Toronto, which exhibited the highest case counts and overall population density. The TPS semivariogram model, a non-rigid representation that passes exactly through points while minimizing surface curvature, performed best in this city (Extended data,²⁹ Section 1 - Figure S2B). We noticed that the separation between Power-based kriging contours was generally larger, in contrast with TPS-based contours which were tightly packed, more bifurcated and frequent. TPS, however, tended to generate elevated estimated counts that did not correlate with known hotspots. Using 5 neighboring locations for interpolation minimized these artifacts without affecting the contour map. However, with 10 neighbors, interpolation accuracy was affected according to a chi-square test (e.g., p < 0.05 for the 103 days in Wave 2). The Power semivariogram did not significantly differ from the actual counts (p-value>0.05). However, interpolation with the TPS model in Toronto was consistently more similar to the actual case counts than the Power model, based on Chi-square test for all hotspots within Waves 1 and 2 (resulting in higher p-values for 85% and 78% over all dates within each wave). We performed kriging in Hamilton, Kitchener/Waterloo, London, Ottawa and Windsor/Essex with the Power model (3 neighbors), and in Toronto with the TPS model (5 neighbors). All other EBK parameters used were defaulted to those provided by the ArcGIS Advanced Geostatistics toolbox (Sector Type: 1 sector, Neighborhood Type: Standard, Search Radius: 1).

We have previously demonstrated that circumscribing the boundaries of analyzed regions by specifying these locations with zero counts improves kriging accuracy, especially at locations proximate to boundaries.³⁰ This approach was used in the present study of interpolated COVID-19 cases for each municipality in Ontario. Municipal borders were derived by combining all PCs associated with the region of interest. The optimal densities and distances of the zero-count buffer circumscribing municipalities by programmatically generating boundaries. Distances between the boundary to the city border were varied from 0.1, 0.5, 2, and 5 km, and densities between locations of 0.05, 0.25, 0.5 and 1 km. Kriging analysis was iterated for each of these defined zero-count borders, and the kriging layers generated were compared to the original kriging results. The location of each kriging contour properly coincided with regions of high case counts, regardless of the shape and density of the zero-count boundary locations. Raising the boundary density consistently increased the COVID-19 case value interpolated of a high-case region while simultaneously reducing the total area of the kriging contours around the region, regardless of semivariogram model used. This effect was less prominent in regions with higher density of cases (i.e., Toronto), and in regions further from the city boundaries (as the kriging contours exhibited consistent shape and magnitude >6 km away from the city boundary). Considering these observations, we chose to use a moderate value for each parameter (0.5 km density and 1km distance from the city boundary). The kriging tool in the Geostatistical Epidemiology Toolbox, however, automatically generates these boundaries 1 km from the selected border with 250 m of distance between each location.

Kriging analysis generates a series of shapefiles containing contour maps that define the interpolated level of COVID-19 cases of a region. An ArcPy script intersected map contours derived from kriging and either FSA or PC boundary files, which associated each kriging-derived case count with these locations. The results of these intersections were displayed as distinct contour distributions (or kriging ‘symbologies’): one for instances with a wider range of inferred case count levels (0-1, 1-5, 5-10, 10-15 … 35-40, 40-50, >50 cases), and another for low densities of inferred cases, which was discriminated by single counts (0-1, 1-2 … 9-10, >10 cases). In compliance with ICES’s reporting criteria, we defined any interpolated value >5 cases (either ≥5-6 or ≥5-10 contour, depending on the symbology type used) as a region of interest, and display PC or FSA designations for only these locations. The intersected area (m²) between contours and the corresponding FSA or PC was determined to compute the size of localized hotspots (e.g., the total area of the FSA or PC considered a region of interest by kriging).

Wave 1

Kriging analysis of COVID-19 cases in this wave identified localized hotspots in Hamilton, Kitchener/Waterloo, London, Ottawa, Toronto, and Windsor/Essex (all localized hotspots are identified in Extended data,²⁹ Section 2). High disease burden FSAs (those flagged by Gi* analysis) were often found to include localized hotspots (>50% of flagged FSAs in cities evaluated; Extended data,²⁹ Section 1 - Table S4). Localized hotspots (regions within contours exceeding interpolated count >5 cases) were most frequently identified within high disease burden FSAs in Ottawa during the first COVID-19 wave (84%; Extended data,²⁹ Section 1 - Table S4). However, kriging also identified 34 additional localized hotspots within FSAs on dates within this period, suggesting that COVID-19 events of interest occur that are missed by Gi* analysis. Due to its high density, kriging analysis of cases in Toronto using the Power semivariogram (used in analysis of all other cities) was found to have poor overlap with Gi* analysis (only 31.6% of Gi*-significant FSAs contained a localized hotspot). The overlap between kriging and Gi* test results was improved using the TPS semivariogram (increased to 60.9%), and this was therefore used for all subsequent kriging analyses of Toronto. PC clustering was uncommon during this wave (<15% of FSAs flagged by Gi* and/or kriging were clustered). The exception was London where 27% of localized hotspots occurred in clustered FSAs [N = 11], and the majority of cases were often localized to a single PC in other cities that were analyzed (Extended data,²⁹ Section 1 - Table S4).

COVID-19 outbreaks frequently affected residents of long-term care residences during this wave. Geostatistical analyses confirmed the ability to identify these outbreaks. For example, an outbreak in a nursing home in Toronto in early April 2020 (M8V; Lakeshore Lodge) was identified by kriging, while the FSA in which it occurred was also significant by Gi* analysis ( Figure 5A). A different outbreak in Toronto in April 2020 (M2K; Extendicare Bayview nursing home) was also detected by kriging analysis, but the FSA in which it was found was not a significant hotspot by Gi* ( Figure 5B; it was significant only using a fixed distance metric of ≥20 km). Similarly, the FSA K1J (Ottawa, ON) was only significant by Gi* using an inverse distance weighted test, despite a mid-April 2020 outbreak in the Rothwell Heights assisted living community which was identified by kriging (Figure 5C). Localized hotspots within the FSA M9P were observed as early as March 31^st, 2020 however the FSA was not flagged by Gi* analysis until April 6^th, 2020 (Extended data,²⁹ Section 1 - Figure S3A). The localized hotspot covers the Village of Humber Heights Retirement Home in Toronto, which was reported to have a COVID-19 outbreak on April 12^th, 13 days after the localized hotspot was first detected by kriging.³¹ In each of these instances, the interpolation of cases by kriging were similar to actual case counts in these FSAs. This was consistent with the concentration of cases within a single PC (a frequent occurrence in the first COVID-19 wave; Extended data,²⁹ Section 1 - Table S4). Hotspots where the case total exceeds the interpolated value is consistent with cases distributed among multiple PCs across an FSA. For example, localized hotspots in FSA M9V were identified sporadically, and the interpolated value does not always match total cases (Extended data,²⁹ Section 1 - Figure S3B). The localized hot spots cover the PC M9V 5B5 starting from April 4^th. Outbreaks in the long-term care home within this PC (Humber Valley Terrance long-term care home in Toronto) were reported on April 16^th, 2020.³² This FSA was not flagged by Gi* analysis until April 19^th (Extended data,²⁹ Section 1 - Figure S3B).

Figure 5. Geostatistical analysis of COVID-19 outbreaks during wave 1.

FSAs exhibiting high-disease burden by Gi* analysis were identified with localized hotspots by kriging or as having significant case clusters by Spatial Autocorrelation. Panels indicate cases in all PCs within FSAs (A) M8V (B) M2K, and (C) K1J, which were amalgamated over 3-day sliding windows and masked (red line indicates actual case counts; privacy consideration require dates with fewer than 5 cases to be masked as zero). Kriging interpolated the number of cases over these FSAs (blue line), which closely approximates true case counts when counts are high. Green asterisks indicate dates where at least one day within a 3-day window was significant by Gi*; purple crosses indicate significantly clustered cases within PCs.

Localized hotspots identified by kriging comprise a narrow range of PCs in each FSA. The area of interpolated contours that overlapped FSAs with significant disease burdens (by Gi* analysis) was compared to the total area of all corresponding FSAs. COVID-19 cases were highly localized, as the highest value contour (where kriging interpolates the greatest number of cases) constituted only a small portion of the total FSA area (Extended data,²⁹ Section 1 - Figure S4). While hotspots found by kriging often coincide temporally with Gi* results (Extended data,²⁹ Section 1 - Table S4), the interpolated hotspots contain fewer cases relative to the entire FSA, since the maxima of kriging contours cover only a portion of the overlapping FSA. The overall area indicated as a region of interest by kriging overlaps a greater fraction of the FSA area, but still covers a fraction of its area (<10% in the majority of cases). Compared to the other COVID-19 waves, however, kriging analysis of cases during this wave tended to span a larger area of each FSA (Extended data,²⁹ Section 1 - Figure S4A). This appears to be consistent with higher transmissibility (which, we speculate, may be related to lower levels of immunity during this period).

COVID-19 case counts were stratified into three categories: biological sex (males vs. females), age (exceeding vs. below 65 years old), and overall health (presence or absence of at-risk, pre-existing comorbidities). Stratified PC-level case counts were evaluated by Cluster and Outlier analysis across all waves ( Table 1). In the first wave, the frequency of high-case outliers (a high-case PC with few cases amongst its immediate neighbors) is biased towards older, at-risk populations ( Table 1A). Across all municipalities, PCs were flagged nearly 5-fold more often in the older population (305 and 62 PCs were identified as a high-case outlier on different dates in this wave, respectively). A similar pattern was observed in at-risk vs healthy populations (321 and 31 high-case outlier PCs, respectively), which is not surprising due to the considerable overlap between the elevated age and groups with at-risk comorbid diagnoses. High-case clusters were identified less often, but also showed this pattern, which is likely related to the high frequency of outbreaks in long-term care homes during this period (Table 1B). Windsor/Essex was an exception to this trend, since the frequencies of PCs designated as high-case outliers were comparable among both age and comorbidity groups (20 vs 17 for older vs younger populations, 19 and 17 for at-risk vs. healthy populations; Table 1A). The frequency of high-case clusters/outliers identified by sex-stratification varied by the region analyzed (e.g., more high-case outliers were identified in the female population in Toronto, in contrast within the male population in Windsor/Essex). This stratified analysis itself does not provide insight into the mode of COVID-19 transmission, that is, the development of a hot spot consisting of a particular subgroup does not indicate that viral transmission occurred exclusively between members of the particular subgroup.

Table 1.

Number of postal codes with high case clusters and outliers by location and wave.

A. Distribution of high-case outlier postal codes in different municipalities
City1 [No. of FSA/ PC/Area2/Pop.3]	Wave4	Estimated total cases5	Number of high-case outlier PCs in stratified populations6
			All cases	Males	Females	≥65 years	<65 years	At-risk	Healthy
(HA) Hamilton [30/ 14,774/1,138/536,915]	1	787	8	3	6	7	1	8	1
	1-2	410	0	0	0	0	0	0	0
	2	8788	100	20	43	31	45	41	24
	2-3	1464	29	6	4	0	25	2	8
	3	9755	150	36	43	0	131	9	84
(KW) Kitchener/Waterloo [18/8,445/200.9/338,205]	1	511	25	14	11	17	4	24	0
	1-2	277	0	0	0	0	0	0	0
	2	4287	31	3	10	12	13	17	6
	2-3	560	11	5	0	0	11	0	8
	3	2733	30	3	0	0	26	0	15
(LO) London [20/10,063/420.6/383,825]	1	363	11	3	3	6	5	6	3
	1-2	213	6	0	2	0	6	0	0
	2	3609	31	3	5	4	21	5	7
	2-3	347	2	0	0	0	0	2	0
	3	4363	28	6	1	0	26	0	10
(OT) Ottawa [47/ 21,963/2,790/934,031]	1	1125	50	10	26	42	2	45	0
	1-2	1327	28	2	10	8	16	11	12
	2	6790	120	25	21	28	84	29	54
	2-3	1178	12	0	2	2	7	2	7
	3	8716	50	8	13	0	48	0	30
(TO) Toronto [101/47,450/630.2/2,732,000]	1	7144	261	81	162	213	33	219	10
	1-2	3148	25	3	10	6	13	10	4
	2	54,460	596	78	178	191	315	216	127
	2-3	7309	89	17	13	1	76	16	40
	3	52,037	428	45	72	3	357	25	181
(WE) Windsor/Essex [24/ 10,758/1,785.1/415,293]	1	768	42	25	13	20	17	19	17
	1-2	1027	35	27	9	3	34	3	30
	2	7804	212	83	91	86	121	116	68
	2-3	530	6	3	0	0	3	3	3
	3	1800	43	11	8	0	41	2	31

1 Map of the southern half of Ontario with two-letter abbreviations of each region of interest indicating their locations within the area. Additional maps below display the border and PC boundaries within each region.

2 Area of region (km²).

3 Population based on 2016 census.

4 Interwaves 1-2 and 2-3 are displayed as ‘1-2’ and ‘2-3’, respectively.

5 Estimated COVID-19 cases in FSAs associated with region (days where cases within FSA were masked were assumed to be 1 case).

6 PCs identified as a high-case outlier on a particular 3-day interval with ≥ 6 cases when analyzing stratified case data.

Interwave 1-2

The first COVID-19 wave ended approximately two months after the imposition of public health prohibitions on mobility, limits on assembly, and social distancing (i.e., a lockdown) in Ontario. June 16^th, 2020 defined the start of interwave 1-2, a period of much lower COVID-19 incidence than the preceding 2 ½ months. Fewer FSAs exhibited high disease burden by Gi* analysis (39.7% of those called in wave 1, despite the duration of interwave 1-2 being longer; Extended data,²⁹ Section 1 - Table S2). Despite the overall decrease in cases province-wide, the frequency of FSAs in the Windsor/Essex region that were designated hotspots by both Gi* and kriging analysis increased relative to the preceding wave (Extended data,²⁹ Section 1 - Table S4). This was due to multiple outbreaks in the Essex County FSAs in Leamington [N8H] and Kingsville [N9Y] (Extended data,²⁹ Section 1 - Figure S5). These outbreaks were publicly documented in guest farm workers, which is reflected by the Cluster and Outlier analysis of stratified case data (a strong bias of high-case clusters and outliers in healthy males <65 years in Windsor/Essex; Extended data,²⁹ Section 1 - Table S4). These same locations were correlated with interpolated kriging contour maxima and often consisted of isolated single PCs (Table 1). However, there were instances where the interpolated count values exceeded the actual case count (i.e., N8H, June 22-26^th, 2020; Extended data,²⁹ Section 1 - Figure S5B) due to concurrent high case counts from neighboring FSAs (i.e., N0P 2J0; Extended data,²⁹ Section 1 - Figure S5C) bleeding into these contours from a single PC.

Wave 2

FSAs with high disease burden were more frequent within the second wave (Extended data,²⁹ Section 1 - Table S2). Even when accounting for its longer duration, there was a >10-fold increase in FSAs with high disease burden based on the Gi* test. The distribution of cases was more widespread across each FSA, as COVID-19 testing became more reliable and pervasive. Localized hotspots were detected less frequently within FSAs with high disease burden, i.e., significant FSAs by Gi* analysis overlapped less often with kriging contours of interest (>5 interpolated cases) in all municipalities except Windsor/Essex (Extended data,²⁹ Section 1 - Table S4). Additionally, interpolated case counts by kriging tended to be less than actual FSA case counts, due to dispersion of cases across more PCs. Most municipalities contained more FSAs with localized hotspots that were not significant by Gi* analysis. Increased case counts across the province during this wave increased the significance threshold for Gi* analysis. This decreased the likelihood that FSAs in municipalities with moderate population densities were designated significant hotspots. For example, FSAs within Kitchener/Waterloo and London were not significant COVID-19 hotspots during wave 2 (Extended data,²⁹ Section 1 - Table S4). Hotspots in these regions were coincident with maximum contour levels interpolated by kriging. PCs had a greater proportion of high-case clusters and/or outliers in younger and healthy populations, in contrast with wave 1 (Table 1). The reduction of COVID-19 outbreaks (i.e., hotspots) in long-term care homes during this period was supported by reduced statistical significance of counts at these locations relative to increased background case counts across the province.

Localized hotspots detected in PCs by kriging were more frequent in Hamilton, London, Ottawa, Toronto and Windsor/Essex during this wave than in the preceding one (Extended data,²⁹ Section 1 - Table S4). For example, an interpolated hotspot in the N5Y FSA (London) with >5 interpolated cases over consecutive 3-day intervals (December 19-21 and December 20-22, 2020; Figure 6) represents a ‘streak’, since kriging analysis interpolated cases over this FSA throughout the previous week (however, it was below the kriging threshold for significance). This hotspot encompasses a 2 building apartment complex that was publicized on December 29 by the local health authorities as a COVID-19 hotspot, two weeks after it was first detected by kriging.³³ While the interpolated case count is lower than the total count in this FSA on this date (N=33, other cases are present at other locations), this hotspot encompasses only 4% of the FSA by area. Other localized hotspots during this wave are depicted (e.g. in Leamington ON; Extended data,²⁹ Section 1 - Figure S6) and catalogued in Extended data,²⁹ Section 2.

Figure 6. Example of a localized hotspot in London Ontario.

(A) Geostatistical analysis of the FSA N5Y (London, ON) identifies a localized hotspot from Dec. 19^th to the 22^nd. Estimates from kriging were lower than actual case counts. (B) The kriging contour of a localized hotspot (yellow) interpolated with ≥5 cases encompasses 4% of the FSA area. Dots indicate PCs at least 1 case, however only those with ≥6 cases are labelled (with PC identifier/case count).

Significant clustering of cases in PCs in close proximity based on Spatial Autocorrelation analysis may be distinct from other tests that identify localized hotspots. Spatially clustered localized COVID-19 hotspots may not be confined within the same region of an FSA. For example, FSA M1P (Toronto) from Nov. 1-3 to Nov. 3-5, 2020 exhibited significant case clustering based on Gi* analysis. However, a localized hotspot identified within this FSA during this time did not overlap this cluster (Extended data,²⁹ Section 1 - Figure S7A). A neighboring FSA M1C showed significant spatially correlated cases within the same timeframe (Nov. 6-8, 2020), however localized hotspots were not detected in this FSA (Extended data,²⁹ Section 1 - Figure S7B). No PC within this FSA had ≥6 cases during this period. This illustrates how combining Spatial Autocorrelation, Gi* and kriging analyses can improve the interpretation of asymmetric COVID-19 case distributions.

Interwave 2-3

The interval between the end of wave 2 and inception of wave 3 was short (~31 days), however by comparison with interwave 1-2, the case counts were considerably higher (averaging 1134 and 154 per day, respectively). Thus, a far greater number of high disease burden FSAs were identified during this period (1.9 and 4.7-fold more often compared to wave 1 and interwave 1-2, respectively; Extended data,²⁹ Section 1 - Table S2). However, most of these occurred in Toronto and Windsor/Essex regions only (Extended data,²⁹ Section 1 - Table S4).

Wave 3

The frequencies of high disease burden FSAs in waves 2 and 3 were comparable upon adjusting for duration of each wave (Extended data,²⁹ Section 1 - Table S2). The regions detected by geostatistical analyses were different. High disease burden FSAs implicated in Hamilton and Kitchener/Waterloo by Gi* analysis were significantly more frequent compared to the preceding wave, while FSAs within Toronto were less often significant (Extended data,²⁹ Section 1 - Table S4). By contrast, localized hotspots within FSAs were generally less frequent than in wave 2, most notably in the Windsor/Essex region which contained 389 localized hotspots in wave 2, but 80 in wave 3 (Extended data,²⁹ Section 1 - Table S4).

High-case clusters and outliers were seldomly identified in older (≥65) populations in wave 3 and interwave 2-3 (Table 1). Similarly, there is a considerable decrease in the number of case clusters/outliers in at-risk populations relative to previous waves. For example, analysis of at-risk populations in Windsor/Essex identified 116 PCs as a high case outlier on a particular date during wave 2, but only 2 in wave 3 (Table 1A). These observations may be related to effectiveness of vaccination, based on the initial prioritization of older and at-risk cohorts for immunization.

The distribution of and relationships between clustered PCs with persistent COVID-19 cases

COVID-19 distributions were evaluated by Cluster and Outlier analysis to identify neighboring PCs with high concentrations of cases. PCs deemed part of a high-case cluster that were clustered over contiguous date ranges were then identified. These PCs were recognized as clustered over the duration of a time “streak”, which by definition, allows for a single day in which the PC was either not clustered or did not meet the minimum case count threshold of ≥6 cases within the 3-day sliding window. PCs that were part of a high-case cluster, but did not meet the case threshold, were not reported. This analysis identified 20 unique PCs with one or more high-case cluster streaks in Hamilton, 5 in Kitchener/Waterloo, 10 in London, 24 in Ottawa, 216 in Toronto, and 19 in Windsor/Essex (Extended data,²⁹ Section 1 - Table S5). Streaks were uncommon in Kitchener/Waterloo and London, as PCs in these cities rarely fulfilled the minimum case threshold. Most streaks were shorter than 4 sliding window periods, indicating that the most frequent streaks were brief due the large number of single day discontinuities from case reporting patterns and ascertainment bias. For these reasons, ≤3-day bins were not considered evidence of a streak. However, 37 streaks longer than 3 days were interpreted to indicate persistent infections at the same locations. Streaks with more than 3 consecutive high-case clusters were uncommon in wave 1, with M6M 2J5 (West Park Long Term Care Center) as the only example of a streak exceeding 3 days. The lack of other locations with persistent infections may be related to inconsistent case reporting which was hampered by limited testing during this period.^4,34

We sought evidence of COVID-19 transmission by identifying pairs of PCs with high case counts that were adjacent to one another by both distance and time. In Figure 7A, we analyzed PCs and their 10 closest neighbors, where streaks were either concurrent, consecutive, or separated by short gaps within the same wave. Concurrent streaks of neighboring PC pairs spanned at least 3 and as many as 29 days in Hamilton, Ottawa, Toronto, and Windsor/Essex (Extended data,²⁹ Section 1 – Table S6). There were a greater number of discontinuous streaks across all municipalities, with gaps separating streaks in neighboring PCs varied from short to long intervals (Extended data,²⁹ Section 1 - Table S7). During the second wave, many pairs of PCs were comprised of neighboring apartment buildings, for example, within Toronto FSA M4H. Other long concurrent paired streaks (exceeding 1 week in duration) occurred in Toronto PCs M1R 1S9 and M1R 1T1, M1L 1K9 and M1L 1L1 (both PC pairs in North York, Toronto), and L0R 1C0 with its neighbors in Hamilton (Extended data,²⁹ Section 1 – Table S6). Pairs of concurrent streaks in close proximity did not occur during wave 1, since streaks were uncommon during this time. However, many neighboring PCs were both found to be clustered both spatially and temporally during interwave 1-2 and waves 2 and 3 (Extended data,²⁹ Section 1 - Table S7). While not concurrent, neighboring PC pairs in rural Essex county were clustered within 10 days of one another during interwave 1-2 (N0P 2L0, N0P 2P0 and N0R 1B0). This is confirmed by publicly reported outbreaks in agriculture facilities of this region during this time. During wave 2, 56 unique pairs of streaks were identified in neighboring PCs in Toronto (18 clustered concurrently, another 14 within a week of one other), 4 pairs in Hamilton (2 concurrent), 1 in London (concurrent), 4 in Ottawa (1 concurrent) and 14 in Windsor/Essex (10 concurrent; all involve PCs in Essex county). There were fewer incidences of streaks during the third wave, which may be related to the effectiveness of vaccination efforts which limit transmission (Figure 7B). During wave 3, we identified an additional 30 pairs of PCs with streaks in Toronto (7 are clustered concurrently, and an additional 4 pairs are offset by 1 week), and 6 in Hamilton (3 concurrent, and another pair offset by 2 days). Interestingly, 9 of the same PC pairs in Toronto clustered during the third wave were also flagged during the second wave, which could be an indication of an undetected reservoir of COVID-19 that persisted in this population during interwave 2-3. These pairs of PCs occurred in the Toronto neighbourhoods of East York (M4H 1J4 and M4H 1J5), Scarborough (M1H 2Y7 and M1H 2E9), North York (M1R 1T1 and M1R 1S9; M6M 5B3 and M6M 5B7), and Etobicoke (M9C 1G7 and M9C 1G6; M9V 3S6 and M9V 4A4; M9V 3S6 and M9V 4M1; M9V 4A4 and M9V 4M1; M9W 6L4 and M9W 6A7).

Figure 7. Contiguous or consecutive pairs of Toronto postal codes with COVID-19 case clusters lasting multiple days.

Pairs of neighboring PCs with high-case count clusters were identified by Cluster and Outlier analysis. Panels show: (A) parameters measured in PC pairs that were elements of significant high-case clusters over multiple consecutive days (allowing for a single day gap), termed “streaks”; (B) characteristics of neighboring pairs of PCs with multi-day streaks in close proximity in Toronto (of ten closest neighbors based on distance between centroids) during waves 2 (top) and 3 (bottom). Wave 1 is not indicated because streaks were infrequent. Each point represents a pair of PCs, with color codes (see legend) indicating ranges of total number of cases in both, the X axis indicates the duration between PC pairs (with negative values representing number of days of overlap between end of first PC streak and beginning of the second), and the Y axis indicates spatial distance between PCs constituting the pair.

We examined evidence for likely community spread of COVID-19 between infected individuals in different pairs of neighboring PCs. This involved determining the distances separating them, the duration between streaks in the respective PCs, the total case counts among PC pairs during these streaks, and then assessing the significance of these characteristics. During waves 2 and 3 in Toronto, paired concurrent or consecutive streaks in such PCs were common (Figure 7B). Frequencies of these PC pairs were negatively correlated with the duration of the gap between the end of a streak in one member of the pair and the beginning of the other (r = -0.57 and -0.50 for waves 2 and 3, respectively). To assess the significance of these results, a multinomial logistic regression model was derived in R (“Streak_Analysis_using_Multinomial_Logistic_Regression_Models.Rmd” which uses the multinom function of the nnet package; available in Zenodo archive²⁹), evaluating (1) the total number of cases in each pair of PCs over the course of both streaks (classified as either <25, between 25 and 49, or ≥50 cases; the response variable of the regression), (2) the duration or gap (in days) separating streaks in a pair of PCs, and (3) the distances between the centroids of both PCs. The duration between streaks was significantly shorter for PC pairs with ≥50 cases relative to those with <25 cases (p = 0.0005 and 0.0141 for waves 2 and 3, respectively, using the Wald two-tailed z-test on the coefficients of the model). The same relationship was also significant when comparing streaks in pairs of PCs with ≥50 cases and those with 25-49 cases (p = 0.0015 and 0.0277 for waves 2 and 3, respectively). Total case counts and time between streaks in a pair of PCs were also inversely correlated (r = −0.41 in both waves 2 and 3; using the CORREL function in Excel), which is consistent with disease transmission between adjacent PCs being more common when cases were more abundant. These parameters were also inversely correlated (r = -0.13 and -0.17, respectively) in Ottawa and Windsor/Essex, but under less stringent criteria that allowed ≥6 cases for only one of the PCs in a pair on a single day. The relationship between distance between PCs and total number of cases was not significant by the regression model in either of these waves. Increasing the number of PC neighbors considered over a larger range of distances did not alter this result (Extended data,²⁹ Section 1 - Figure S8). However, upon dividing groups of neighboring PCs thresholded according to average quartile values (distance [182 m] and average gap between streaks [17.5 days]) distance between neighbors was significant when combined with the degree to which they coincided temporally. Below average duration between streaks and distance separating PCs were 2.8-fold more likely to exhibit ≥50 cases than all other PC pairs (p = 0.002 by Mantel-Haenszel chi square test). Transmission of higher levels of COVID-19 between pairs of PCs is supported by abbreviated intervals between streaks and their proximity to one another.

Multiple paired PC streaks in Toronto were tightly grouped within a single neighborhood over a short period of time (Extended data,²⁹ Section 1 - Table S7). Directional network diagrams were created to illustrate the connectivity and temporal relationships between these and other paired high-case cluster streaks in close proximity (Extended data,²⁹ Section 1 - Figures S8 and S9 for waves 2 and 3, respectively). The preponderance of paired streaks were bimodal (9 paired streaks where only two neighbors interact [Extended data,²⁹ Section 1 - Figure S9D], of which 17 occurred in wave 3 [Extended data,²⁹ Section 1 - Figure S10A]), and in close succession (≤7 days between streaks). Multi-node networks involving ≥3 postal codes occurred in 3 distinct neighborhoods in wave 2, and another three neighborhoods within wave 3. These neighborhoods were often comprised of multiple adjacent apartment buildings, and rarely occurred in low-density residential areas. The PCs where streaks tended to coincide were often proximate to facilities where communities may have congregated (grocery stores, schools, parks, or religious sites). We investigated whether the participation of many PCs in a multi-node network was a result of higher case counts during a streak (which would suggest an increased probability of transmission between these locations). However, bimodal networks of two PCs with recurrent high case counts were not rare (e.g., >100 cases within a 9-day streak in M1R 1T1), so other factors are necessary to explain why some PCs form multi-node networks. Furthermore, the average numbers of cases within streaks in bimodal networks and multi-node networks were comparable (18.8 and 17.0 cases per streak, respectively).

In some instances, PCs began to accrue cases prior to their inclusion in significant clusters of PCs or reaching the minimum numbers of cases to be designated a COVID-19 hotspot (M4H 1K1, M4H 1K2, M4H 1K4, M4H 1L1, M4H 1L6, and M4H 1L7). These PCs were nevertheless elements of the M4H FSA streak. Similarly, M4H 1L3 and M4H 1L4 were elements of a cluster containing M4H 1J4 on Nov. 8-10, 2020, at least 2 weeks prior to the occurrence of the corresponding multi-day streaks. Consecutive streaks within a set of clustered PCs in Etobicoke (M9V 3S6, M9V 3Z8, M9V 4A4, M9V 4A9, M9V 4M1, M9V 4P1, and M9V 5G9; centroids were all within 0.63km of one another) occurred between mid-October and late December 2020 (Figure 8B). M9V 5G9 is separated from the other PCs by a wooded area, but all are close to a major thoroughfare, Kipling Avenue. A striking directional acyclic network of COVID-19 cases involved 7 adjacent PCs whose streaks occurred concurrently and were occurred nearly consecutively along Bathurst Street, Toronto (M2R 1Z2, M2R 1Z7, M2R 1Z8, M2R 1Z9, M2R 2A1, M2R 1Z6, M2R 2Y2; within 0.57km of one another; Figure 8C). Curiously, the order of occurrence of COVID-19 infections among these PCs consistently followed a southerly direction along Bathurst Street.

Figure 8. Multi-node directional networks of closely situated clustered postal codes.

Directional acyclic graphs (left) organize adjacent high-case clustered PC streaks within the same FSAs ordered according to their occurrence. Panels indicate PCs within (A) M4H and (B) M9V in Toronto during wave 2. Each connection (or ‘edge’) linking two PCs represents a pair of clustered streaks occurring within 30 days of each other (all significant PC pairs are indicated in Extended data,²⁹ Section 1 - Table S7). Labels on each edge indicate the duration (in days) separating adjacent PC pairs and the number of cases which occurred in the combined streaks (“dates.cases”). Negative values, when present, indicate the number of overlapping days of concurrent streaks in pairs of PCs. Arrows indicate the temporal order of these paired streaks (from earlier to the later streak). Maps (right) show the physical locations of each PC within the corresponding networks.

The multi-node networks during wave 3 consisted of only three postal codes (M3C, M9V, M9W; Extended data,²⁹ Section 1 - Figure S10A). Increasing the number of neighbors (by relaxing distance constraints, from 10 to 30 PCs) for these PC networks results in multi-node networks that were previously identified as bimodal (M4H, M4C, and M4X; Extended data,²⁹ Section 1 - Figure S10B). This procedure also resulted in four additional bimodal pairs of PCs in wave 2 (Extended data,²⁹ Section 1 - Figure S9D); however, the total number of PCs comprising the multi-node networks in wave 2 was unchanged (Extended data,²⁹ Section 1 - Figures S8A-C).