2.1 Time evolution and early patterns of spread

We identified 385,558 cases between November 15^th 2021 and January 6^th 2022, of which 227,286 were likely Omicron. From 1^st May 2021 to 7^th September 2021 we identified 269,838 cases, of which 229,073 were likely Delta. The remaining cases in these periods (those with no S-gene result, or a different result) are excluded. The start date for each of these periods is the first date from which there are consistent rises in cases that are likely the new variant.

Omicron cases had a doubling time (the time taken for newly reported daily cases to double) of 2.9 days over the first 28 days, compared to 6.2 days for Delta (Fig A in S1 Text). Over half of all DZs had reported an Omicron case in the wave within 29 days, whereas for Delta this took 39 days (Fig B in S1 Text).

The reproduction number R_t consistently rose for Omicron, peaking at above 2 for nearly all local authorities 28 days in to the outbreak, and only consistently falling below 1 after 50 days (Fig C in S1 Text). Reproduction numbers for Delta are less consistent between LAs; while the number generally remains above 1 for most LAs in the period, there is no coherent peak at the start of the wave.

In the intermediate period during which Omicron became dominant and Delta declined, the age distributions by variant differed (Fig D in S1 Text). Taking the mid-points of the five-year age brackets, the mean ages of the Delta-type cases was 3.9 years lower than the Omicron-type cases (31.8 years compared to 35.7 years). A Student’s t test shows this difference to be statistically significant (t = −52.2, p < 0.001). This was the case from relatively early on when Omicron accounted for at least 5% of cases. However, the median ages are equal (both 32.5 years), as in the Omicron-type cases there is a trough in those aged 0–14, with fewer than 50% of cases in this age group Omicron, but then a peak in the 20–29 age group.

2.2 Case distribution and model fit

Fig 1, shows the distribution of COVID-19 cases for the Omicron and Delta waves broken down by age, sex, prior cases (serving as a proxy for prior immunity from infection), deprivation and health board. Omicron case rates were highest in younger adults, peaking at 90 cases/1,000 in ages 20–24. There was only a small difference in rates between men and women. Case rates were much lower amongst those that had tested positive for COVID-19 previously. Fig 2 shows case rates per DZ. Geographically, case rates fall with increasing rurality, most notably in Orkney, Shetland and the Western Isles (all island communities). The trend with respect to multiple deprivation decile is bimodal, with higher rates towards the highest and lowest deciles.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

Summary of 227,286 Omicron COVID-19 cases in Scotland between November 15^th 2021 and January 6^th 2022 (blue, filled), and 229,073 Delta cases from 1^st May 2021 to 7^th September 2021 (green, filled). The full population (N = 5, 465, 169) is broken down by age range, prior case status (whether a person had previously reported a COVID-19 case prior to that specific wave, and when), deprivation (of place of residence, per the SIMD decile, with 1 the most deprived), rurality (of place of residence, per the census Urban/Rural Classification) and location (at the level of Scottish health board). Cases are given per 1,000 people in that group (with subpopulation N recorded on the axis labels). The corresponding case rates as fit by our models are superimposed. Note that the subpopulations in the prior case status plot change across waves, due to being at different points in time.

https://doi.org/10.1371/journal.pcbi.1011611.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. COVID-19 cases in Scotland over the Delta period (A) as compared to Omicron (B), with focus on the Greater Glasgow region (C, D).

Each point indicates the population-weighted centroid of a DZ, with the colour representing the number of cases reported. Base maps obtained from Natural Earth [39].

https://doi.org/10.1371/journal.pcbi.1011611.g002

The fit case rates from our random forest regression models are overlaid onto Fig 1. We achieve a good fit to these larger-scale trends. The model slightly under-fits the age ranges 15–24, where case rates were the highest overall. Variable importance outputs are presented in Fig G in S1 Text, with node purity and accuracy loss.

Fig 3A shows model performance at DZ level, comparing observed cases to fit cases. Beginning with Omicron cases, our full model explains 70% (fit: 71%, test: 62%) of local variation in the case distribution (R-squared for case numbers, aggregated at a DZ level), with a poorer fit for cohorts with very high case counts. A “reduced” random forest model informed by population and population density alone explained 59% (fit: 60%, test: 55%) of variation. A model informed by only population/deprivation rank explained 53% (fit: 53%, test: 51%), and one informed by only population/age explained 48% (fit: 48%, test: 51%). Fig 3A shows further deviation of the data-fit slopes away from the diagonal for these “reduced” models.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Performance of different models.

(A) comparing observed cases to fit cases at DZ level. Each point represents a DZ. Points deviating from the diagonal indicate DZs with less accurate fits. The full model is compared with performance of reduced models informed with only population, and one of either age, overall deprivation rank, or population density. Also shown is residual clustering as measured by the Moran’s I statistic, at different physical distances (B) and network-based distances (C). Higher values represent higher autocorrelation between model residuals, when comparing DZs sitting within a given locus. DZs are defined as nearest neighbours of one another if they share a boundary.

https://doi.org/10.1371/journal.pcbi.1011611.g003

Considering now earlier Delta cases from 1^st May to 7^th September 2021, the geographical distribution (Fig 2) is visually similar, with a concentration of high case rates in the denser “central belt”. Cases skewed slightly younger (Fig 1), with the highest rates within ages 15–19. The distributon with respect to deprivation decile remains bimodal, with higher rates in both the most and least deprived DZs. Model performance was similar, explaining 72% (fit: 73%, test: 61%) of DZ-level variation.

Fig 3B and 3C shows for both the Delta and Omicron models, autocorrelation of residuals (as measured by the Moran’s I statistic, Section 4.5) within 1km is 0.35, falling to 0.15 at 5km, and 0.05 at 50km. The reduced models exhibit much higher residual autocorrelation, with the density-only model performing best, but persisting over larger distances (see Fig F in S1 Text for a map view of residuals).

2.3 Accumulated local effects

Fig 4 shows the accumulated local effects (ALEs) of all explanatory variables in the model (see Section 4.4 for definition).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Accumulated local effects across all explanatory variables.

For each variable, the x-axis represents the range of values of that variable in the data, and the y-axis (note scale differences for population, age, sex and prior case status) is the ALE for that variable value. The overall magnitude of the ALE represents the relative size of the effect.

https://doi.org/10.1371/journal.pcbi.1011611.g004

Population, age, sex, and prior case status have ALEs that follow the empirical distributions observed in Fig 1; ALEs are strongly positive for ages between 15–40, and those that had never reported a case before.

Beyond these variables, Fig 4 shows that features such as low population density, high vaccination uptake, a low mean household size, and a low rate of negative LFD test reporting are protective. We note that for vaccination uptake, the protective value at zero is likely an artefact arising from cohorts with ages 0–9 that were not eligible.

The effects for many variables associated with social deprivation such as the ratio of working age people with no qualifications and the rate of income deprivation (see Section B.2 in S1 Text for full descriptions) are weaker. This is consistent with the small degree of deprivation-level variation seen in Fig 1.

The directionality of the ALEs remain broadly consistent across both waves. Some risk factors were more pronounced in the Delta model, including in mean hosehold size, population density and the proportion of individuals belonging to a black or minority ethnicity. Conversely, cohorts with very high student populations were associated more strongly with high case rates in the Omicron fit.

3.1 Risk factors

We presented the accumulated local effects (Fig 4), revealing broad indicators for higher or lower case rates, and how they changed between waves. It is difficult to fully disentangle whether a difference was caused by a change in control measures, or a change in virus strain. Nonetheless, our analyses provide some important insights.

To begin, high mean household size emerges as a risk factor, consistent with the high secondary attack rates for SARS-CoV-2 [40, 41], and increased risk of inter-household transmission relative to contacts outside of the home [42]. That this, and high population density are both stronger risk factors for Delta may reflect the stronger NPIs at this tme increasing the proportion of within-DZ or within-household transmissions.

High vaccine uptake (amongst those eligible) is also protective, more so with Delta, consistent with higher rates of immune breakthrough with the Omicron variant as compared to Delta [43–45]. We do not know the specific vaccination status of those in the test data, however, and linked data may show a stronger protective effect.

For Delta, a high proportion of individuals of black and minority ethnicity is a stronger risk factor. In the UK, this is also a risk factor for severe COVID-19 outcomes [46–48] but without detailed, linked data, it is difficult to firmly establish drivers for a heightened risk during the Delta wave. Differences may emerge from known variations in vaccination uptake [49] and occupation [50] (thus ability to work from home or effectively physical distance), and the relative impacts of those factors changing across the two waves.

Finally, living in a deprived community was suggested from early on [51] and has since also emerged as a risk factor for severe COVID-19 disease [52–57]. However, the corresponding ALEs for the variables associated with deprivation are small. Deprivation effects may be captured by proxy with other variables that correlate with deprivation such as age [58] and vaccine uptake [59, 60].

3.2 Testing frequency

The low case rate variation with deprivation (Fig 1) contrasts with observed inequalities over severe outcomes [22–25, 61], suggesting that those living in more deprived communities experience a higher inherent case-hospitalisation rate. We suspect that a lower proportion of case ascertainment, however, may also be a factor.

An important and unique variable in our model is the rate at which negative LFD tests were reported throughout the period. We found high rates of negative test reporting to be a risk factor. This suggests a variation in case ascertainment across different demographics, which may in turn lead to skews in the observed case distribution [22, 62, 63].

Further work (Fig H and Table A in S1 Text) shows that up to February 2023, the rate of LFD testing and positivity varied substantially across deprivation (quintile 1: 3.6 tests/person, 4.61% positive; quintile 5: 6.7 tests/person, 3.57% positive) as well as sex (M: 3.7 tests/person, 4.82% positive; F: 7.0 tests/person 3.30% positive). If demographic differences in testing behaviour correspond to differences in case ascertainment, the profile of all infections may then be biased from reported cases, and testing rates may be obscuring the true patterns of infection over sex and deprivation.

In addition, the magnitude of the risk factor (as seen in the ALE, Fig 4) plateaus beyond a certain rate (>∼1 test/person in each period). This hints at a deeper relationship between true incidence, the frequency of testing (and whom amongst the population is taking those tests), and the proportion of infections that are ascertained.

Our model is unique in including negative test reporting, and has revealed strong differences between different demographics that may bias the profile of cases. Beyond the work presented here, further analysess of reported cases need to be considered with these strong skews in testing behaviour in mind.

3.3 Conclusion

The COVID-19 data studied here are remarkable in terms of volume and resolution, and has allowed us to assess a national-level epidemic at extremely fine scale. However, regardless of resolution, cases only partially represent the full underlying pattern of infection. Variations in testing frequency and known trends in severe outcomes suggest that the distribution of infections may have been very different to that of reported cases. By incorporating trends on cases, testing behaviour, and severe outcomes more closely linked to infection (hospitalisation, ICU admission and mortality), it may be possible to build a much more comprehensive retrospective picture of how infections were distributed amongst the population.

Importantly, while our access to such finely-grained data was exceptional, it can be expected that such data are likely to become more common in the future, and may become available in real time. As such, our demonstration of the utility of such data points the way to an important approach to improving data analysis supporting control policy response to infectious disease emergencies in the future.

4.1 Preparation of case data

We use COVID-19 testing data from Public Health Scotland’s electronic Data Research and Innovation Service (eDRIS) system, dated from July 14^th 2022. The data include individual tests by type (polymerase chain reaction (PCR) or rapid lateral flow device (LFD)), test result (positive, negative, void, inconclusive), test date, S-gene test result if known (positive, dropout, inconclusive), age, sex, and residing data zone (DZ, a census area typically comprising 500–1,000 individuals). De-identified IDs link repeat tests by the same individual. We reduce the raw test data to cases by removing duplicate tests by the same individual within 60 days (taking the date of the first PCR positive as the case date, or the first LFD in the absence of any PCR). These metadata—in particular the DZ, specifying location to within an area as small as 0.1km² in densely populated areas—therefore identify cases at a fine spatio-temporal scale. Data on vaccine administrations are also provided by eDRIS.

This analysis considers the BA.1 sub-variant of the Omicron lineage only. The sub-variant BA.2/B.1.1.529.2 later replaced BA.1, becoming dominant in Scotland from around 25^th February 2022. This variant, like Delta, has an S-gene positive test signature. However by the end of the period studied the BA.2 variant was only being identified in fewer than 1% of fully sequenced cases in the UK [64], and here we assume all remaining S-gene positive cases to be Delta.

Prior to January 6^th 2022 in Scotland, positive LFD tests (typically taken at home) required PCR confirmation. Approximately 90% of cases in this period have a definitive S-gene result. A policy change then dropped this PCR requirement [65], after which cases with S-gene results fell to about 50% by February 2022 (per eDRIS data).

For Omicron cases, we gather from the data S-gene dropout cases between 15^th November 2021 and 6^th January 2022, and for the Delta outbreak, S-gene positive cases between 1^st May and 7^th September 2021 (choosing this end date to have a similar number of cases in each set). We exclude cases that have a different, or no S-gene result.

Using the linked historical tests, we label cases based on whether the individual had either: never tested positive before; had tested positive in the last six months prior to the start of that wave, or; last tested positive over six months prior to the start of that wave. We denote this the prior case status, as a proxy for infection-based immunity.

Finally to prepare the cases data to be fit, we group individuals that have the same age range, sex, residing datazone, and prior case status, terming these subsets of individuals cohorts. As an illustrative example, a cohort may be a population of 38 males aged between 50–54 residing in a given datazone “X”, that have never tested positive for COVID-19 before, among whom 9 Omicron COVID-19 cases were identified. This is the highest practical resolution we can acheive using the eDRIS case data, and our model (Section 4.3) fits case counts at this resolution.

4.2 Time series analysis

4.3 Model

Our statistical model is designed to explain variation in COVID-19 case numbers as prepared in Section 4.1, and identify risk factors amongst a broad range of variables, using random forest regression. We fit models to the distribution of Delta and Omicron cases respectively, allowing for comparison of risk factors across the two waves.

4.4 Accumulated local effects

To identify risk factors amongst the explanatory variables used to inform the model, we calculate the accumulated local effects (ALEs) of each variable. The ALEs describe how the model fit value changes, in response to changing one variable value in isolation, averaged over many different entries in the data [72]. In this context, ALEs indicate whether a variable value is associated with fewer or more cases in general over the data. If the ALE is greater than zero, the fit cases generally increases given that variable value.

4.5 Moran’s I autocorrelation statistic

To probe geographical variation in cases not explained by the model, we measure the Moran’s I autocorrelation [73, 74] on the residuals (the difference between the data and fit value), relating to their physical location. We compare local DZ-aggregated residuals over physical distances (from 1–100km), as well as network distance (number of nearest neighbours apart). For a set of N residuals y_i, the Moran’s I is a measure of autocorrelation: with the mean of all residuals, and w_i,j is an associated weight of the pair of observations (i, j), with w_i,i = 0. To measure the autocorrelation between residuals within a separation d (either a physical or network-based distance) of one another, we set w_i,j = 1 if dist(i, j) ≤ d, and 0 otherwise. Fully correlated residuals would have I = 1, whereas I = 0 would indicate no correlation.

This measure characterises how effective our models are at explaining geographical variation, and with different distances d shows over what length scales residual autocorrelation persists.