Global sensitivity analysis in epidemiological modeling.

Operations researchers worldwide rely extensively on quantitative simulations to model alternative aspects of the COVID-19 pandemic. Proper uncertainty quantification and sensitivity analysis are fundamental to enrich the modeling process and communicate correctly informed insights to decision-makers. We develop a methodology to obtain insights on key uncertainty drivers, trend analysis and interaction quantification through an innovative combination of probabilistic sensitivity techniques and machine learning tools. We illustrate the approach by applying it to a representative of the family of susceptible-infectious-recovered (SIR) models recently used in the context of the COVID-19 pandemic. We focus on data of the early pandemic progression in Italy and the United States (the U.S.). We perform the analysis for both cases of correlated and uncorrelated inputs. Results show that quarantine rate and intervention time are the key uncertainty drivers, have opposite effects on the number of total infected individuals and are involved in the most relevant interactions.

Introduction

Research teams worldwide are actively involved in supporting decision-makers with forecasts on alternative aspects of the COVID-19 pandemic (Kaplan, 2020). A variety of modeling techniques, ranging from agent-based simulations to generalizations of the classical family of compartmental models are in use (Currie, Fowler, Kotiadis, Monks, Onggo, Robertson, et al., 2020, Robertson, 2019). Much attention is paid to models that support the assessment of the effectiveness of intervention measures (Ferretti et al., 2020). The type of interventions and their timing differ from country to country and, in some cases, from region to region within the same country.

The lack of knowledge about aspects such as the time-lag between when an intervention is issued and when it actually takes place, the effectiveness of quarantine measures and the epidemiological characteristics of the pandemic itself, creates uncertainty in model predictions. The recent manifesto of Saltelli et al. (2020) calls for the use of methods that allow researchers to perform proper uncertainty quantification and explain how output uncertainty can be apportioned to its sources.

At the same time, an investigation of the literature (see Section 2) shows that sensitivity analysis and uncertainty quantification are often performed in an ad-hoc fashion or even neglected: there is a notable gap between the methods actually employed in operations research (OR) pandemic modeling and the techniques made available by recent literature. For instance, the most widely employed methods fall into the category of one-at-a-time sensitivity methods, a finding in line with the investigation of Saltelli et al. (2019). However, if a model response deviates sharply from additivity, analysts need to be careful in attributing a general validity to results of one-at-a-time sensitivities, because their indications become strongly dependent on the point at which the inputs are fixed (Saltelli & Annoni, 2010). At the same time, scenario analysis is performed solely to evaluate the model at a few locations of the input space, and to draw conclusions from qualitative comparisons of the results, not accompanied by the simultaneous calculation of sensitivity measures. This makes the analysis far from a thorough uncertainty quantification. Overall, the lack of more systematic approaches may create problems in result interpretation and communication. Recent developments, however, show that it is possible to gain insights on input importance, trend analysis and interaction quantification when performing a scenario analysis and that a better use of one-at-a-time sensitivities can lead to indications about interactions.

In this work, we aim to bridge this gap, proposing a goal-oriented approach to support OR modeling in the context of the COVID-19 pandemic. We frame the analysis within three sensitivity settings: factor prioritization, trend analysis and interaction quantification, and use a combination of methods developed in the machine learning literature for explainable artificial intelligence and in the OR literature of global sensitivity analysis. Specifically, for factor prioritization, we rely on the combination of variance-based and moment independent sensitivity measures (Borgonovo, Plischke, 2016, Saltelli, Tarantola, 2002). This choice allows us to address both dependent as well as independent inputs. To compute total indices (one of the commonly used variance-based sensitivity indices), we employ Jansen intuition (Jansen, 1999) that allows their estimation as variance of one-at-a-time sensitivities — see also the recent work of Owen & Hoyt (2021). Knowledge of these sensitivity measures provides insights about the internal structure of the model, helping analysts to understand the overall relevance of interactions. However, these indices do not reveal the sign of interactions, nor do they indicate the group of inputs involved in the most relevant interactions on a local scale. Then, we perform further interaction analysis using the design of experiments for higher order interactions (Kleijnen, 2005, Wu, 2015). Overall, the design that we employ is a bridge between the intuition at the basis of Tornado diagrams (Eschenbach, 1992, Howard, 1988, Wang, Dyer, Hahn, 2017) and screening designs such as the Morris method (Morris, 1991). In particular, the availability of one-at-a-time sensitivities at randomized locations can be used to obtain insights on the direction of change of the simulator.

We propose to complement this information with individual conditional expectation (ICE) plots that provide global indications about how the quantity of interest varies when a given uncertain input varies over its support (Goldstein, Kapelner, Bleich, & Pitkin, 2015). We use these methods to examine a generalized SIR model introduced in Wang et al. (2020) and Peng, Yang, Zhang, Zhuge, & Hong (2020) to describe the early progression of the pandemic in China; the model is modified to accommodate the effect of interventions. We fit the model to data corresponding to the following time series: Italy from February 24th to April 20th, Spain from February 26th to April 22nd, Germany from March 2nd to April 27th, and the U.S. from March 4th to April 29th (The dates are chosen so that each spanned period contains 56 days). Following the framework of Berger & Smith (2019), we then consider uncertainty epidemiological as well as intervention related inputs; we run 1,000,000 Monte Carlo simulations to assess the corresponding uncertainty in model predictions. We then exploit these model runs to compute sensitivity measures that provide insights for the above-mentioned settings.

Among other findings, it emerges that intervention time is the most important input in Italy but is less important in the U.S., where intervention policies were state-based and lacked a centralized intervention date for a national lockdown. Also, intervention-based inputs are involved in the most important interactions, with a prevailing negative interaction sign for the interaction between quarantine rate and intervention time (see Section 4 for greater detail). To test the stability of these indications, we repeat the experiments considering input correlations. We register a difference in factor prioritization, with the sum of first order variance-based indices greater than unity when the model is fitted to data for Italy and a switch in ranking, with one of the correlated variables becoming more relevant than in the uncorrelated case for both Italy and the U.S. Results concerning trend and interaction quantification remain virtually unchanged.

The reminder of the manuscript is organized as follows. Section 2 reviews relevant literature on sensitivity analysis in OR epidemiological modeling. Section 3 presents the sensitivity goals, and the methods we use. Section 4 presents the numerical test case, illustrating model calibration, data collection and sensitivity analysis results. Section 5 discusses the main findings of the work and highlights areas for future research. Supplementary Section 1 presents additional results for Spain and Germany.

Related literature: SIR models in operations research and sensitivity analysis

The Operations Research/Management Science (OR/MS) literature in pandemics is vast and has addressed aspects ranging from disaster control to optimal vaccination. We refer to works such as Altay & Green (2006), Dasaklis, Pappis, & Rachaniotis (2012) and Dimitrov & Meyers (2010) for comprehensive reviews. Without claiming exhaustiveness, in this section we provide a concise overview on models used to quantitatively support investigations, with a focus on how sensitivity analysis has been performed.

Models used to support OR/MS investigations range from differential equations/system dynamics (e.g., Kaplan, Craft, Wein, 2003, Rahmandad, Sterman, 2008)1 , to agent-based simulation (Pichler, Pangallo, del Rio-Chanona, Lafond, & Doyne Farmer, 2020), to network-based models (see the reviews of Dimitrov, Meyers, 2010, Robertson, 2019), to discrete event simulators, to hybrid simulators (Currie et al., 2020). In several applications, forecasts from mathematical epidemiological models are fed into optimization models, playing the role of simulating infectious disease transmission dynamics. Optimal policies are identified by building corresponding mathematical programs; the optimization objectives involve minimization of the vaccination cost (Yarmand, Ivy, Denton, & Lloyd, 2014), of the proportion of infected population (Du, Sai, Kong, 2020, Duijzer, van Jaarsveld, Dekker, 2018, Rachaniotis, Dasaklis, Pappis, 2012) or of the total number of fatalities (Ren, Órdoñez, & Wu, 2013), or maximization of both humanitarian and economic performance (Büyüktahtakın, Des-Bordes, Kıbış, 2018, Du, Sai, Kong, 2020, Torabi G, Ali Torabi, Shakouri G, 2018, Tebbens, Thompson, 2009).

SIR models (Kermack, McKendrick, 1991a, Kermack, McKendrick, 1991b, Kermack, McKendrick, 1991c) and their generalizations (generalized SIR, henceforth) are frequently used (for more details see Hethcote, 2000, Murray, 2003). For example, Kaplan, Craft, & Wein (2002) propose an early application of a generalized SIR model to support intervention policies against smallpox bioterror attacks. Berman, Gavious, & Menezes (2012) use a generalized SIR model to support optimal resource allocation in emergency response following a smallpox bioterrorist attack on an airport. Tebbens & Thompson (2009) analyze alternative policies for long-term disease eradication simulating the progression of a hypothetical pandemic with a two-disease SIR model. Rachaniotis et al. (2012) use a generalized SIR model to simulate the evolution of A(H1N1)v influenza in the context of determining the optimal scheduling of mobile medical resources across alternative geographical regions. Yarmand et al. (2014) rely on a generalized SIR model to support the determination of optimal 2-phase vaccine allocation strategies, and Enayati & Özaltın (2020) rely on a similar model for the minimization of vaccine doses in the context of an influenza pandemic. Shamsi G et al. (2018) simulate a pandemic progression using an SIR model to support the evaluation of an option contract for ensuring provisioning from two vaccine suppliers. Beyond the OR/MS applications, SIR models have been applied to describe the spreading of alternative diseases. Sharareh, Sabounchi, Sayama, & MacDonald (2016) and Büyüktahtakın et al. (2018) use generalized SIR models to support policies against the EBOLA pandemic, Rimbaud et al. (2018) in the context of the Sharka (a plant virus) infection, and Du et al. (2020) in the context of the cholera infection.

Coronavirus-type infections are modeled in Dye & Gay (2003) and Lipsitch et al. (2003), who apply generalized SIR models to forecast the spread of the 2003 coronavirus-severe acute respiratory syndrome (SARS). With regard to the COVID-19 pandemic, Fang, Nie, & Penny (2020), Lin et al. (2020), Kucharski et al. (2020) and Tang et al. (2020) use generalized SIR models to describe the early phases of its evolution.

An important step in the use of pandemic simulators both within and outside OR applications is model calibration. The model is usually calibrated comparing the values of a quantity of interest after data collection against the forecasts of the same quantity simulated by the model (see Dye, Gay, 2003, Lipsitch, Cohen, Cooper, Robins, Ma, James, et al., 2003 for the SARS outbreak, and Dehning, Zierenberg, Spitzner, Wibral, Neto, Wilczek, et al., 2020, Ferguson, Laydon, Nedjati-Gilani, Imai, Ainslie, & Baguelin, et al. for the COVID-19 pandemic). Target outcomes include the number of infected individuals (Yarmand et al., 2014), the number of deceased (Ferguson, Laydon, Nedjati-Gilani, Imai, Ainslie, & Baguelin, et al., Lourenço, Paton, Ghafari, Kraemer, Thompson, & Simmonds, et al.), the basic reproduction number () (Berger, Herkenhoff, Mongey, 2020, Tanner, Sattenspiel, Ntaimo, 2008) or the epidemic peak (Tang, Wang, Li, Bragazzi, Tang, Xiao, et al., 2020, Yarmand, Ivy, Denton, Lloyd, 2014). The process involves the maximization of a performance measure or the minimization of a loss function, such as the coefficient of determination () (Rahmandad & Sterman, 2008), the root mean square error (RMSE) (Chang et al., 2021), the mean absolute error (MAE) (Fang et al., 2020) or the Akaike information criterion (AIC) (Bertozzi, Franco, Mohler, Short, & Sledge, 2020). In OR applications, we also register two additional (and opposite) situations: on the one hand, works in which the calibration is implicit, because researchers borrow models already calibrated in the literature and use these models to perform further investigations — see Enayati & Özaltın (2020); Kaplan, Craft, Wein, 2002, Kaplan, Craft, Wein, 2003; Lin et al. (2020); Ren et al. (2013); Tebbens & Thompson (2009). On the other hand, works in which inputs are elicited are either directly assigned by the modelers or derived from experts’ opinions, as in Shamsi G et al. (2018). This happens especially when models involve poorly known compartments, e.g., a network compartment in Rahmandad & Sterman (2008), or a mobile medical team in Rachaniotis et al. (2012). Even when carried out for models with known compartments, the calibration process does not eliminate uncertainty in the model forecasts: Not all the inputs that enter into the model may be subjected to calibration, on other inputs we might face problems of lack or sparsity of data. This leads to the need for a proper uncertainty quantification. At the same time, a thorough sensitivity analysis is needed to allow analysts to appreciate the key drivers of variability and how the forecasts depend on the uncertain inputs. In the remainder of this section, we investigate how sensitivity and uncertainty analyses are performed in the OR-related literature.

Regarding sensitivity analysis, the three most widely applied methods are one-at-a-time input variations, one-way sensitivity functions and scenario analysis. One-at-a-time sensitivity analysis investigates the changes in the model output when analysts vary one input at a time from a base case to one or more sensitivity cases. This is the technique at the basis of Tornado diagrams (Howard, 1988). It is used in Ren et al. (2013) to study changes in the number of fatalities under an optimal resource allocation strategy. The inputs (e.g., resource level and initial transmission rate) are varied from a low to a medium to a high value. A similar approach is used in Büyüktahtakın et al. (2018) where the sensitivity of the objective function at the optimum is tested for individual variations of migration and intervention related inputs. Similarly, Berman et al. (2012); Bertozzi et al. (2020); Lin et al. (2020); Wu, Leung, & Leung (2020) vary one-at-a-time a set of epidemiological input parameters to test their impact on the pandemic evolution.

A second popular approach is based on one-way sensitivity functions (van der Gaag, Renooij, & Coupe, 2007). Kaplan, Craft, Wein, 2002, Kaplan, Craft, Wein, 2003 study the marginal behavior of the number of fatalities varying five model inputs, that comprise the number of initially infected individuals and of available vaccinators, over predetermined ranges. Shamsi G et al. (2018) use several one-way sensitivity plots to study the variation of several outputs, such as the social cost and percentage of recovered individuals, responding to changes in vaccine delivery time. Tebbens & Thompson (2009) use one-way sensitivity functions to study the impact of population size and population structure on policy decision rules, such as the incremental cost-effectiveness ratio and the time until disease elimination.

While in a one-at-a-time sensitivity approach all input parameters are fixed but one, in a scenario analysis one can inspect the model response at several locations of the model input space. This technique is employed in Duijzer et al. (2018), who compare the value of the objective function (number of infected individuals) at the optimum varying vaccination time and vaccine efficacy rates simultaneously. Ferguson et al. (2020) discuss how peak hospital intensive care units (ICU) demand varies considering alternative combinations of mitigation strategies and policy triggers. Enayati & Özaltın (2020) study the changes in optimal vaccination policies varying simultaneously the number of geographic regions and vaccine efficacy. Lourenço et al. (2020) evaluate how the epidemic forecast changes by varying simultaneously the proportion of vulnerable population and the basic reproduction numbers. A similar approach is used in Bertozzi et al. (2020) and Tang et al. (2020).

A minority of works employ more complex approaches. Differential indicators (elasticities) are used in Upadhyay, Pal, Kumari, & Roy (2019), and the expected value of perfect information in Yarmand et al. (2014). A few works employ a global sensitivity analysis: Nsoesie, Beckman, & Marathe (2012) perform uncertainty quantification followed by calculation of the Pearson correlation coefficients between the model output and the model inputs; Sobol’ variance-based sensitivity indices are more recently applied by Rimbaud et al. (2018) and Alexanderian, Gremaud, & Smith (2020) to identify the most influential model inputs.

Regarding goals, we find the following. In the majority of works, the analysts’ goal is factor prioritization, i.e., the identification of the most influential input (Rachaniotis, Dasaklis, Pappis, 2012, Upadhyay, Pal, Kumari, Roy, 2019). A further frequently addressed objective is the stability of the model response, especially in the context of model calibration. To illustrate, Kaplan et al. (2003) compare one-way sensitivity functions obtained by the exact and the approximated models to check if the approximation accurately mimics the original mechanism; Enayati & Özaltın (2020) perform scenario analysis to calibrate their mathematical programming model with respect to geographic areas and age groups, and to validate the choice of the numerical optimization algorithm. Almost no work considers the relevance of interactions and, when sought, information on direction of change is mixed with factor prioritization.

Overall, this analysis reveals that sensitivity analysis and uncertainty quantification are most often neglected or performed unsystematically in the context of OR epidemic modeling, with the risk of communicating partial (if not misleading) insights about the model behavior and of under-exploiting the modeling efforts. The reminder of this work aims at reducing the gap. We start by proposing a goal-oriented approach.

A goal-oriented approach

This section illustrates the goal-oriented approach we introduce in this work. We refer to Appendix A for the mathematical definitions of the sensitivity measures we employ (Table 1 ) and additional quantitative details.

Table 1

Summary of SA insights and sensitivity methods used in this work.

SA insights/name	Symbol (Equation)	Dependent Inputsa
Prioritization
First order indices	S_i(6), ηi2(8)	Not affected
Total indices	T_i(9)	Interpretation affected
Kuiper-based	βiKu(24)	Not affected
Trend Analysis
Partial dep. func.	h_i(x_i)(13)	Not affected
One-way sensitivity functions	g_i(x_i;x_−i)	Not affected
First order Newton Quotients	νik→k+1(21)	Not affected
Interaction Quantification
Mean Dimension	D_g(12)	Not applicable
Higher order Newton quotients	νi,jk→k+1(22)	Not affected

Note: when inputs are independent, all above sensitivity methods are well-defined.

We consider that the analyst is dealing with a black-box input-output mapping , whose output depends on uncertain inputs , . In sensitivity analysis, the extraction of insights has been made formal by the notion of sensitivity analysis setting (Saltelli & Tarantola, 2002). A setting is a way of framing the sensitivity quest in such a way that the answer can be confidently entrusted to a well-identified measure (Saltelli et al., 2008, p. 24). (Henceforth, we shall use the terms goal/setting interchangeably.)

The first goal we address is factor prioritization. In a factor prioritization analysts are interested in determining the inputs that most contribute to the variability of the output. For this setting, we would suggest relying on an ensemble of sensitivity measures, because alternative sensitivity measures consider different aspects of the output response, and one needs to match the quantity of interest with the sensitivity measure (see Borgonovo, Hazen, Jose, & Plischke, 2021b, among others). Moreover, alternative sensitivity measures possess different properties (see Appendix A for details). An ideal ensemble could be the triplet composed of first order variance-based sensitivity indices, in Table 1, the total order sensitivity indices, , and a moment independent sensitivity measure. (In our numerical experiments, we shall rely on the index based on the Kuiper distance between distributions.) The rationale is that ranks inputs based on their individual contributions to the output variance, delivers importance in consideration also of the input interaction with other inputs, and ranks inputs based on their influence on the entire distribution of the output and complies with Renyi’s postulate D of measures of statistical dependence (Renyi, 1959): a null value of implies that the model output is independent of . A further motivation for the use of the total indices comes from their link with Tornado diagrams, one of the earliest tools developed for determining factor importance. Tornado diagrams are a graphical representation of one-at-a-time sensitivity measures. In the literature review, we have seen that operations research investigators frequently evaluate models fixing the inputs at pre-selected scenarios; for simplicity, consider a base case and a sensitivity case . Then, the finite-change indices represent the effect of varying input from the base to the sensitivity case. These main effects are the sensitivity measures of Tornado diagrams — The symbol denotes the set of all inputs excluding (see Appendix A for formal definitions). Thus, using Tornado diagrams the analyst ranks before if . However, this statement holds only for the chosen scenarios and , unless these two scenarios encompass all the possible realizations of the inputs in line with the critique in Saltelli & D’Hombres (2010). A natural way to remedy this shortcoming is to compute these sensitivity measures at several locations in the input space. Taking the variance of these finite change sensitivity measures allows one to obtain the structural contribution of the inputs to the variance of the model output (see Appendix A for further details).

A second insight frequently sought by analysts to increase transparency of black-box input-output mapping is the partial behavior of the quantity of interest as a function of one or more inputs (trend analysis, henceforth). Recently, several studies have appeared on this insight in machine learning, especially in relationship with interpretability (Guidotti et al., 2018). In terms of sensitivity analysis settings, we can formulate the goal of these investigations as: In trend analysis, analysts are interested in determining how the model output behaves as a function of one or more inputs . A first candidate sensitivity measure for this setting is the Newton quotient (Eq. (21)). In fact, if is monotonic in then will either be null or will consistently bear the same sign (positive if increasing, negative otherwise). Then, if the analyst registers a change in sign of the on a set of non-null measures, she can conclude that the model is not monotonic with respect to . Note that knowledge of immediately yields the corresponding Newton quotient . Thus, the evaluation of finite change indices at randomized locations provides a natural bridge between factor prioritization and trend identification. Recently, however, several indicators of trend have been developed in the machine learning literature. Particularly successful have been indicators such as partial dependence functions (Friedman, 2001), individual conditional expectation (ICE) plots (Goldstein et al., 2015) and accumulated local effect plots (Apley & Zhu, 2020). In this work, we shall make use of the so-called ICE partial dependence plots, which consist of graphical representation of a partial dependence function together with one-way sensitivity functions. We recall that the graph of a partial dependence function ( in Table 1) represents how the model output varies in expectation as the uncertain input of interest varies from its lowest to its highest value. Each one-way sensitivity function ( in Table 1) visualizes how changes as a function of given that the other inputs are fixed at a particular value. We recall that one-way sensitivity functions are widely used in the management sciences, where they are at the basis of graphical representations such as spider and spaghetti plots (Castillo, Gutierrez, Hadi, 1996, Eschenbach, 1992). One advantage of partial dependence as well as one-way sensitivity functions is represented by the fact that they are monotonicity consistent independently of whether inputs are correlated. Monotonicity consistency implies that if the input-output mapping is monotonic in , then the graphs of any one-way sensitivity function and of the partial dependence function are monotonic (Borgonovo, Baucell, Plischke, Barr, & Rabitz, 2021a). The implication is that the missed monotonicity of any of these graphs immediately indicates that the input-output mapping is non-monotonic (Appendix A report further details).

A third setting is interaction quantification. We can formulate this setting as follows: In an interaction quantification analysts are interested in determining whether interactions play an important role in the model response, what are the most relevant interactions, and whether interactions are synergic (positive) or antagonistic (negative). To gain information on item , natural candidate sensitivity measures are the differences , which, under independence, suggest the portion of due to interactions. We also recall that insights on the presence of interactions can be obtained by ICE and partial dependence plots (see Appendix A). In fact, it can be proven that if a model is separately additive in , then the one-way sensitivity functions are non-intersecting (they all differ by a constant). Thus, if the graphs of alternative one-way sensitivity functions intersect, the model is not additive. These indicators, however, do not reveal specifically what the most important interactions are, nor the sign of interactions. To obtain this information, one needs to compute higher order interaction indices. With higher order global sensitivity indices, analysts gain insights on the magnitude of higher order interactions, while with second order Newton quotients, (Table 1), they also gain information about their sign.

The presence of input dependences may affect the interpretation/properties of a sensitivity measure. The third column of Table 1 presents a synthetic summary for each method. The first order indices and are not impacted in their interpretation as expected reduction in the model output variance for fixing . The total order indices computed as variance of first order finite change indices modify their interpretation under input dependence (see Appendix A for further details). Moment independent sensitivity measures are not impacted, neither in terms of their interpretation nor mathematical definition. None of the selected indicators for trend identification are affected by the presence of dependences: The sign of is unchanged, and ICE plots as well as partial dependence functions are monotonicity consistent, as proven recently in Borgonovo et al. (2021a). In interaction quantification, the indication provided by (the not-normalized version of) the total indices estimated as variance of the first order finite change indices is unaffected by the presence of interactions (see Appendix A). The indications of second order Newton quotients are not affected by the presence of dependences. However under dependence, the differences are no more meaningful.

Let us now come to estimation. In our case, consistent and efficient estimators have been proposed in the literature for all the techniques we use, and for several of them corresponding packages are freely available (We provide additional details on estimation, computational cost and available packages in Section 3 of the Supplementary Material). Primarily, the choice of the design/package is subordinated to the goal of the analysis; however, computational time and resources may impose constraint on the analyst’s selection in realistic application. Clearly, there is a trade-off among number of insights, accuracy in estimation, and the time and resources required for the analysis. Then, the strategy in this case is to exploit the same model runs to gain insights on alternative settings simultaneously. Because we are considering a sensitivity analysis within an uncertainty quantification, the nominal minimal cost is , the size of an available Monte Carlo sample. Post-processing this sample to estimate variance-based, moment independent sensitivity measures, as well as partial dependence functions, the analyst obtains indications for the three settings of factor prioritization, trend analysis and interaction quantification. If computational resources afford a strategy based on the randomized calculation of first order finite change indices, then the analyst can compute all sensitivity measures in Table 1, with the exception of second order Newton ratios. To quantify these indices as well, additional model runs are needed; it is then the analyst’s choice to decide whether/how to carry out additional simulation experiments, in consideration of computational burden.

Application: sensitivity for COVID-19 modeling

In this section, we illustrate the insights that can be obtained by applying a systematic sensitivity analysis in the context of a generalized SIR model developed to forecast the expansion of the early phases of the COVID-19 pandemic. In Section 4.1 we present the generalized SIR model; in Section 4.2 we present the data and the model calibration; in Section 4.3 we discuss results for uncertainty quantification; and in Section 4.4 we present results for sensitivity analysis concerning key drivers of uncertainty, interaction quantification and trend analysis.

The generalized SIR model

We start from the version used to simulate the COVID-19 outbreak in China in the work of Peng et al. (2020). Note that a modification of the model is used in Zhang et al. (2020) and López & Rodó (2020). The model belongs to the family of susceptible-exposed-infectious-recovered (SEIR) models.

In total, the model considers seven states: susceptible , insusceptible , exposed (but not yet infectious) , infectious , quarantined , recovered and deceased (Fig. 1 ). The insusceptible compartment is introduced in Peng et al. (2020) to model the effect of potentially introduced public health measures, such as social distancing policies and contact tracing that can counterbalance the initial spread of infection. For simplicity, natural birth and death processes are not considered so that a constant population () is assumed, with .

Fig. 1

The generalised SIR model; adapted from Peng et al. (2020). The bold part in the equations corresponds to the classical SEIR model.

The left panel in Fig. 1 reports the flowchart with the links among compartments, and the right panel displays the corresponding differential equations. The bolded equations are present in the classical SEIR model. As per the first two equations in the right panel, the shift from susceptible to insusceptible is associated with the protection rate , and from susceptible to exposed with the infection rate . The third equation indicates that exposed individuals may enter the infective state at rate . The fourth equation indicates that infected individuals become quarantined at rate . The fifth equation models the fact that individuals who leave the infective state are either recovered or deceased. The cure () and mortality rate () in the last three equations are time dependent. Specifically, the cure rate is written as , and the mortality rate is written as . The model assumes that recovered individuals cannot become re-infected.2 The Matlab implementation of the model is freely accessible at https://github.com/ECheynet/SEIR.

Data collection and model calibration

We calibrated the model on openly available data for the COVID-19 pandemic from Italy, Germany, Spain, and the U.S. For the sake of space, we report results for Italy and the U.S. — results for Spain and Germany are in the supplementary material. For Italy, we collected data for the period (February 24th to April 20th, 2020), from Italy’s Department of Civil Protection publicly available repository https://github.com/pcm-dpc/COVID-19. For the U.S., we collected data from the COVID-19 Data Repository by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU CSSE COVID-19 Data) (Dong, Du, & Gardner, 2020) at https://github.com/CSSEGISandData/COVID-19, for the period March 4th to April 29th, 2020. The dates are chosen such that both periods comprise 56 days from the initial date. We first discuss the calibration when the model is fitted to data concerning the evolution of the pandemic in Italy (for Italy, henceforth, for brevity).

We divide the fitting into two periods, pre- and post-intervention. For the pre-intervention period, we calibrate the model with data from February 24th to March 8th, and use the remaining data for the post-intervention period. We fit the parameters in periods by maximizing the model coefficient of determination () between the actual and forecasted numbers of quarantined, recovered and deceased individuals. Results are reported in Fig. 2 and Table 2. Fig. 2 displays the forecasted trajectories in Italy (continuous lines) against the real observations (circles), including the number of quarantined, recovered, deceased and total infected individuals (note that the total number of infected individuals is the sum of the number of quarantined, recovered and deceased individuals). Table 2 shows that the average coefficient of model determination equals for pre-intervention data, and in the post-intervention phase.

Fig. 2

Italy: Model fitting pre/post-intervention; points are real data, lines are the fitted model predictions.

Table 2

Estimated model inputs for Italy and the United States.

	Days	R^2	α	β	γ	δ	λ_0,λ_1	κ_0,κ_1
Italy
Pre-intervention	24-Feb to 08-Mar (13 days)	0.9542	0.0000	1.1801	0.4633	0.5985	0.0437,0.1161	0.0162,0.0461
Post-intervention	09-Mar to 20-Apr (42 days)	0.9980	0.1098	2.0000	0.0704	0.3750	0.0167,2.0000	0.0240,0.0432
The United States
Pre-intervention	04-Mar to 16-Mar (12 days)	0.9441	0.0309	1.3466	0.5414	0.4835	0.0133, 0.0056	0.0098, 0.0736
Post-intervention	17-Mar to 29-Apr (43 days)	0.9968	0.0100	1.1905	0.9993	0.9149	0.0076, 0.1714	0.0086, 0.0279

Table 2 also reports the resulting parameters for Italy: the post-intervention parameter estimates are significantly different from the pre-intervention ones. For instance, the protection rate () increases by over 10,000%; the average latent time () by about .

Proceeding in a similar fashion for the U.S. data, we obtain the results in rows 5 to 8 of Table 2. The pre- and post-intervention models for the U.S. register average greater than 0.94. Also here we find a significant difference between the pre- and post-intervention parameter values.

Overall, our calibration exercise shows a reasonable fit of the model with the available data. However, the analysis carried out thus far is purely deterministic and does not take into account uncertainty around the fitted values. Indeed, Peng et al. (2020) assign ranges to these inputs and carry out a one-at-a-time sensitivity analysis, setting the inputs to the extreme of their ranges. As mentioned, this type of analysis provides only a limited exploration of the model input space.

Uncertainty quantification

In this section, we consider variability in model predictions that originates when the uncertain inputs are allowed to vary simultaneously over their entire supports. The uncertain inputs are the protection rate (), the infection rate (), the average latent time (), the average quarantine rate (), the number of initially infected individuals () and the intervention time (Interv.).3

To assign corresponding distributions, we proceed as follows. For the post-intervention input parameters , , and , we consider the same ranges assigned in previous literature (see Peng et al., 2020). For the number of initially infected individuals () and Intervention time (Interv.), we assign the discrete ranges reported in Table 3 . All input parameters are assumed to be uniformly distributed. To test the stability of these findings, we perform experiments in the presence and absence of correlations. We hypothesize a positive correlation between , the protection rate and , the quarantine rate. The rationale is that the higher the protection rate, the higher the efficacy of quarantine measures should be. However, this needs to be regarded as a thought experiment. In fact, an exploration of related literature has not allowed us to find works in which SEIR input parameters have actually been correlated. (Some works consider correlations at a higher level: For example, Keeling, Rand, & Morris (1997) and Dottori & Fabricius (2015) consider the probabilistic dependence between susceptible individuals and their infective neighbours, that is, a correlation between and , which in our case, are model outputs.)

Table 3

Uncertainty quantification: model input distributions. The post-intervention parameters , , and are assigned distributions over the support assigned in Peng et al. (2020). The distributions are centred at the corresponding values in the fourth and last rows of Table 2, respectively. The intervention time is assumed to be uniformly distributed over a one-week range after the nominal issuance date. For experiments with correlations, a 0.5 correlation between the protection () and the quarantine rate () is hypothesized.

	α	β	γ^−1	δ	I_0	Intervention time (Interv.)
Distribution	Uniform	Uniform	Uniform	Uniform	Discrete uniform	Discrete uniform
Support	α±10%	β±10%	γ^−1±10%	δ±30%	I_0±20%	{Intervention Day+z}, with z=0,1,…,7

We then proceed with generating two Monte Carlo samples of size , one for the independent and one for the dependent case, and run the SEIR model. We consider as output the total number of infected individuals at the end of the time series (, henceforth). Time wise, the entire analysis (simulator runs plus calculation of the sensitivity measures) takes about two hours per sample on a personal computer with 64-bit Windows 10 operating system, Intel® Xeon® processor (3.60 GHz) and 64 GB RAM.

Panels (a) and (b) in Fig. 3 present results for Italy, with independent and dependent inputs, respectively. The dotted line is the actual trajectory of the number of infected individuals as per the recorded data. The shaded areas correspond to the prediction interval (from the 2.5th to the 97.5th quantiles) for the entire trajectory from February 24th to April 20th. Note that the actual recorded data fall in the prediction intervals at all dates. Panels (c) and (d) in Fig. 3 display snapshots of the empirical distributions of on the end date (April 20th for Italy), for the cases of independent and dependent inputs, respectively. The empirical densities are right-skewed in all cases (details on the means, medians and standard deviations are reported in Panels (c) and (d) of Fig. 3). With correlations (Panel (d)), the estimated mean, median and standard deviation increase by about or less than , while the skewness increases by .

Fig. 3

Uncertainty quantification: Empirical densities for the number of infected individuals in Italy on April 20th and the U.S. on April 29th; dark blue circles in Panels (a) and (b) are the actual number of infected individuals; vertical thick blue lines in (c)–(f) correspond to actual recorded values; vertical light blue lines in (c)–(f) correspond to the 2.5th and 97.5th percentiles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We perform a similar analysis for the U.S. The period is now March 4th to April 29th, 2020. Panels (e) and (f) report the snapshots of the empirical distributions of on the end date. The presence of correlation has a negligible influence in the median and skewness, with changes smaller than . However, the mean and standard deviation increase by and , respectively. Overall, for both countries, model forecasts are consistent with the recorded data, with the actual number of infected individuals (dark blue vertical line in Fig. 3) falling in the prediction interval [2.5th quantile, 97.5th quantile] (light vertical blue lines in Fig. 3).

Uncertainty quantification is essential to ensure transparency in communications between scientists, stakeholders and the public (Bayarri, Paulo, Berger, Sacks, Cafeo, Cavendish, et al., 2007, Dunson, 2018, IHME, Saltelli, Bammer, Bruno, Charters, Di Fiore, Didier, et al., 2020). Results from uncertainty quantification need to be critically examined before making decisions based on a simulator’s forecasts. To illustrate, suppose that a panel of experts includes one or more quantitative measures of dispersion to decide whether or not to accept a model forecast. For instance, the panel may choose to proceed with recommendations if the coefficient of variation (the ratio between standard deviation and mean value of ) is smaller than 50% and if, simultaneously, the 95% prediction interval is smaller than the mean. Then, consider predictions for US as in Fig. 3. The coefficient of variation at the end of the period equals approximately 200% and the 95% prediction interval is about eight times higher than the mean. Therefore, the forecasts would not meet the uncertainty quantification acceptability criteria. The experts would then need to take further actions. A first step is understanding where to allocate further modeling efforts and data collection. Methods for global sensitivity analysis provide guidance in this respect, as we discuss in the next subsection.

Global sensitivity analysis results

Factor prioritizationFig. 4 displays the values of the three selected global sensitivity measures (total order indices , first order indices and — see Section 3 and Appendix A for definitions) for Italy and the U.S. Panels (a) and (b) report results for Italy, under independence and dependence, respectively. Panel (a) in Fig. 4 shows that all sensitivity measures concur in identifying intervention time (Interv.) as the key driver of prediction variability for Italy. This input is followed by quarantine rate (), initially infected individuals () and infection rate (); protection rate () and average latent time () play a minor role. On a relative scale, intervention time is to 6.18 times more important than quarantine rate, depending on the sensitivity measure.

Fig. 4

Importance of the six inputs according to the global sensitivity measures. The bars over report the total order global sensitivity indices of the uncertain inputs.

With correlations, Panel (b) shows a change in ranking, with protection rate () becoming the third most influential input parameter, both according to (or ) and . At the same time, the estimates of vary only slightly. This confirms that the increase in relevance of is attributable to its statistical correlation with , and not to an increase in its structural relevance. For the U.S., Panel (c) shows that the most important input is the quarantine rate (), followed by the infection rate () and Interv. Specifically, is about 2.38 to 5.86 times more important than , and 6.79 to 275.30 times more important than Interv. This is a notable difference with respect to the results in Panel (a), and it seems to be in line with the fact that the U.S. adopted a state-based approach with a gradual imposition of intervention measures in the first phase of the pandemic.4 With correlations (Panel (d)), the importance of the protection rate () increases, with this input ranking second and becoming third, according to both and . The estimates of and the corresponding ranking, instead, remain basically identical. Note the similarity of these results to the ones of Panel (b).

We now consider a practical use of these insights. Let us come back to the use of global sensitivity analysis within an uncertainty quantification, and assume that, given the information on the key-uncertainty drivers, the analyst gathers additional information on the two most important inputs for the independent case in the U.S., namely and and consider that the data collection yields a change in the input distributions from uniform to a on the same support. With this new assignment, the mean values of the input parameters remain the same, but their variances are reduced by about . Repeating the uncertainty quantification with these new distributions and comparing to the results in Section 4.3 (see supplementary Section 2 for additional details), the coefficient of variation now decreases from 2.08 (Uniform) to 1.78 (Beta), and the prediction interval drops from 7.82 (Uniform) to 5.90 (Beta) times the mean. Then, the predictions would still not be included in the panel decision-making process.

Trend analysis For this setting, we employ both partial dependence and one-way sensitivity functions. We use the R package ICEbox (Goldstein, Kapelner, Bleich, & Kapelner, 2017) — see Section 3 in Supplementary Material for additional details — on a subset of size , consisting of randomly selected points from our original samples, with correlated and uncorrelated inputs.

Fig. 5 displays the individual conditional expectation (ICE) plots of the inputs for Italy (Panels (a) and (b)) and the U.S. (Panels (c) and (d)). Each panel contains a scatterplot of the model output against the input of interest (dots); the highlighted thick line is the graph of a partial dependence function ; each grey line is a one-way sensitivity function () that visualizes how the forecast of changes as a function of given that the other inputs are fixed at a particular value.

Fig. 5

ICE and partial dependence (highlighted thick line) plots for Italy (Panels (a) and (b)) and for the U.S. (Panels (c) and (d)). The results are obtained using the ICEbox package (Goldstein et al., 2017) and a subsample of 10,000 Monte Carlo runs.

To illustrate, the first graph in Panel (a) shows how varies as a function of the protection rate () for Italy with independent inputs: decreases slightly as increases. The last graph shows a strong and increasing dependence of on the intervention time (Interv.). For instance, a delay of seven days would cause an increase in the mean of of around (This increase is calculated by comparing the mean conditional on the intervention being effective on March 9th, 2020, 14 on the horizontal axis, to the mean conditional the intervention being effective on March 16th, 2020, 21 on the horizontal axis). We also register a non-negligible impact of the quarantine rate (fourth graph in the Italy panel), with decreasing as increases. Note that the impact is milder than that of Interv., consistently with the factor prioritization results. Panel (b) reports results for trend analysis in the correlated inputs case. Findings are similar to those in Panel (a); the trend is the same in all graphs as it is expected by the monotonicity consistency of the selected trend indicators. Panels (c) and (d) report results for the U.S. They show that the quarantine rate has a stronger effect on than for Italy, an effect that remains decreasing. Intervention time still has an increasing effect, although weaker than for Italy. Also for the U.S., insights remain unaffected by the presence of correlations.

Interaction quantification In the independent input case, the quantity provides the portion of the output variance explained by interactions. The two pie charts in Fig. 6 display this quantity in black against the individual contributions in different colors. The graphs show that interactions contribute to about of the variance of for the Italy data fitted model. For the U.S., interactions account for of the variance of . To corroborate these findings, we also compute the mean dimension of the model ( in Table 1; see Appendix A for quantitative details). We register a mean dimension of 1.17 for Italy and 1.34 for U.S., confirming the stronger effect of interactions when the model is calibrated with the U.S. data. The values of these mean dimensions in any case may signal that we would not expect interactions of a higher order than two to play an important role. To gain a more precise understanding of this point, however, the analyst needs to perform further analysis aimed at investigating higher order interactions. The fast execution time of the SEIR model at hand allows us to study interaction of all orders at randomized locations.

Fig. 6

The portions of the output variance due to interactions (black) versus individual contributions (in different colors), according to the variance-based global sensitivity measure.

We use a sample size of 20,000, replicating a full factorial design and decomposing the output changes via a finite difference operation (see the code discussed in Section 3 of the Supplementary Material). This design allows us to compute a total of 20,000 second, third, fourth, fifth and sixth order interaction indices (). The availability of these indices yields thorough information about the interactions in this model.

Panels (a) and (b) in Fig. 7 report results for Italy and Panels (c) and (d) for the U.S. The first graph in each panel displays the top 25 interaction effects among the inputs, ranked by the average magnitude of the corresponding sensitivity indices over the 20,000 realizations. The horizontal axis displays the interaction orders. The vertical bars evidence the relevance of pairwise interactions, followed by third order interactions, with fourth order interactions playing a minor role. No interactions of order five or six are registered. This result agrees with the indications of the mean effective dimension that point at a stronger relevance of lower order interactions.

Fig. 7

Interaction effects in Italy ((a) and (b)), and the U.S. ((c) and (d)). Graphs in the first column present the top 25 averages of the magnitudes of interaction effects of all orders among the inputs. Panels in the last three columns of each row display histograms of Newton quotients for the top three interactions: orange for positive effects and green for negative effects. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

More specifically, results reveal that, for both Italy and the U.S., the strongest interactions systematically involve pairs of most important inputs. Specifically, in Italy, the highest interaction is between intervention time (Interv.) and quarantine rate (), followed by the interactions between intervention time and infection rate (), and intervention time and number of initially infected individuals (). In the U.S., the strongest interaction is registered between the infection rate () and the quarantine rate (), followed by the interaction between infection rate () and average latent time ().

The remaining three graphs in each panel of Fig. 7 indicate the sign of the three highest interactions identified in the first graph of the row. The graphs report histograms of second order Newton quotients calculated at random locations (see Appendix A and Eq. (22) for detailed definitions). Locally positive interactions are colored in orange, and negative interactions are in green. To illustrate, for Italy, the second and third graphs of Panel (a) show that the interactions between quarantine rate () and intervention time are locally negative (antagonistic), while the interactions between the number of initially infected individuals () and intervention time are locally positive (synergistic) in all scenarios.

The negative value of the interaction between intervention time and quarantine rate evidences an overlap between the effects of the two interventions, as quarantine measures are a subset of lockdown measures. Finally, we repeat the interaction analysis for the correlated input case as well. Results are very similar, and we do not report them for the sake of space.

Final remarks

OR/MS investigations heavily rely on quantitative models to support emergency and disaster management in the context of the COVID-19 pandemic. Our work moves in the direction advocated by Murdoch, Singh, Kumbier, Abbasi-Asl, & Yu (2019) and Saltelli et al. (2020) to provide analysts with state-of-the-art methods that allow them to thoroughly analyze the behavior of the model at hand, as well as to obtain a rigorous uncertainty quantification. We have proposed an approach that exploits synergies among different methods to yield insights simultaneously on alternative sensitivity settings. We have seen that applying variance-based, moment independent and ICE partial dependence plots on the same dataset generated via a Monte Carlo simulation provide analysts with indications on factor prioritization, trend analysis and interaction quantification; a design based on the randomization of finite change sensitivity indices yields, in addition, total variance-based indices and insights about direction of change at randomized locations.

For the specific SEIR model and data used in our application, uncertainty quantification has displayed information about the variability of model predictions. An expert can then make an informed decision on whether these predictions are reliable enough to support actual decision-making. A global sensitivity analysis has revealed that such variability is driven by the intervention-related parameters in the model, rather than by the epidemiological parameters. The trend analysis suggests that a delay in the implementation of intervention measures can notably increase the expected total number of infected individuals. For interaction quantification, results show a strong interaction between intervention time and quarantine rate for Italy, and between quarantine and infection rates for the U.S. We have repeated the experiments in the case of input correlations, obtaining similar results.

From a modeling viewpoint, the strong impact of quarantine rate and intervention time suggest that these variables need to be subject to further scrutiny, either in terms of data collection or of modeling efforts with the goal of improving the model output predictions and reducing the overall uncertainty. Repeating the analysis after having hypothesized a reduction in their variance, leads to a notable reduction in output uncertainty. From a policy-making viewpoint, the fact that at least one input related to non-pharmaceutical interventions stands out as important for all examined countries signals the relevance of contact prevention measures. Consider a decision between light or strong non-pharmaceutical interventions: Suppose that in a particular region, results suggest quarantine rate as the most influential factor. The decision maker would be encouraged to increase the rigour of contact prevention measures, but without resorting to a full lockdown. The converse would happen if lockdown was the most important variable. Consider the case of Italy in Panel (a) of Fig. 5. The partial dependence plots show that every additional day of delay in lockdown causes an average increase of about 28,600 infected individuals at the end of the period. In order to obtain such a drop with an increase in the quarantine efficacy (), we would have to increase by over 40%, i.e., beyond its currently assigned range. Then, a lockdown is the only viable path.

While we have applied the approach to a generalized SIR model, applicability is not restricted to this type of model and can be extended to other types of simulators used to forecast or study the COVID-19 pandemic. However, modifications of the method are needed for certain families of models, such as agent-based simulators, for which some of the elements subject to a sensitivity analysis might not be numerical. The method is also directly applicable to models characterized by a fast execution time. This becomes a limitation when time-consuming simulators are of concern. In this case, analysts can still apply the method, but either with a reduced number of model runs, or, if the original model is well fitted by an emulator (Kleijnen, 2015), by employing the fast-running emulator in place of the original model. These aspects cannot be comprised within the scope of the present manuscript and are part of future research.