key: cord-0450925-lyv9y8g7
authors: Borgonovo, Emanuele; Lu, Xuefei
title: Is Time to Intervention in the COVID-19 Outbreak Really Important? A Global Sensitivity Analysis Approach
date: 2020-05-04
journal: nan
DOI: nan
sha: 3f13edf13cdb5709fadd85c3c420567edae637b7
doc_id: 450925
cord_uid: lyv9y8g7

Italy has been one of the first countries timewise strongly impacted by the COVID-19 pandemic. The adoption of social distancing and heavy lockdown measures is posing a heavy burden on the population and the economy. The timing of the measures has crucial policy-making implications. Using publicly available data for the pandemic progression in Italy, we quantitatively assess the effect of the intervention time on the pandemic expansion, with a methodology that combines a generalized susceptible-exposed-infectious-recovered (SEIR) model together with statistical learning methods. The modeling shows that the lockdown has strongly deviated the pandemic trajectory in Italy. However, the difference between the forecasts and real data up to 20 April 2020 can be explained only by the existence of a time lag between the actual issuance date and the full effect of the measures. To understand the relative importance of intervention with respect to other factors, a thorough uncertainty quantification of the model predictions is performed. Global sensitivity indices show that the the time of intervention is 4 times more relevant than quarantine, and eight times more important than intrinsic features of the pandemic such as protection and infection rates. The relevance of their interactions is also quantified and studied.

The COVID-19 epidemic initially developed in the Chinese region of Wuhan [1] before becoming a pandemic. As of public records on 21 April 2020, 2,478,634 people have been affected internationally, and more than 26% of the world population is involved in some form of shelter-in-place or lockdown order. In the public political debate and in the scientific debate as well, great attention is posed on the effectiveness of intervention measures. Time-wise, Italy has been one of the first Countries, after China and Iran, significantly hit by the pandemic and then forced to resort to the implementation of strictly increasing quarantine measures by an exponential increase in the COVID-19 infection. The situation in Italy up to the first half of March is well described in [2, 3] . Embracing the recommendations of the scientific experts (see also the considerations in [2] concerning the exponential growth of the epidemic at the beginning of March 2020), since March 08 and March 09 2020 the local (Lombardy) and national (Italy) governments have imposed severe containment measures (lockdown, for simplicity, henceforth). The measures have, indeed, impacted the progression of the epidemic and avoided the exponential growth. On 20 April 2020 Italy counts a total of 181,228 confirmed individuals, with 24,114 deaths and a number of currently infected individuals equal to 108,237. [4] poses the question of how country-based mitigation measures may influence the course of the COVID-19 epidemic world wide. We have seen alternative attitudes internationally, with some countries intervening with strict lockdown measures within a short amount of time since the insurgence of the epidemics, and other adopting a more stage-phased approach.

Questions that emerge are, then, whether the effect of the intervention can be quantitatively observed from the data and whether its effect is immediate or delayed. In fact, it is natural to expect a time-lag between the moment an intervention is decided and the time in which it is actually taking in place. A the same time, the intervention time is not be the sole factor explaining a change in the pandemic time progression. The quarantine, the infection rate, the protection rate, the average latent time are factors determining an epidemic evolution that are expected to play a role. We then wish to quantify the relative importance of the intervention time with respect to these other factors.

To answer these questions, we resort to mathematical modelling coupled with global sensitivity analysis. Forecasting the expansion of the COVID-19 epidemics is subject of intensive research investigation [5, 6] . We investigate the behavior of the pandemic combining predictions from a generalized SEIR model recently used in [7, 8, 9] to describe the COVID-19 pandemic in China with machine learning approaches to increase interpretability of the findings [10] . We fit the generalized SEIR model to data of the pandemic in Italy covering the period 24 February to 20 April. We then run an uncertainty analysis with 1,000,000 Monte Carlo simulations and apply global sensitivity measures to quantify the relative importance of the factors mentioned above. We employ a replicated finite-difference decomposition method to study the intensity and sign of factor interactions.

Time Delay of Intervention Effects. We start with results for the effects of the intervention time. 

React on 09-Mar Figure 1 shows the following. The graphs represents the dynamic evolution of the pandemic comparing real data (circles) against the generalized SEIR model forecast (continuous line). The first graph report the results obtained for the case in which the intervention is modeled as fully active in the first day. We register a poor fit (R 2 = 0.3861), with the forecast largely underestimating the number of totally infected individuals (the root mean squared error (RMSE) is RM SE = 1.7416 × 10 4 ). The second graph shows reports the comparison in the case the intervention is assumed to have effect 1 day later than the official issuance date. We register better fit, confirmed by an increase of R 2 from 0.3861 to 0.5726 (the RMSE is now 1.4611 × 10 4 , with a decrease of 16%). Then, graphs (c) to (f) show the effect of simulating the intervention with one additional day of delay per graph. Note that by simulating the intervention as starting 5 days later than the actual issuance date, the fitting performance improves notably. We register R 2 = 0.9912 and RM SE = 2.1365 × 10 3 , with a decrease of 87.73% with respect to graph (a).

Due to the fact that intervention time is not the only factor influencing the forecast, we also tried to tune the forecast by changing other parameters such as I 0 and δ and combinations, while leaving the intervention time at the issuance date. Our search did not lead to satisfactory results, signaling the intervention time as the main source of discrepancy. To confirm the relevance of the intervention time, we performed further analyses via global sensitivity measures. Figure 3 shows that time to intervention (Interv.) is the key-driver in the variability of model predictions according to all the sensitivity measures used, i.e., variance-based total indices (T i ), first order indices (S i ) and the distribution-based sensitivity measure (β Ku ) (see the supplementary quantitative appendix A for further details on the definition of these indices and the motivation for their simultaneous utilization.) It is followed by quarantine rate, number of initially infected individuals (I 0 ) and infection rate (β). On a relative scale time to intervention is (depending on the sensitivity measure) 4.70 to 6.18 times more important than the quarantine rate. In turn, the quarantine rate is 1.50 to 2.20 times more important than the initial number of affected individuals, which, in turn, is 1.30 to 1.60 times more important than the infection rate. Protection rate, infection rate and average latent time play an even minor role. Figure 4 displays the conditional regression plots estimated from a subsample of size 500,000. Interaction Quantification Figure 5 displays the magnitude of the interaction effects among the factors. The highest recorded interactions (the first six bars are individual effects) are among δ, I 0 and the intervention time, namely the three most important factors. In terms of overall relevance, interaction effects account for less than 8% of the variance of the affected individuals. The mean dimension is 1.1683 (see Appendix A for the mathematical definitions). This means interactions have a low overall effect, although the effect is not entirely negligible. In fact, the highest interaction (between intervention time and δ) is 17.86% of the highest individual effect. Aside this interaction, we register a non-negligible interaction between intervention time and initial number of infected individuals (I 0 ). Note that these three factors are also the ones identified as most influential by the global sensitivity measures.

We have made inference on publicly available data of the COVID-19 pandemics in Italy through application of a SEIR model and the use of statistical techniques for sensitivity analysis. The investigation has evidenced a time-lag between the moment interventions have been issued and the moment in which their effect has actually taken place changing the pandemic path. Indeed this result seem to reflect a physical mechanism. While shelter-in-place measures are officially established on a given date, it is unlikely that they have an immediate effect. For instance, individuals might not be able to interrupt all contacts immediately, they do not have instantaneous access to masks or other protecting equipment. Moreover, some activities cannot shutdown immediately due to technical reasons, and certain industrial compartments have to remain open to keep delivering essential services to the population. ). The first six bars concern the individual factor effects, the next 15 represent the average magnitudes of the second order effects, the next 20 bars refer to third order effects, etc.. factors relate to social distancing measures. Not only, but combined with the results of the time to intervention analysis, the findings indicate that the earlier intervention measures are introduced, but also the more rapidly (strongly) they can be implemented, the more effective they are in flattening the pandemic curve.

From a modeling viewpoint, the strong impact of the intervention time suggests indeed that this variable needs to be taken into account to increase the prediction accuracy of SEIR models.

Regarding limitations of the present analysis, further work can be done to determine the distributions of the factors more precisely and, potentially, to include statistical correlations among the factors. Also, the analysis has been focused on the total number of infected individuals at 56 days from the initial observations. However, the analysis can also be made time dependent, varying the final time (e.g., using an analysis at 50 days or at 100 days). Of course, the longer the time horizon, the larger the uncertainty in the forecast we expect. Also, we have chosen one specification of the SEIR models and alternative specifications may be used.

The paper moves in the direction advocated by [10] . The uncertainty quantification is directly applicable to other SIER models and/or even to other type of models used to forecast or study the COVID-19 pandemics. Also, while we have analyzed data from Italy, the analysis can be replicated with data coming from other Countries and regions. The approach yields coherent insights not only regarding the most important inputs but also about the role of interactions and direction of change. This enriches the spectrum of insights gained from the model and the data with respect to current practice, with analyses that are limited only to varying one factor at a time. The analysis helps modelers in understanding key-factors of uncertainty increasing transparency. Collecting information on these factors might reduce uncertainty in the forecasts improving modelling efficacy and helping result communications to the public (see the recent critique in [14] ).

We consider a representative of the SEIR family of models. These models are a cornerstone in epidemiological studies [11, 12] and have been used in the very recent studies of [7, 8] regarding the COVID-19 pandemics. We consider the version used to simulate the COVID-19 outbreak in China in the work of [9] . This generalization allows a new quarantined state that includes the effect of soft measures before lockdown. In total, the following seven different states are considered: susceptible S(t), insusceptible P (t), exposed (but not yet infectious) E(t), infectious I(t), quarantined Q(t), recovered R(t), deceased D(t). The generalised SEIR model. Adapted from [9] .

The model includes the following parameters: i) protection rate: α; ii) infection rate, β; iii) average latent time: γ −1 ; iv) average quarantine rate: δ; v) time-dependent cure rate coefficients: λ 0 and λ 1 (the time dependent cure rate is written as by λ (t) = λ 0 (1 − exp (−λ 1 t))); vi) time-dependent mortality rate coefficients: κ 0 and κ 1 (the mortality rate curve is written as κ (t) = κ 0 exp (−κ 1 t)). For simplicity, the birth and natural death processes are not considered so that a constant population is assumed. The above model has been implemented in MATLAB [13] and can be openly accessed at https://github.com/ECheynet/SEIR.

We collected data from Italy's the Department of Civil Protection publicly available repository https://github.com/pcm-dpc/COVID-19. Figure 6 displays the number of total, quarantined, recovered, deceased individuals as a function of time, as per the publicly available data from 24 February till 20 April. The vertical lines display the dates of 09 March 2020, at which the lockdown measures were officially deployed. Fitting the generalized SEIR model discussed above on the data from 24 February to 08 March (13 days) for Italy, we obtain the pre-lockdown parameter estimates. Then, fitting the same model on the data from 09 March to 20 April (42 days), we obtain the after-lockdown parameter estimates reported in Table 1 . Note that the after-parameters estimates are significantly different from the pre-lockdown estimates. For instance, the protection rate α increases by over 2,000 times; The average latent time increases by over 6.5 times; the quarantine rate δ decreases by 37%. Uncertainty quantification To further understand the results, we have let the model undergo an uncertainty quantification in order to be able to compute global sensitivity measures. We used the parameter distribution assignments displayed in Table 2 (see also [9] ) and generated a random Monte Carlo sample of size N = 1, 000, 000. Table 2 : Uncertainty quantification: the after lockdown parametric values of α, β, (γ) −1 and δ are assigned distributions as per the data of [9] . The intervention time is assumed to be a uniformly distributed random variable in a one-week range after the actual issuance date. Independence is assumed.

The intervention day Distribution Discrete uniform Discrete uniform Support I 0 ± 20% {09-Mar+z}, with z = 0, 1, ..., 7

From this sample, we have quantified the uncertainty about the forecasted number of infected individuals in Italy from 24 February to 20 April (56 days). Also, we have calculated importance indices that go under the name of global sensitivity measures to determine the quantitatively the relative importance of the inputs [See Supplementary Appendix A]. Specifically, we have used the first and total order global sensitivity measures based on variance contribution and a sensitivity measures based on the distance between cumulative distribution functions. The rationale of the choice is as follows. Regarding importance, the total order as well as the distribution-based sensitivity measure possess the properties that they are null if and only if the model output is probabilistically independent of the input of interest. Thus, their null value reassures us that the input is, indeed, not influential. Regarding interaction quantification, the difference between the first and the total order indices provides information about the relevance of interactions, giving insights on whether a factor's influence is due to its individual action or due to interactions with other factors. This insight, however, is an overall insight and does not reveal which interactions (if any) are important. Towards this understanding we use a sample of size 20,000, replicating a full factorial design and decomposing the output changes via a finite difference operation (see Appendix A fpor further details). The scheme allows us to compute a total of 20,000 first order effects for each input, 20,000 second order interaction indices for all pairs, and 20,000 interaction indices of order 3, 4, and 5 and 6, enabling interactions to be thoroughly studied.

From the original sample, we also computed conditional regression curves. These curves express the expected number of infected individuals as a function of each of the uncertain parameters. The graphs lead information about whether the increase in a parameter leads to an increase (decrease) of the number of infected individuals.

This section is divided into two parts. In the first part, we address the definition of the global sensitivity measures used in this work. In the second part we ketch the principles of their estimation and list the corresponding subroutines.

For broad overviews on sensitivity analysis, we refer the reader to [15, 17] . We concisely overview the methods used here. Let y = g(x), where g : X → R is a multivariate input-output mapping, x = (x 1 , x 2 , ..., x n ) be the vector of the model inputs, with x ∈ X and y the output of interest. Let also (X , B(X ), P) be the measure space of the model inputs and let x 0 , x 1 ,..., x k ,x k+1 ,..., x N be a collection of points in X sampled from the input probability measure. Let y 0 = g(x 0 ), y 1 = g(x 1 ),...,y N = g(x N ) be the corresponding values of the model the output. Then, let ∆y k = y k+1 − y k denote the variation in the model output registered with the inputs shift from x k to x k+1 . We then have a sequence of finite changes ∆y k , k = 0, 1, 2, ..., N − 1. We can decompose each change in 2 n effects, writing:

where

In eq. (2), the point x k+1 second order contributions, quantifying the residual interaction between x i and x j ; a similar interpretation is shared by higher order terms. Dividing these second order effects by (x k+1

, we obtain the second order Newton ratios that provide indications about whether interactions are locally positive (synergistic) or negative (antagonistic).

Sampling the model N -times, we obtain N − 1 replicates of the finite change indices, φ k→k+1 i,j,...,m . These sensitivity measures can be used, on the one hand, to obtain detailed regional information on interactions. Assume the inputs are independent. Let V Y the variance of Y induced by uncertainty in the model inputs. As proven in [20] , we can write

where

... (4) and the functions g i (x i ), g i,j (x i , x j ) obtained by taking appropriate conditional expectations [20, 17] . In equation (3), V 2 i is the portion of the variance attributed to X i alone, V 2 i,j the portion of the variance attributed to the residual interaction between X i and X j . In the literature, one considers the normalized version

called first order Sobol' indices. One also defines the total sensitivity index of x i , here denoted by T i [21] . The index is written as:

and represents the total portion of the variance of Y associated with x i . We note that by a result due to [22] , the sum of the total order indices equals the mean dimension of g in the superimposition sense.

There is a bridge between finite change sensitivity indices φ k→k+1 i and the total index T 2 i . In fact, let Φ i denote the random variable whose realizations are the φ k→k+1 i . Then, from the theory of the Sobol' estimation design [23] [24] [25], we get

That is, by calculating the variance of the available first order effects and halving we obtain an estimate of the total order indices for each input. These estimates, compared with estimates of the first order indices provides information on how much each input is involved in interaction with the remaning ones. In particular, we note the notion of mean effective dimension provides a quantitative indication about the relevance of interactions. This notion has been introduced in [26, 22] . The intuition is as follows. The terms S i,j,...,r = V i,j,...,r V Y are positive and sum to unity. Thus, they can be regarded as forming a probability mass function. Actually, they place a mass on the random variable T whose realizations all the possible combinations of indices. Then, one defines the dimension distribution of g in the superimposition sense as the distribution of the cardinality of |T |. The mean dimension is then [22] :

respectively. A mean effective dimension in the superimposition sense equal to unity indicates the absence of interactions. [22] also proves that the mean effective dimension is equal to the sum of total effects. That is, we have:

T i . We recall that first order variance-based sensitivity measures fall in the so-called common rationale of global sensitivity measures [27] . That is, they are sensitivity measures of the type

where d(·, ·) is a distance or divergence between distributions, P Y is the marginal probability measure of Y and P Y |Xi is the conditional probability measure of Y given X i . Depending on the choice of d(·, ·), one obtains alternative ways of measuring importance. To illustrate, let F Y and F Y |Xi denote the cumulative distribution function of Y and the conditional distribution function of Y given X i . Selecting the Kuiper distance between cumulative distributions functions, we get [28] β Ku

We note that a null value of T i and β Ku i indicates that Y is independent of X i . First order variance-based sensitivity measure (S i ) do not possess this property and their null value is not necessarily a signal of the fact that Y does not depend on X i . However, first order variance-based sensitivity measures remain well defined when inputs are dependent, similarly to β Ku i , while problems associated with the interpretation of T i arise.

The selected sensitivity measures for our calculations are the first order main effects φ k→k+1 i , the higher order effects φ k→k+1 i,j,...,m computed at randomized locations. The implementation of the decomposition in eqs. (1) and (2) is implemented in the surboutine finitechanges.m, which is available at https://github.com/LuXuefei/ FiniteChanges. Then, such decomposition is performed at N randomized locations appropriately generated through random Monte Carlo generation implemented in the code rand.m and randsample.m, available at MATLAB.

From the randomized first order effects φ k→k+1 i , it is then possible to obtain estimates of the first order indices of all inputs from

The estimation of the sensitivity measures S i and β Ku i is based on the so-called given data approach. We refer to [29] and [27] for a detailed treatment. We just recall the key intuition here. Let X i a subset of the reals that represents the support of model input x i . Then, consider a partition of X i into M bins, such that ∪ M m=1 X m i = X i and X m i ∩ X r i = ∅ for r, m = 1, 2, ..., M , r = m. The intuition of a given-data estimation is then that of substituting the condition X i = x i (that is X i is exactly at x i ) with the bin condition X i ∈ X m i . Intuitively, this condition asks that X i is in a suitably defined interval around x i . With this intuition, one obtaines the estimate

where p m,i is equal to the probability that X i belongs to bin X m i , and where F Y and F Y |Xi∈X m i are the marginal empirical distribution function of Y and the conditional empirical distribution function given that X i ∈ X m i . Then, one links the partition size (M) to the sample size (N ), so that as N increases the parition size tends to zero and the limiting partition is X i itself. [27] prove that, under mild assumptions, (the continuity of the operator d(·, ·) in both its arguments and the Riemann-Stieltjes integrability of d(F Y , F Y |Xi )), ξ i converges to ξ i as the sample size increases; that is, ξ i is an asymptotically consistent estimator of ξ i . These calculations are implemented in the subroutine betaKS3.m. Finally, the conditional regression functions r i (x i ) = E[Y |X i = x i ] have been estimated using the subroutine cosi.m. Both subroutines are developed by Elmar Plischke and available at https://zenodo.org/record/885332#.XpnB44gzZGM (see [30] and [31] for details).

A new Coronavirus Associated with Human Respiratory Disease in China

COVID-19 and Italy: what next? The Lancet

Critical Care Utilization for the COVID-19 Outbreak in Lombardy, Italy: Early experience and Forecast during An Emergency Response

How will country-based mitigation measures influence the course of the COVID-19 epidemic? The Lancet

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

With COVID-19, modeling takes on life and death importance

S Flasche, and Centre for Mathematical Modelling of Infectious Diseases COVID-19 working group. Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The Lancet Infectious Diseases

Phase-adjusted estimation of the number of Coronavirus Disease

Epidemic analysis of COVID-19 in China by dynamical modeling

A Short Comment on Statistical versus Mathematical Modelling

A Contribution to the Mathematical Theory of Epidemics

Mathematics of infectious diseases

Generalized SEIR Epidemic Model (fitting and computation

What 5 Coronavirus Models Say the Next Month Will Look Like. The New York Times

Global Sensitivity Analysis -The Primer

Introduction to Uncertainty Quantification

Sensitivity Analysis: An Introduction for the Management Scientist

High dimensional model representations

Sensitivity analysis with finite changes: An application to modified EOQ models

The Jackknife Estimate of Variance

Importance Measures in Global Sensitivity Analysis of Nonlinear Models

The Dimension Distribution and Quadrature Test Functions

Sensitivity Estimates for Nonlinear Mathematical Models

Statistical inference for Sobol pick-freeze Monte Carlo method

Rabitti Simulator Inputs Screening: from Elementary Effects to Mean Dimensions

Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension

A Common Rationale for Global Sensitivity Measures and their Estimation

Invariant Probabilistic Sensitivity Analysis

Global Sensitivity Measures from Given Data

An effective algorithm for computing global sensitivity indices (EASI)

How to Compute Variance-Based Sensitivity Indicators with