Assessing the Effect of School Closures on the Spread of COVID-19 in Zurich

Abstract

The effect of school closure on the spread of COVID-19 has been discussed intensively in the literature and the news. To capture the interdependencies between children and adults, we consider daily age-stratified incidence data and contact patterns between age groups which change over time to reflect social distancing policy indicators. We fit a multivariate time-series endemic–epidemic model to such data from the Canton of Zurich, Switzerland and use the model to predict the age-specific incidence in a counterfactual approach (with and without school closures). The results indicate a 17% median increase of incidence in the youngest age group (0–14 year olds), whereas the relative increase in the other age groups drops to values between 2% and 3%. We argue that our approach is more informative to policy makers than summarising the effect of school closures with time-dependent effective reproduction numbers, which are difficult to estimate due to the sparsity of incidence counts within the relevant age groups.

1 INTRODUCTION

Public health legislation authorises officials to implement disease control measures such as closing schools. The usefulness of school closures has been seen in certain infectious disease outbreaks but not in others. School closure for infectious disease control may have additional health and wider societal effects as school is not just education but also serves an important social function. Determining the usefulness of school closures is therefore of great value to policy makers. We know school closures have an effect on social mixing (Luca et al., 2018), and so are often used as a first-line approach for control of epidemics, particularly when no pharmaceutical prophylaxis is available. However, it would seem this effect might not be large for the early coronavirus disease (COVID-19) outbreaks (European Centre for Disease Prevention and Control, 2020). For this reason, we are interested in examining the alternative scenario where school closure to combat COVID-19 was not introduced among school-aged children.

To examine true social distancing interventions implemented, we fit a multivariate endemic–epidemic (EE) model to the observed case data from the Canton of Zurich, Switzerland incorporating social contact patterns. This is a proof-of-concept study of the feasibility of using time-varying contact matrices in EE modelling approaches as a means to support data-driven policy making. Our model includes an age structure which has been highlighted as important in models focusing on school closures (Jackson et al., 2014). We then predict from the fitted model given data until 17 March 2020 (the first day following declaration of a state of emergency) with assumed time-varying contact weights implemented (Figure 1) and (the alternative scenario) with time-varying contact weights ignoring changes to contacts due to school closures (see supporting information for a visualisation). We refer to our alternative scenario as a counterfactual (where the term is used in the philosophical meaning of the term—what has not happened or is not the case—rather than in a counterfactual analysis of causation). The changes in the alternative scenario affect the contacts of the youngest age group (0–14 year olds) in school; the age group that covers both compulsory and non-compulsory education. We then investigated the difference in the number of expected cases under the two scenarios to evaluate the usefulness of school closures.

2 METHODS

Finkenstädt and Grenfell (2000) showcased how to formulate a time-series susceptible–infected–recovered model through a case study of endemic measles infections and compared their methodology with established results from compartmental modelling, linking mathematical and statistical modelling. They highlighted a need to incorporate epidemic dynamics in statistical models; a need which was addressed by the EE modelling framework (introduced in Held et al., 2005).

2.1 Endemic–epidemic modelling

The EE framework is a time-series analysis-based method for infectious disease surveillance data. It can be derived from a mechanistic model of disease transmission (Bauer & Wakefield, 2018; Höhle, 2016; Wakefield et al., 2020), linking it with other modelling approaches. The multivariate formulation used in this work is the age-stratified EE model (Meyer & Held, 2017) where COVID-19 cases $Yat$ for age group a on day t are given by

where $λat$ is the mean and $ψa>0$ the overdispersion parameter of a negative binomial distribution, where the limiting case $ψa→0$ represents the standard Poisson assumption. Overdispersion is sometimes termed k in the infectious disease literature when using negative binomial distributions to examine superspreading (e.g. Endo et al., 2019; Lloyd-Smith et al., 2020).

The EE model is a two-component infectious disease model where the mean $λat$ is decomposed additively into an endemic (ν) and epidemic (ϕ) component. Age-specific population fractions $ea$ enter as known offsets in the endemic component, whereas the epidemic component depends on contact matrix weights $ca,a′,t$ representing transmission between age groups a and $a′$ on day t (for respiratory disease contacts are proxy for transmission events), and $ul$ is the discrete-time serial interval distribution (Bracher & Held, 2022). We use a shifted (normalised) Poisson distribution with weights $ul∝κl−1/(l−1)!exp(−κ),κ>0$ with a maximum lag of $lmax=7$⁠, as this has shown to be useful in other analyses of daily COVID-19 data (e.g. Grimée et al., 2022; Ssentongo et al., 2021). The transmission weights $ca,a′,t$ (entries of the contact matrix) are known but the serial interval lag distribution $ul$ (represented by the parameter κ) is not.

Various models were considered for the endemic and epidemic components of the model (see the supporting information for an overview of all models considered) and the best fitting model was determined based on the Bayesian information criterion (BIC) and model convergence. We always included information on public holidays (as contacts may differ on those days) and daily testing rates to account for possible temporal changes in underascertainment. In addition to this, we also included daily temperature, linear time trends and sine–cosine waves (a smooth non-linear trend not picked up by other parts of the model) as potential covariates in both components.

2.2 Time-varying contact matrices

Total contact matrices are a weighted linear combination of setting-specific contact matrices. We chose to use the Mistry et al. (2021) synthetic contact matrix compartments (Figure 2) in our work as they were the only ones we found provided with uncertainty estimates in the weights. Additionally, they were created with the European setting in mind, providing additional realism. We create Zurich-specific policy indicators following the methodology from Hale et al. (2021). As we are interested in reductions of contacts, we incorporate the policy indicators such that they take values between 0 and 1, where a higher value is a situation with less social distancing policy in place (Figure 3). This ensures that no metric used to incorporate policy changes increases contacts, thereby creating an artificially inflated baseline. While policy for households exists, artificial contact matrices are based on data on household size rather than surveyed contacts in the household setting and so we do not apply policy indicators to this setting, essentially applying a constant indicator of 1 to household contacts. We increase the lower bounds of workplace and other indicators to be 0.1 rather than 0 to reflect essential workers and necessary contacts, such as food shopping (see the supporting information for details on policy indicator construction). The results of applying these indicators to the compartments of the contact matrix and combining using the Mistry et al. (2021) weights are showcased in Figure 1. We additionally reduced contacts in school settings on school holidays in both versions of the time-varying contact matrices (for the two scenarios) as school holidays are known to cause a drop in contacts (Eames et al., 2012).

2.3 Counterfactual analysis

To examine the effect of school closures, we used a version of the time-varying contact matrix $ca,a′,t$ which does not have reductions applied to contacts in the school setting among the youngest age group, only regular school holidays. We predicted the course of the epidemic under two scenarios (true disease control measures as implemented; scenario A and adjusted to not have school closure; scenario B). The predictions are obtained using the methodology described in Held et al. (2017), Appendix A. To analyse the alternative scenario (scenario B), we adjust the transmission weights (time-varying contact matrix) in the EE model and calculate the predictive mean vector for the adjusted path forecast. That is, after fitting our EE model (described above), we predicted the epidemic from 17 March 2020 with the two options for time-varying matrices. To compare the two scenarios, we considered the absolute and relative increases in predicted cases between scenarios B and A. We expect both to be positive as scenario B (alternative scenario) is a deviation from the true non-pharmaceutical measures which were implemented and has a lower level of disease control in place, hence more contact/transmission opportunities.

2.4 Incorporating uncertainty

The importance of acknowledging uncertainty in COVID-19 modelling was highlighted by Davey Smith et al. (2020). We take into account both external (contact matrix weights) and internal (model coefficient estimates) parameter uncertainty in the counterfactual analysis. We used the weights reported by Mistry et al. (2021) to estimate a set of plausible values for the different contact matrix compartments to obtain total contacts across settings. We sampled the weights from normal distributions with the estimated means and standard errors of the reported weights used to create the time-varying contact matrix, allowing us to incorporate the corresponding external uncertainty. In particular, we simulated 1000 weights and created 1000 versions of the time-varying contact matrices.

Additionally, for each of 1000 versions of the EE model with those time-varying contact matrices, we sampled all entries in Table 1 using a 40-dimensional multivariate normal distribution based on the estimated parameter vector and the associated variance–covariance matrix. While Figure 4 and Table 1 show the results for using the Mistry et al. (2021) weights and the parameter estimates as given, our main results in Table 2 incorporate both external (contact matrices) and internal (model) uncertainties.

3 RESULTS

When EE models are too complex, some of the coefficients diverge; this is seen in large standard errors. Table 1 in the supporting information lists all the BIC values for models considered (including divergent models), but we selected the best model which had no diverging effects. We experienced identifiability problems when including temperature in the models. The final model has 40 parameters and log-linear predictors given by:

and

where $Tt$ is the daily testing rate at time t. We calculate the amplitude and phase shift of each sinusoidal wave based on the sine/cosine coefficients (following Held & Paul, 2012), see the supporting information for details.

The selected model has the lag distribution and fit shown in Figure 4. We see that the model captures the patterns observed in the data well, apart from the youngest age group (0–14 year olds), though this is likely an artefact of the low number of cases overall in that age group. The shape of the serial interval distribution (Figure 4, left) is right-skewed with a sharp peak early on. The peak found in this work is earlier than expected based on the literature and so we have conducted a sensitivity analysis with other estimates from EE models (Grimée et al., 2022; Ssentongo et al., 2021), see supporting information for details. The model coefficients are listed in Table 1 and the model fit is shown in Figure 4 (right). In Table 1, ν and ϕ denote the endemic and epidemic components, respectively. The coefficients (Table 1) show the expected pattern with a strong day of the week effect. We also see a seasonal pattern in particular in the epidemic component (amplitude and phase) which is counterbalanced by a positive time trend. There is considerable variation in the endemic and epidemic intercepts between age groups. Testing rate has a positive and significant effect in the endemic component whereas the effect on the epidemic component is less pronounced.

The result of the counterfactual analysis is given in Table 2. The distributions of our predicted counts are very skewed, which is why we provide median and percentiles as point and interval estimates rather than means and standard deviations. We see that the number of cases never increases more than 17% and—as would be expected—that most of the increase is found among the youngest age group (0–14 year olds). The largest increase in case counts is found in the age group 45–65, who were not considered a vulnerable group at the time. Based on the median, we see 128 additional cases (a 3% increase) with 9 additional cases in the oldest age group (which might be considered of most concern; those aged 80 and over). However, the 90% quantiles are considerably larger (329 additional cases overall and 25 cases in the oldest age group, respectively). Similar patterns of burden (which age groups have the highest case counts under the two scenarios) are seen in both scenarios. The large uncertainty seen is not so surprising as the data from Zurich is rather sparse. If we had instead applied the method to a larger population, such as the entire federation of Switzerland, precision would be increased.

We also compared the temporal dynamic in scenario B versus A. An increase in cases is observed in age group 0–14 already in April but only later in the other age groups. The next age group to experience an increase in cases is 25–44; the parents of age group 0–14. Thus, the dynamic follows the contact patterns and underlines the importance of including both age-stratified data and contact matrices in models. See the supporting information for an illustration and description of the patterns.

4 CONCLUSION

Understanding the effect of disease control interventions is useful for preparedness for future epidemics. We found it possible to utilise time-varying contact information in the EE modelling framework and fit an appropriate model to COVID-19 case surveillance data. We adjusted the contact matrix based on a timeline of COVID-19 events focused on Zurich. Our focus was school-based social distancing measures, as these might be considered a priority for policy makers when choosing exit strategies or phasing out of measures. In our counterfactual analysis we did not assume interventions were applied equally to the entire population but we acknowledge that changes may affect other parts of the population through contact patterns between age groups and so the effects of school closure are not to be considered in isolation.

In this work, we have assumed that the model parameters do not change when considering scenario B (the alternative scenario). However, some natural experiments were occurring where some schools are closed and other schools are kept open during the same time period in the same school district, see for example Berger et al. (2020) who examined whether increased testing can be used in conjunction with keeping schools open. This may allow us to re-evaluate that assumption in the future. Vlachos et al. (2021) have estimated the effect of school closures on the spread of SARS-CoV-2 virus (the agent causing COVID-19) among parents and teachers in Sweden, where lower-secondary schools (pupils aged 14–16) remained open during the first wave. Their results indicate that keeping lower-secondary schools open had minor consequences for the overall transmission, in line with our results.

The estimation of age- and time-dependent multiplication factors to address additional causes of underreporting (Noufaily, 2020), such as the presence of asymptomatic COVID-19 infections are being considered in our ongoing work. We are examining the option of using reporting rates as the basis of such multiplication factors. The rates are based on adjusting case fatality rates for delays between hospitalisation and deaths in the vein of Nishiura et al. (2009) and Russell et al. (2020). Such an approach allows us to address underascertainment (i.e. capture asymptomatic cases which are not likely to be reported as well as those symptomatic cases that are reported). We expect underreporting to likely be more pronounced among children because they may have more asymptomatic cases (causing them to have fewer cases in the data). Our model will need suitable amendments once such reporting rates are available.

Flaxman et al. (2020) studied the effect of non-pharmaceutical interventions on COVID-19 in Europe based on overall time-varying reproduction numbers (⁠$Rt$⁠). School closures are an intervention directly affecting only children and adolescents. Mortality in these age groups is notably very low, so our study aimed to quantify the indirect effects of school closures in other age groups based on age-stratified incidence data and appropriate contact information. We argue that risk communication strategies regarding disease control initiatives for COVID-19 should be based on age-specific (rather than overall) effective reproduction numbers, but those are difficult to estimate if data are sparse. This seems of particular importance if the interventions being considered are specific to certain groups of the population rather than overarching.

When a single summary indicator is reported to the public, uncertainty of its estimation should be included. For this reason, we believe our work to be particularly useful as we capture a lot of the statistical uncertainty inherent to our work. While the effective reproduction number is an attractive summary measure due to the threshold property (>1 indicates continued epidemic growth), some of the nuances behind its calculation will be lost to summarisation. These include the fact that $Rt$ can be approximated in different ways depending on the data at hand and model considered. This echoes Becker (2015) and recent comments made by the European Union’s Joint Research Centre (Annunziato & Asikainen, 2020; Gostic et al., 2020), as well as the aforementioned comment on uncertainty by Davey Smith et al. (2020). Additional context should be provided as setting-specific disease control challenges in which the estimate was obtained need to be communicated as well as $Rt$⁠, for example healthcare surge capacity, cases in nursing homes compared to cases in schools, and the impact of superspreading events although community transmission is low. We suggest other indicators, such as predicted case counts (as found here) are used alongside $Rt$ to provide broader context as it is not fully informative in isolation. For additional discussion, see the supporting information.

In essence, the EE framework allowed us to examine questions of interest to policy makers with the data available at the time but remains flexible enough that in future additional information could be included, for example the inclusion of vaccines as another model covariate, underreporting through adjustments to the case counts, and mask usage could also be included in our model as a different type of reduction to contacts as a proxy for its effect on transmission events. In this work, we found the results are somewhat sensitive to the choice of transmission weights. For this reason, we highlight that while the Mistry et al. (2021) weights have school as the most important setting, this seems to be in in contrast to a recent review by Mousa et al. (2021) which found that most contacts occur outside of school settings. As this work is a first-order approach, further analysis will also consider sensitivity analysis of contact matrices. In our future endeavours, we intend to construct an EE model to examine similar questions at national (Switzerland) level, which involves further intricacy as many of the decisions (as well as dates of holidays) are determined at cantonal level rather than federal level.

CODE AND DATA

Code and publicly available data used in this manuscript can be accessed via https://gitlab.switch.ch/suspend/COVID-19-school-ZH. Project updates and associated analyses are available via https://suspend.pages.switch.ch/project. For more information and extended descriptions of our work, see the supporting information.

Supporting Information

Additional supporting information may be found in the online version of the article at the publisher’s website.

ACKNOWLEDGEMENTS

This work would not be possible without the financial support of the Swiss National Science Foundation (https://data.snf.ch/covid-19/snsf/196247) and the confidential data provided by Gesundheitsdirektion Zürich. They did not play any role in the formulation of this manuscript and the views expressed in this manuscript are those of the authors. The additional members of the SUSPend modelling consortium (in anti-alphabetical order by surname) are: Jon Wakefield, Sebastian Meyer, Niel Hens, Mathilde Grimée, Deborah Chiavi and Johannes Bracher (https://suspend.pages.switch.ch/project).

REFERENCES