Using Population Based Kalman Estimator to Model COVID-19 Epidemic in France: Estimating the Effects of Non-Pharmaceutical Interventions on the Dynamics of Epidemic

Hospitalization Data

Hospitalization data are extracted from the SI-VIC database (Système d’Information pour le suivi des VICtimes), a governmental system created in 2016 in order to identify and follow victims in exceptional circumstances (e.g. terror attacks). Since March 18^th 2020, the SI-VIC database provides, to Santé Publique France, the daily number of hospitalized COVID-19 patients, at multiple geographical scales. In this work, we focus on the 12 French non-insular regions.

An entry in the SI-VIC database corresponds to a patient hospitalized in connection with COVID-19. Specifically, it requires the presence of at least one of two criteria: i) a biologically confirmed diagnosis of COVID-19 (e.g. RT-PCR positive test result), or ii) a chest CT scan suggestive of COVID-19. In this analysis, we rely on the daily incident number of hospitalizations , and the total number of individuals hospitalized daily (Y ^H). The two data series are displayed in Figure 1 over a period of 391 days (from March 2, 2020 to March 28, 2021). Of note, in order to be able to compare the magnitude of epidemic in each region, we standardized the data by the size of the population in each region: direct interpretation would be the number of incident or prevalent hospitalization for 100,000 inhabitants.

• Download figure
• Open in new tab

Fig 1.

Top: standardized total prevalent number of individuals hospitalized daily for 100,000 inhabitants . Bottom: standardized daily incident number of hospitalizations for 100,000 inhabitants .

Non-Pharmaceutical interventions

Timing and modalities of the various NPI implemented in France over the course of the epidemic have been gathered from the French government’s action summary website. In France, public health interventions have been highly multi-pronged. In our analysis, we considered the following summarized NPI occurring during the first year of the epidemic in France: i) first lockdown (with two phases of relaxation/reopening as described bellow), ii) second lockdown (with one phase of relaxation/reopening as described bellow), iii) 8PM curfew, iv) 6PM curfew, v) school closure, vi) bar and restaurant closure, vii) barrier gestures (including all mandatory sanitary protocols: physical distancing, hand washing, part-time remote work and, mask wearing in public spaces). Of note, NPI such as travel bans, enhanced testing, contact-trace-isolate, were ignored to ensure identifiability, as they were either i) with a magnitude difficult to quantify, or ii) enforced in complete overlap with other NPI. In addition, partial interventions at a sub-regional level were not considered as an implementation of a measure. This resulted in a rather similar profile of interventions across regions as most of them were applied simultaneously in the 12 regions of interest. Figure 2 summarizes all 10 considered NPI over time.

• Download figure
• Open in new tab

Fig 2.

Major NPI at the regional level

The first and the second lockdowns have been disconnected due to their different modalities resulting in different behaviors, and thus possibly different impact on transmission. For example, during the first lockdown (from March 17, 2020 to May 11, 2020), the entire population had to work from home – the only exceptions were for workers in vital sectors, such as medical, security, or food sectors – and outings could not exceed one hour in a one kilo-meter perimeter. Whereas during the second lockdown (from October 29, 2020 to December 15, 2020), working on-site was authorized when working from home was not feasible, and outings were limited to only 3 hours in a 20 kilometers perimeter. In addition, the end of the first lockdown was gradual and separated by the government in three official phases (phase 1: May 11, 2020 to June 2, 2020, phase 2: June 2, 2020 to June 22, 2020 and phase 3: after June 22, 2020) with many evolving measures as the authorized distance of travels, the reopening of cultural places (i.e., museums), the reopening of non-essential stores, etc. A government campaign to raise awareness about barrier gestures started at the end of the first lockdown. Masks and hand washing were mandatory in many places such as public transports, schools and companies. We assume that barrier gestures start on May 11^th 2020, assuming that most of mandatory measures with a potential strong impact started to be implemented at this date. Concerning the second lockdown, a reopening of the non-essential stores happened 2 weeks prior to the end of the lockdown prior to Christmas. To account for this, we separate the second lockdown in two phases: a full lockdown until November 28^th 2020 and a reduced lockdown afterwards.

School closures were documented as of the holiday schedule, which can vary between regions. In France, schools were also closed during the first lockdown. Furthermore, the reopening of schools from May 11, 2020 (end of first lockdown) to July 4, 2020 (end of term) was very progressive with phased school and level-reopening and attendance rising slowly back to normal. We averaged the estimated school attendance to be 30% of school capacities during that transition. Complete closure of all of bars and restaurants happened twice over the study period: first a few days before the implementation of the first lockdown in March 2020 (as for schools closure), and along the beginning of the second lockdown in October 2020. The measure was lifted progressively in all regions after the end of the first lockdown in June 2020. At the end of the second lockdown, the measure was not lifted considering the epidemic situation was not good enough. It was still in place at the end of the study period.

Curfew measures started to be implemented October 17, 2020 in the 12 regions of interest in several major cities and in Île-de-France, with a curfew from 9PM to 6AM. It was extended on October 22^nd to 54 departments (sub-regional administrative unit). It was suspended during the second lockdown from October 30^th to December 15^th 2020, when a new one was implemented at 8PM at a national level. From January 2-12, 2021, the starting hour was changed to 6PM progressively at the departmental level until January 16^th when it switched to 6PM for all 12 regions of interest. Finally on March 20, 2021, the starting hour was changed again nationally to 7PM. To model the measures we considered only two variables, grouping curfews starting at 9PM or 8PM and those starting at 6PM or 7PM, and we set the value to 1 only when a whole region was under the curfew. Assuming that curfews and lockdowns induce different behaviors, curfews are not considered to be included in lockdowns, even if it is forbidden to go out in the evening during lockdowns. Similarly, 6PM curfew and 8PM curfew were considered as different interventions instead of nested ones as they are likely to induce different behaviors.

VoC of SARS-CoV-2

Some variants of the SARS-CoV-2 virus have been classified by national and international health authorities as VoCs as they impact transmissibility or virulence, or decrease the effectiveness of measures World Health Organization. Since January 2021, French health authorities have conducted surveys to estimate the prevalence of three VoCs: 20I/501Y.V1 (Alpha), 20H/501Y.V2 (Beta) and 20J/501Y.V3 (Gamma). The survey of VoC Delta started after this study. We therefore included the cumulative proportion of cases infected by one of these VoCs as a potential covariate explaining the transmission. We used data from two cross-sectional “flash” surveys performed on January 7-8, 2021 and January 27, 2021, see Santé Publique France; Santé Publique France and the weekly estimation of VoC circulation provided by the SI-DEP database at the regional level from February 12, 2021 to March 28,2021.

Between January 8, 2021 and February 12, 2021, the estimated proportion of the sum of the three VoC have increased from a national mean of 3.3% to a national mean of 46%. To complete the missing data, we assume that the proportion before January 8, 2021 equals to 0% and that the evolution is linear between January 8 and January 27, 2021, and then between January 27 and February 12, 2021. Of note, a logistic and an exponential growth were also tested with no significant changes in conclusions (results not shown). With no data reported for Bourgogne-Franche-Comté region on January 27, 2021, only one slope was estimated on the time window. This variable is denoted V oC_i(t) in what follows. Representations of the VoC proportion in each region over time are given in Section 1 of the Web Supplementary Material.

Vaccination process

Vaccination started in France on December 27, 2020. Three COVID vaccines were at that time authorized: BNT162b2 mRNA (Pfizer), ChAdOx1 nCoV-19 (AztraZeneca) and mRNA-1273 (Moderna). To take into account the vaccination process, we used the VAC-SI database which provides the cumulative percentage of the population who have been vaccinated with at least one dose of vaccine over time. Vaccination was first targeted to the elderly 75+ years of age. The proportion of vaccinated individuals, denoted here within as V, ramped up to about 12% by the end of study period. Representations of the vaccination coverage in the population over time are given in Section 1 of the Web Supplementary Material.

The mechanistic model

We model the evolution of the COVID-19 epidemic using an extended SEIR type model Wang et al. (2020), called a SEIRAH model, adapted from Wang et al. (2020); Prague et al. (2020) where the population of size N is divided into 5 compartments: susceptible S, latent exposed E, symptomatic infectious I, asymptomatic/pauci-symptomatic infectious A, hospitalized H, removed R (i.e. both recovered and deceased), see Figure 4. The number of vaccinated people denoted by V is assumed to be known, see Section 2.2. The dynamics of such model is given by

where α, r_E, D_E, r_I, D_I, D_Q, D_H are time-independent parameters described in Table 1 while b is a function of time modeling the disease transmission rate.

• Download figure
• Open in new tab

Fig 4.

SEIRAH model representation – adapted from Wang et al. (2020); Prague et al. (2020)

The observation model

The two quantities Y ^H and relate to the solutions of System (1) respectively as for all regions i = 1, …, 12, for all observation time in days j = 1, …, 391: and , in which and represent normally distributed constant measurement errors.

Effective reproductive number and attack rates

When individuals are homogeneous and mix uniformly, the effective reproductive ratio R_eff(t) is defined as the mean number of infections generated during the infectious period of a single infectious case at time t. In this model, the effective reproductive ratio can be written as a function of model parameters (see Section 3 of the Web Supplementary Material for details), namely When neglecting the deaths, the proportion of infected individuals assuming no waning immunity – also called attack rates – among the population in each region at a given date is given by:

Framing the problem

In order to sustain the modeling choices behind the disease transmission rate b, and more generally the selected inference procedure, we used Kalman-based estimation strategies Simon (2006); Bensoussan (2018) applied to the epidemic model (1). We define the global transmission rate t ↦ b(t), which accounts for the proportion of susceptible removed from the system thanks to vaccination: We propose to estimate b and ultimately to retrieve b using an estimation of vaccinated people V, see Section 2.2. In each region, we then introduce a dynamic equation for b of the form where νⁱ consists in a Wiener process such as for all t, s ≥ 0, νⁱ(t) − νⁱ(s) ∼ 𝒩 (0, (t − s)σ_ν) with is known and constant in all regions and gⁱ is a function describing the evolution of the global transmission rates bⁱ in each region.

After discretization using forward Euler time-scheme with small-enough time-step δt, we end up with a discrete-time dynamical system applied to the variable x = (E, I, R, A, H)^T ∈ ℝ⁵, for each region i = 1, …, 12: where f accounts for the dynamics of (E, I, R, A, H) in (1) while S is reconstructed after-wards using S = N − (E + I + R + A + H) in each region. In the time discrete system, now represent independent random variables, normally distributed with 0 mean and a variance equal to . Moreover, if constant parameters have to be estimated, the vector gathers all of them. For the estimation in the following, we transform the variable z to account for biological constraints. Remarking that all state variables are positive and bounded by N, the total population size of the region, we transform the state variable using x ↦ logit(x/N). We also apply a similar transformation to b ↦ logit(b/max_b). Of note, state and transmission rate variables are more likely to have Gaussian distributions in the transformation space. After calibration, we fixed max_b = 1.5 and checked that other values did not substantially changed the results (result not shown).

Population approach

To perform the estimation, we rely on an extension of the classical Unscented Kalman filter (UKF) Julier and Uhlmann (1997, 2002); Simon (2006). The particularity in this application, as described in Section 2.4, is that multiple series of data are observed jointly in multiple regions as we observe multiple realizations of the same epidemic. In order to take into account in our Kalman estimation that parameters in different regions are correlated as in a population approach, we follow a recently proposed population-based Kalman formulation Collin, Prague and Moireau (2020). As in mixed-effect models Lavielle (2014), each initial uncertainty variable is assumed to be randomly distributed around a common population intercept with a Gaussian distribution of unknown covariance Q₀, namely: When we treat the population intercept as the empirical mean over the population members in the construction of the objective function, we recover a classical filtering problem Simon (2006) on the aggregated variable . The only difference is the formulation of the initial covariance prior which couples the observations across regions and can be written as: where ⊗ indicates a Kronecker product, is a prior of Q₀ and M a small penalization matrix that guarantees the overall invertibility of . As a consequence, the matrix is not block diagonal with respect to the region i and thus all the region dynamics are coupled. The resulting time-discrete Kalman estimator couples all the regions together to produce a population-based estimation. Note that in such strategy, forcing a variable to be constant in the population is possible by simply choosing such that is very small with respect to Tr(M). Conversely, Tr(M) small with respect to in a large population N_r »1 will encourage each region to remain decoupled from each others. Starting from a given prior knowledge, our Kalman implementation uses the available measurements , to compute recursively in time the following estimates: and Then exploiting that gather the augmented state , 1 ≤ i ≤ 12 of the 12 regions, a simple post-processing over the regions provides estimations of b and the state variable. We refer to Collin, Prague and Moireau (2020) for further details about this Kalman-based population approach.

Estimation strategy

As we want to inject as few information as possible on the shape of the transmission rate b, we will first consider that the Wiener process ν defined in Equation (3) is a time-dependent function which thus encompass the complete dynamics of b. This corresponds to choosing g ≡ 0, meaning we have no prior knowledge on the evolution of transmission rates, i.e. on the NPI effects. However, to avoid overfitting, our objective is to disentangle the latent trajectory of b from possible noise in the data which can be observed in Figure 1. We adopt a 3-step approach described below consisting in smoothing the trajectory of b:

1. Estimate a reasonable prior for the initial transmission rate bbefore the start of any NPI. We use data before the first lockdown (10 days available) and assume the transmission rate for t = 1 …10 days to be a constant. In other words, we apply the population Kalman filter estimation described above with θ from Equation (4) being reduced to b_init.

2. Estimate the shape of bwith prior on initial value but without any prior on the dynamics. We set the initial value of bⁱ for times t = 1 …10 days to and take g(t) ≡ 0 such that Equation (3) rewrites db(t) = dν(t). We apply the population Kalman filter estimation described above. Of note, the parameter vector θ from Equation (4) is now empty. As there is no information on b, the model error is very important and could lead to an over-fitting of the data. We then build a prior for the dynamics of b by fitting a parametric shape based on sum of logistic functions on the weighted average trajectories of b over all regions with least-square method. Then, we build a prior for the dynamics of b by fitting a parametric shape based on sum of logistic functions on the weighted average trajectories of b over all regions with least-square method. We choose logistic functions as they seem well adapted to model variations due to for example stiff lockdowns or smooth unlocks. We set the number of logistic functions from the number of main changes of variations of the weighted average trajectories of b. This function will represent a prior g on the dynamics.

3. Estimate the shape of bwith prior on initial value and informative prior on the dynamics. We set the initial value of bⁱ for times t = 1 …10 days to and take db(t) = g(t)dt + dν(t) such as in Equation (3). We apply the population Kalman filter estimation described above. As the dynamics of b still includes a modeling noise, the shape of b could be different from the prior transmission rate defined from Step 2. This is the result of this final analysis that is used for further description of the effect of NPI.

Mixed effects model

Using the function b obtained with population Kalman filter as described in Section 2.4, one can determine the impacts of NPIs, seasonal weather conditions, and VoCs on the transmission rate . We followed the example of Flaxman et al. Flaxman et al. (2020) - in considering interventions effects to be multiplicative. There-fore, we fitted a linear mixed effect model on the logarithm transformation of the transmission. The model equations are given in Section 4.2 of the Web Supplementary material. It consists of an intercept, which represents the average transmission of COVID-19 across all regions without NPI without VoC circulation and for an average French weather condition (W=0), the sum of effect of the 10 NPIs described in Section 2.1, the effect of weather and the effect of the VoC percentage. It is known that the transmission of SARS-CoV-2 is very different indoors and outdoors Bulfone et al. (2021). Thus, we also added an interaction effect between bar and restaurant closure and the weather, accounting for the opening of outdoor sitting areas. Of note, this is the only interaction tested to avoid overfitting. We finally added random effects to account for heterogeneity between regions. We added a random intercept and random slopes for the effect of first, second lockdown and 6PM curfew as they may be highly variable between regions. We assume a full covariance matrix for the random effects, allowing effects to be correlated together in each region. In particular, we believe that they may be impacted by various factors not accounted for in the epidemic model, such as the density, the age distribution or the urbanization of the region.

Note that we added a 7-day delay to the lockdowns, which can be interpreted as a necessary time to allow people to organize themselves for adaptation (implementation of working from home, childcare etc …). This choice has been motivated by the observed 7-day delay in the transmission rates obtained by Kalman filters and will be discussed. To ease interpretation of the estimated effects, parameters were transformed back and expressed as a decrease or an increase in transmission in percentage by applying the function x ↦ 100(e^x − 1). Classical 95% confidence intervals were obtained using 100(e^{x+/−1.96SE(x)} −1), were SE is the standard error obtained from the regression.

Interpretation of the impact of seasonal weather conditions

To facilitate understanding, we calculated some NPI effects during the summer and winter periods. To do so, we replace in the results W by its summer (resp. winter) average from June 21^st, 2020 to October 21^st, 2020 (resp. before June 21^st, 2020 and after October 21^st, 2020) over all regions.

Basic reproductive number

The intercept of the regression presented above represents the mean transmission rate over all regions when there is no NPI in place, no VoC and the weather condition is taken to be the average weather condition over a year in France. Thus, inserting it in Equation (2) directly provides the basic reproductive number and 95% confidence intervals can be computed using the standard error of this parameter.

Basic reproductive number

After adjusting for average weather over a period of one year in France, we estimate the national average basic reproductive number at 3.10 [2.95 ; 3.26].

Effects of NPI on transmission rate

Table 2 summarizes the estimation of the effect of NPI on the transmission rates derived from the fixed effects of the model provided in Section 2.5. Regression residuals, fixed and random effects are given in Section 4.2 of the Web Supplementary Materials. Figure 7 shows the corresponding individual fits. The very restrictive first lockdown (respectively the less restrictive second lockdown) decreases the transmission rate by ∼ 80% (respectively ∼ 55%). The impact of school closure is around 7% and the bar and restaurant closure around 10%. The barrier gestures decrease the transmission by ∼ 46%.

• Download figure
• Open in new tab

Fig 7.

Transmission rates. Estimated using the population-based Kalman filter. Individual fits for each region. Random effects on log(b_init), b_l1, b_l2 and b_{curf 6PM}.

One can see that the values between the curfew at 6PM and the curfew at 8PM (around 30%) are very close and not statistically different (3% [-2% ; 8%]). We demonstrate that all NPI in this analysis decrease the transmission and have an effect significantly different from zero.

Effect of weather

The seasonal weather conditions variable is unitless, so the interpretation of its estimated effect should be made by comparing a value to a reference, which is set to be the average weather in France. We find that transmission is on average significantly increased by 10% during the winter period, and significantly decreased by 22% during the summer period, compared to average weather. Figure 8 (top) presents the estimated effect of seasonal weather conditions over the period of study, by region.

• Download figure
• Open in new tab

Fig 8.

Top: Estimated effect of seasonal weather conditions and its 95% confidence band for the 12 regions of interest over the study period using the global average value over the period as reference for comparison. Middle: Estimated effect of bar and restaurant closure for the 12 regions of interest over the study period. In red, the main effect of -10%. In black and grey, the effect with the interaction with weather conditions and its 95% confidence band. Bottom: Estimated effect of VoC appearance and its 95% confidence band from January 1, 2021.

The estimated interaction between the closure of bars and restaurants and the seasonal weather conditions is statistically significant (p=0.037) but is also complex to interpret. We show that although closing bars and restaurants always has a significant effect, decreasing the transmission, it is slightly more effective in winter (11% [8%;14%] decrease) rather than in summer (8% [4%;11%] decrease). Figure 8-Middle presents the estimated effect of bar and restaurant closure for the 12 regions of interest over the study period and its interaction with weather conditions. The weather conditions impacts the estimation, with a stronger effect on transmission reduction in Northern regions compared to Southern ones, and also a stronger effect is observed in the Spring 2021 due to more favorable weather.

Effect of VoC

Concerning the effect of the three VoCs, our results show that an increase from 0 to 100% of the proportion would increase the transmission rate by an estimated 22% [15% ; 28%]. Figure 8 (bottom) presents how the increased transmissibility varies over time due to VoC proportion.