Shapley Effect Estimation using Polynomial Chaos

Adrian Stein Tarunraj Singh [email protected] Department of Mechanical and Aerospace Engineering, University at Buffalo (SUNY),
Buffalo, NY 14260-4400, USA

Abstract

This paper presents an approach for estimating Shapley effects for use as global sensitivity metrics to quantify the relative importance of uncertain model parameters. Polynomial Chaos expansion, a well established approach for developing surrogate models is proposed to be used to estimate Shapley effects. Polynomial Chaos permits the transformation of a stochastic process to a deterministic model which can then be used to efficiently evaluate statistical moments of the quantity of interest. These moments include conditional variances which are algebraically mapped to Shapley effects. The polynomial chaos based estimates of Shapley effects are validated using Monte Carlo simulations and tested on the benchmark Ishigami function and on the dynamic SEIR epidemic model and the Bergman Type 1 diabetes model. The results illustrate the correct ranking of uncertain variables for the Ishigami function in contrast to the Sobol indices and illustrates the time-varying rank ordering of the model parameters for the dynamic models.

keywords:

Uncertainty Quantification , Sensitivity Analysis , Sobol , Shapley , Polynomial Chaos Expansion

1 Introduction

Uncertainties are ubiquitous and are broadly classified as Aleatory (irreducible) and Epistemic (reducible). Uncertainty Quantification (UQ) subsumes uncertainty characterization which corresponds to a parametric or non-parametric representation of the uncertainty, and uncertainty propagation which refers to the process of synthesizing the uncertainty in the output due to uncertainties in the input of a model. UQ also includes Global Sensitivity analysis which endeavours to rank the uncertain inputs based on their contributions to the uncertainty in the output.

Quantifying the uncertainty of systems or processes has gained remarkable interest driven by interest in model reduction, decision making under uncertainty, reliability engineering and explainable machine learning [1]. Local sensitivity analysis is a decades-old approach which evaluates the variation in a system’s output to small perturbations about nominal parameters. Meanwhile, Global Sensitivity Analysis (GSA) endeavours to rank the importance of all input variables over their domain of operation. It should be noted that GSA requires prior knowledge about the distribution of the uncertain input variables. GSA began gaining popularity because of the profound growth of computational power during the end of the last century, where the two most common GSA approaches were (1) Morris method [2] and (2) Sobol Analysis [3]. The latter, Sobol analysis, is a variance-based analysis and requires the input variables to be independent. Sobol and Gershman [4] introduced a derivative-based approach which integrated over the uncertain domain served as the global sensitivity metric. This was modified by Kucherenko [5] and called the Derivative-based Global Sensitivity Measures (DGSM). At the same time Sobol’ and Kucherenko revealed a relationship between the Sobol indices and the DGSM[6], which is useful for high dimensional systems [7]. Nandi et al. [8] proposed using a sampling based approach for numerical integration called Conjugate Unscented Transform to more efficiently evaluate Kucherenko’s [5] DGSM metrics.

The global sensitivity indices are generally evaluated via Monte Carlo (MC) sampling or quasi-Monte Carlo sampling, and most recently, via Polynomial Chaos Expansion (PCE). The MC makes use of a “pick-and-freeze” method [3][9] while a PCE surrogate model can be developed using an intrusive or non-intrusive approach. In the intrusive method, one expands the governing equations using separation of variables and applying the Galerkin projection is called the intrusive method. The non-intrusive approach is a sample based approach and being a black-box approach makes coding rather easy [10] [11]. However, for systems with many input variables, it is well know that the PCE can be a computational burden [12]. An even deeper connection between GSA and PCE was made by Sudret [12], who showed that the Sobol indices can be algebraically mapped to the polynomial chaos (PC) -coefficients. Motivated by the computational burden of calculating DGSM measures via MC or quasi MC, Sudret and Mai [13] found a way to compute DGSM measures with PC-coefficients as well.

GSA has been applied in disparate domains since the recognition that it should be an integral part of mathematical modeling [1]. Harman and Johnston [14] described in detail how intrusive PCEs are applied on dynamical system, where they chose the Susceptible-Infected-Recovered (SIR) disease model, and compared it to MC results, while performing a Sobol analysis. A more complex disease model is the Susceptible-Exposed-Infected-Recovered (SEIR) disease model, which was used by Jensen et al. [15]. Based on Danish data for the spread of Covid-19, they performed a Sobol analysis via PCE. To compute Sobol indices faster, Blatman and Sudret [16] proposed a sparse PCE which was used by Lin et al. [17] in a 3D multiphysics model for lithium-ion batteries. In the control engineering field Singh et al. [18] used PCE for the design of robust input shapers. Sobol analysis has also been used in for the analysis of energy consumption of battery electric vehicles [19], Bayesian networks [20], type 1 diabetes modelling [21], state of charge estimation [22] or thermal runaway [23] of lithium-ion batteries [24], underwater gliders [25] or in NOx level sensitivity analysis in premixed burners [26].

Borgonovo [27] introduced a non-moment based global sensitivity measure which compares the distance between two probability density functions (pdf) as the basis of ranking the relative importance of uncertain input variables. Nandi and Singh [28] also introduced a variety of non-moment based global sensitivity measures based on statistical distances including the Wasserstein, Hellinger and Kolmogorov distances. They proposed using polynomial chaos surrogate models to alleviate the computational burden of estimating statistical distances. Both of these papers and one by Chun et al. [29] confirmed the incorrect ranking of the three uncertain variables of the Ishigami function generated by the Sobol indices. The variance based Sobol indices orders the variables in descending order of importance as “1,2,3”, while the non-moment based metrics ranks them as “2,1,3”. This discrepancy can be attributed to the fact that the Ishigami function is characterized by a bi-modal probability density function which the variance based Sobol metrics are unable to rank order correctly.

Shapley value is a concept proposed by Lloyd Shapley in 1953 [30] for a fair distribution of credit amongst players based on their value to the coalition, in a cooperative game setting. The Shapley effect quantifies the average marginal contribution of each player by considering all possible coalitions that the player could be a part of. One of the advantage of using Shapley effects as proxies for global sensitivity measures is that the uncertain input variables can be correlated [31]. The downside of using Shapley effects for a large number of uncertain variables, is the computational cost involved in evaluating the marginal cost of a player in a large selection of coalitions [9] [32] [33] [34] [35].

Owen [36] noted a unique characteristic of the Shapley effects, that is Shapley effects lie between the first-order and total Sobol indices. Iooss and Prieur [35] while investigating dependent input variables found that the “sandwich” effect doesn not always hold true. Owen articulated that Sobol indices are easy to calculate and could be used to bound the Shapley effects for independent variables.

Shapley effects have recently been used in machine learning for feature selection, model interpretability, and model reduction by pruning [37]. The domain of explainable or interpretable Artificial Intelligence is exploiting the use of Shapley effects to help comprehend the decisions of outputs of machine learning models since one of the shortcomings of machine learning based models is their opaqueness due to their black box architecture. A classic example of the use of Shapley effects is to allocate landing fees for a proposed runway; based on aircraft size which correlates to length of required runway [38]. Shapley effects have been used in hydraulic models [39], biomanufacturing process risk [40], power flow models for islanded microgrids [41], drug sensitivity [42], earth fissures [43], and fatigue blade load estimations [44]. Some researchers compared Sobol indices and Shapley effects such as: [45], [46], [47], [35]. Most researchers uses MC sampling to calculate the Shapley effects, such as: [48] [49] [50] [33] [34] [25], to name a few.

Recentely, new ideas have entered the GSA community. One of these new ideas, by Gayrard et al. [51], was to perform stochastic processing with independent increments as inputs. Another, by Song et al. [34], established the connection between the first-order, total, and Shapley effects using the concept of semivalues. Benoumechiara et al. [52] used a Gaussian Process (Kriging model) to calculate Shapley effects faster than with MC, where they presented their findings using violin plots. With a slightly different parametric choice for the Ishigami function than in [29] [27], Benoumechiara et al. [52] and Plischke et al. [32] ranked the order as “2,1,3” derived from Shapley effects. Benoumechiara et al. [52] acknowledged the relationship between Sobol and Shapley for linear Gaussian models and highlighted Sudret’s work [12] on computing Sobol indices from PC-coefficients. Their final remark states: It would be interesting to have such a decomposition for the Shapley effects. because even using Kriging to estimate Shapley effects is computationally expensive, especially in higher dimensions. Our paper is positioned likewise, by establishing a relationship between PC-coefficients and Shapley effects. We will present two static models; which are the Ishigami function in the scheme of [52] [32] and a user-defined function with four uncertain input variables. To illustrate the proposed approach on dynamical models, the SEIR disease model and the Bergman type 1 diabetes model are considered. We highlight that this paper is only focused on independent input variables with a hypercube representing the uncertain domain.

This paper is structured as follows: Section 2 describes the variance-based sensitivity analysis, which includes a description of Sobol analysis, Shapley value theorem, PCE, mapping between PC-coefficients, and Sobol indices, as well as the Borgonovo metric. Section 3 presents the static problems including the Ishigami function and a user-defined function. Section 4 illustrates the GSA on the SEIR disease model and on the Bergman type 1 diabetes model. Our concluding remarks are presented in section 5.

2 Review of Global Sensitivity Metrics

This section reviews well established global sensitivity metrics including Sobol, a variance based metric and the Borgonovo delta metric, a non-moment based metric. It also proposes the use of the Shapley effect as a metric and illustrates how PC-coefficients can algebraically map to Shapley effects.

2.1 Sobol indices

Assume a function is given as:

\displaystyle y=f(x),\>\>\>x\in\mathcal{R}^{n}

(1)

where $y$ is a scalar output and $x$ represents the input parameters, presumed to be real numbers. The ANOVA [53] decomposition can be written as [12] :

\displaystyle f(x_{1},...,x_{n})=f_{0}+\sum^{n}_{i=1}f_{i}(x_{i})+\sum_{1\leq i<j\leq n}f_{ij}(x_{i},x_{j})+...+f_{1,2,...,n}(x_{1},...,x_{n})

(2)

where $f_{0}$ is a constant, $f_{i}$ depends purely on $x_{i}$ and $f_{i,j}$ purely on $x_{i}$ and $x_{j}$ . Furthermore, the following condition holds:

\displaystyle\int^{1}_{0}f_{1,2,...,s}(x_{1},x_{2},...,x_{s})dx_{k}=0,\>\>1\leq k\leq s

(3)

which means that the input variables are independent. It can be stated that:


$\displaystyle f_{o}$	$\displaystyle=\int_{\mathcal{R}}f(x)dx=\operatorname{\mathbb{E}}\left[y\right]$	(4a)
$\displaystyle f_{i}(x_{i})$	$\displaystyle=\int_{\mathcal{R}}f(x)dx_{\sim{i}}-f_{0}=\operatorname{\mathbb{E}}\left[y\|x_{i}\right]-f_{0}$	(4b)
$\displaystyle f_{ij}(x_{i},x_{j})$	$\displaystyle=\int_{\mathcal{R}}f(x)dx_{\sim{i},{j}}=\operatorname{\mathbb{E}}\left[y\|x_{i},x_{j}\right]-f_{i}(x_{i})-f_{j}(x_{j})-f_{0}$	(4c)

where $\sim{i}$ denotes the integration over all variables except $i$ and implies that $x_{i}$ is a fixed input variable. It is clear that the total variance of $f(x)$ can be calculated as follows:

\displaystyle\mathbb{V}=Var\left(f(x)\right)=\int_{\mathcal{R}}f^{2}(x)dx-f_{0}^{2}

(5)

where it is possible to decompose the total variance as:

\displaystyle\mathbb{V}=\sum^{n}_{i=1}\mathbb{V}_{i}+\sum_{1\leq i<j\leq n}\mathbb{V}_{ij}+...+\mathbb{V}_{1,2,...,n}

(6)

where we define the set $s$ as $s=1,2,...,n$ . Similar to the Eqns. (4a)-(4c) we can represent the variances as:

	$\displaystyle\mathbb{V}_{i}$	$\displaystyle=\mathbb{V}\left(\operatorname{\mathbb{E}}\left[y\|x_{i}\right]\right)$		(7)
	$\displaystyle\mathbb{V}_{ij}$	$\displaystyle=\mathbb{V}\left(\operatorname{\mathbb{E}}\left[y\|x_{i},x_{j}\right]\right)-\mathbb{V}_{i}-\mathbb{V}_{j}.$		(8)

The Sobol indices can now be calculated as:

\displaystyle S_{i_{1},...,i_{s}}=\frac{\mathbb{V}_{i_{1},...i_{s}}}{\mathbb{V}}

(9)

where it is known that

\displaystyle\sum_{i}^{n}S_{i}+\sum_{1\leq i<j\leq n}S_{ij}+...+S_{1,2,...n}=1.

(10)

The total Sobol indices can be calculated as:

\displaystyle S_{T_{i}}=\sum_{\kappa_{i}}S_{1,...,s}

(11)

where $\kappa_{i}=\{(i_{1},...,i_{s}):\exists k,1\leq k\leq s,i_{k}=i\}$ .

2.2 Shapley value concept

The Shapley effect of each variable is defined as [35]:

	$\displaystyle Sh_{i}$	$\displaystyle=\frac{1}{d!}\sum_{u\subseteq N\setminus\{i\}}\|u\|!(d-1-\|u\|)!\left(c(u\cup\{i\})-c(u)\right)$		(12)
	$\displaystyle\leftrightarrow Sh_{i}$	$\displaystyle=\frac{1}{d}\sum_{u\subseteq N\setminus\{i\}}\binom{d-1}{\|u\|}^{-1}\left(c(u\cup\{i\})-c(u)\right)$		(13)

where $d$ is the number of variables, $c(...)$ is the worth or value of the coalition. $c(u\cup\{i\})$ is a coalition including variable $i$ and $c(u)$ excluding variable $i$ . Shapley effects are shown as importance measures representing the marginal contribution of variable $i$ to a coalition $u$ . The Shapley value concept has the following features [35]:

1.

Relevant when contribution of each player is unequal
2.

Cooperation among players is beneficial rather than working independently
3.

Shapley effects cannot be negative
4.

Shapley effects sum-up to $1$

These features lead to the assertions that:

1.

Coalition $u$ is at least the empty set $\varnothing$
2.

Examining Eqn. (13), it is clear that coalition $|u|$ cannot have more than $d-1$ players, because it is already a subset and it is excluding player $i$ from the set, so there can be at most $d-1$ coalitions of that type.

Iooss and Prieur [35] proposed a variance based metric for the worth of a coalition:

\displaystyle c(u)=\frac{Var\left(E[Y|X_{u}]\right)}{Var(Y)}.

(14)

By assigning uncertain variables as players in a cooperative game setting, Shapley effects could be used as a Global Sensitivity Metric.

2.3 Polynomial Chaos Expansion

A PCE has the form of [10] [18]:

\displaystyle y(\zeta)=\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta)

(15)

where $y$ is the output variable, $a_{i}$ are the PC-coefficients and $\Phi_{i}$ are the basis functions. $P$ is the degree of the PCE. The orthogonality condition of the basis functions is given by [10]:

\displaystyle\left<\Phi_{i}(\zeta),\Phi_{j}(\zeta)\right>=\int_{\Omega}\Phi_{i}(\zeta)\Phi_{j}(x)w(\zeta)d\zeta=g_{i}^{2}\delta_{ij}

(16)

where $g_{i}$ are the coefficients which remain from inner product and $\delta_{ij}$ is the Kronecker delta. $w(x)$ is the pdf of $\zeta$ . The Wiener-Askey scheme provides the appropriate polynomial basis functions for standard random variables as illustrated in Table LABEL:table:Wiener_Askey_scheme.

Table 1: Wiener-Askey scheme for polynomial chaos expansion.

k_{l}

and

k_{u}

are the user-defined lower and upper bounds respectively.

Distribution of $\zeta$	Polynomial family	Support
Gaussian	Hermite	$\left(-\infty,\infty\right)$
gamma	Laguerre	$\left[0,\infty\right)$
beta	Jacobi	$\left[k_{l},k_{u}\right]$
uniform	Legendre	$\left[k_{l},k_{u}\right]$

The PC-coefficients can be determined using an intrusive or non-intrusive approach. The intrusive approach requires evaluation of inner products and might not result in closed form expressions for some nonlinear functions. The non-intrusive approach can be considered a black-box approach which requires solving a least-squares problem to determine the coefficients. In this paper we use the intrusive approach where the coefficients are determined by exploiting Galerkin projections. Consider Eqn. (15) which can be approximated as:

\displaystyle y(\zeta)=f(\zeta)\approx f_{PC}=\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta).

(17)

The left-hand side of Eqn. (17) represents the equations which are known and need to be approximated with a $P$ order polynomial with coefficients $a_{i}$ as shown on the right-hand side of Eqn. (17). Assuming $\zeta$ is uniformly distributed and selecting for instance $P=2$ results in the PCE:

\displaystyle y(\zeta)=f(\zeta)\approx f_{PC}=a_{0}\Phi_{0}(\zeta)+a_{1}\Phi_{1}(\zeta)+a_{2}\Phi_{2}(\zeta).

(18)

The Galerkin projection of Eqn. (18) results in:


	$\displaystyle\int f(\zeta)\Phi_{0}(\zeta)w(\zeta)d\zeta=g_{0}^{2}a_{0}$		(19a)
	$\displaystyle\int f(\zeta)\Phi_{1}(\zeta)w(\zeta)d\zeta=g_{1}^{2}a_{1}$		(19b)
	$\displaystyle\int f(\zeta)\Phi_{2}(\zeta)w(\zeta)d\zeta=g_{2}^{2}a_{2}.$		(19c)

From Eqns. (19a)-(19c), the PC-coefficients $a_{0}$ to $a_{2}$ can be solved by dividing by the constant factors $g_{0}$ to $g_{2}$ . For problems with multiple random variable, the basis functions are determined by the tensor product of the one-dimensional basis functions as illustrated in Table LABEL:table:PCE_mixed_variables for two random variables both of which are uniformly distributed.

Table 2: Polynomials for two uncertain variables (both uniform distributed) in a system for a degree

P=2

	$1$	$\zeta_{1}$	$(3\zeta_{1}^{2}-1)/2$
$1$	$1$	$\zeta_{1}$	$(3\zeta_{1}^{2}-1)/2$
$\zeta_{2}$	$\zeta_{2}$	$\zeta_{1}\zeta_{2}$
$(3\zeta_{2}^{2}-1)/2$	$(3\zeta_{2}^{2}-1)/2$

The polynomial $\Phi_{ij}$ denotes the polynomial basis function, where $i$ and $j$ are the degrees of the $\zeta_{1}$ and $\zeta_{2}$ variables respectively. For this reason, we can write the PCE of the system as:

	$\displaystyle y(\zeta_{1},\zeta_{2})=f(\zeta_{1},\zeta_{2})\approx f_{PC}=a_{00}\Phi_{00}(\zeta_{1},\zeta_{2})+a_{10}\Phi_{10}(\zeta_{1},\zeta_{2})+a_{01}\Phi_{01}(\zeta_{1},\zeta_{2})...$
	$\displaystyle...+a_{11}\Phi_{11}(\zeta_{1},\zeta_{2})+a_{20}\Phi_{20}(\zeta_{1},\zeta_{2})+a_{02}\Phi_{02}(\zeta_{1},\zeta_{2}).$		(20)

Multiplying all combinations of $\Phi_{ij}$ with a maximum order of $P=2$ and integrating over $\zeta_{1}$ and $\zeta_{2}$ results in $6$ equations in 6 unknowns which permit determining the PC-coefficients $a_{ij}$ .

2.4 Mapping between PC-coefficients and Sobol Indices

The mean and variance of the random output expressed via the PCE can be calculated as follows [12]:

	$\displaystyle\bar{y_{PC}}$	$\displaystyle=\operatorname{\mathbb{E}}\left[f_{PC}(\zeta)\right]=\int_{\mathcal{R}}\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta)\underbrace{\Phi_{0}(\zeta)}_{=1}w(\zeta)d\zeta=a_{0}$		(21)
	$\displaystyle\mathbb{V}_{PC}$	$\displaystyle=Var\left(\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta)\right)=\int_{\mathcal{R}}\left(\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta)\sum^{P}_{i=0}a_{i}\Phi_{i}(\zeta)-a_{0}^{2}\right)w(\zeta)d\zeta=\sum^{P}_{i=1}a_{i}^{2}\operatorname{\mathbb{E}}\left[\Phi_{i}^{2}(\zeta)\right]$		(22)

Consider a univariate polynomial family and denote it as $C_{P}(x),P\in\mathbb{N}$ , where $P$ is the polynomial degree. $C_{P}$ could for instance represent the family of Legendre Polynomials. When more than one variable in a system is uncertain, this idea can be generalized to a multivariate polynomial of order $M$ (number of variables) and degree $P$ , where the multiplication of univariate polynomials creates the multivariate polynomial of degree $P_{m}$ . However, the resulting degree $P_{m}$ must be less or equal than $P$ . Therefore, we define:

\displaystyle\alpha=\{\alpha_{i}:i=1,...,M\},\>\>\>\>\alpha_{i}\geq 0,\>\>\>\>\sum^{M}_{i=1}\alpha_{i}\leq P

(23)

which introduces a set, where $\alpha_{i}$ is a natural number and corresponds to the degree of the i-th univariate polynomial with respect to the variable $x_{i}$ . Then, the multivariate polynomial $\Phi_{i}$ can be expressed as the multiplication of $M$ univariate polynomials:

\displaystyle\Phi_{i}=\Phi_{\alpha}(x_{1},...,x_{M})=\prod^{M}_{i=1}\Phi_{\alpha_{i}}(x_{i}).

(24)

The number of polynomials $n_{C}$ satisfying Eqn. (23) is given by [12]:

\displaystyle n_{C}=\frac{\left(M+P\right)!}{M!P!}.

(25)

Let us define $\psi_{i_{1},...,i_{s}}$ which consists of $n_{C}$ $\alpha$ -tuples where $i_{1},...,i_{s}$ corresponds to a set of variables which are chosen. Therefore, certain tuple elements are zero.

\displaystyle\psi_{i_{1},...,i_{s}}=\forall\alpha:\begin{cases}\alpha_{k}>0\>\>\>\text{,when}\>\>\>k\in\{i_{1},...,i_{s}\}\\ \alpha_{k}=0\>\>\>\text{,when}\>\>\>k\not\in\{i_{1},...,i_{s}\}\end{cases}

(26)

where $k=1,...,M$ . Eqn. (26) states that any $k$ is non-zero if it lies in $\{i_{1},...,i_{s}\}$ . Now, the ANOVA decomposition of $f_{PC}$ can be derived as:

	$\displaystyle f_{PC}=a_{0}+\sum_{1\leq i_{1}\leq M}\sum_{\alpha\in\psi_{i_{1}}}a_{\alpha}\Phi_{\alpha}(x_{i_{1}})+\sum_{1\leq i_{1}<i_{2}\leq M}\sum_{\alpha\in\psi_{i_{1},i_{2}}}a_{\alpha}\Phi_{\alpha}(x_{i_{1}},x_{i_{2}})+...$
	$\displaystyle...+\sum_{1\leq i_{1}<...<i_{s}\leq M}\sum_{\alpha\in\psi_{i_{1},...,i_{s}}}a_{\alpha}\Phi_{\alpha}(x_{i_{1}},...,x_{i_{s}})+...+\sum_{\alpha\in\psi_{1,2,...,M}}a_{\alpha}\Phi_{\alpha}(x_{1},...,x_{M}).$		(27)

Note that each term $\sum_{\alpha\in\psi_{i_{1},...,i_{s}}}a_{\alpha}\Phi_{\alpha}(x_{i_{1}},...,x_{i_{s}})$ reads like: The term $a_{\alpha}\Phi_{\alpha}(x_{i_{1}},...,x_{i_{s}})$ gets only considered if the tuple $\alpha$ purely depends on the variables $x_{i_{1}},...,x_{i_{s}}$ . Accordingly, the PC based Sobol indices can be calculated as [12]:

\displaystyle S_{i_{1},...,i_{s}}=\sum_{\alpha\in\psi_{i_{1},...,i_{s}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}.

(28)

Since the PC-coefficients can now be mapped to the Sobol indices, the question arises if the Shapley effects can also be calculated from the PC-coefficients. One could say that we just need to calculate the worth $c(\{...\})$ of all coalitions. For a system with $n$ unknowns, the Sobol indices can mapped to the worth $c(\{...\})$ of all coalitions as follows:


$\displaystyle c(\{i_{1}\})$	$\displaystyle=\sum_{\alpha\in\psi_{i_{1}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}$	(29a)
$\displaystyle c(\{i_{1},i_{2}\})$	$\displaystyle=\sum_{\alpha\subset\psi_{i_{1},i_{2}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}+\sum_{\alpha\in\psi_{i_{1},i_{2}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}$	(29b)
	$\displaystyle\>\>\vdots$
$\displaystyle c(\{i_{1},...,i_{s}\})$	$\displaystyle=\sum_{\alpha\subset\psi_{i_{1},...,i_{s}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}+\sum_{\alpha\in\psi_{i_{1},...,i_{s}}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}$	(29c)
	$\displaystyle\>\>\vdots$
$\displaystyle c(\{1,2,..,M\})$	$\displaystyle=\sum_{\alpha\subset\psi_{1,2,...,M}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}+\sum_{\alpha\in\psi_{1,2,...,M}}\frac{a^{2}_{\alpha}\operatorname{\mathbb{E}}\left[\Phi^{2}_{\alpha}\right]}{\mathbb{V}_{PC}}\approx 1.$	(29d)

Note that $\alpha\subset\psi_{i_{1},...,i_{s}}$ stands for all possible combinations of proper subsets of $\{i_{1},...,i_{s}\}$ whereas $\alpha\in\psi_{i_{1},...,i_{s}}$ stands for the intersecting set of $\{i_{1},...,i_{s}\}$ . For a set $\{i_{1},...,i_{s}\}$ we can derive $2^{s}-1$ proper subsets, because the empty set is a proper subset of every set except for the empty set itself. The worths presented in Eqns. (29a)-(29d) can be directly substituted into Eqn. (12).

If one maps PC-coefficients to Sobol indices, the Shapley effects can be calculated as:

\displaystyle Sh_{i}=S_{i}+\frac{1}{2}\sum_{1\leq i<j\leq M}S_{ij}+\frac{1}{3}\sum_{1\leq i<j<k\leq M}S_{ijk}+...+\frac{1}{M-1}\sum_{1\leq i<...<s\leq M}S_{i,...,s}+\frac{1}{M}S_{1,2,...,M}.

(30)

Figure 1 shows the difference between calculating Shapley effects with Eqn. (14) and PC-coefficients. This sheds light on the relation of Sobol indices and Shapley effects. Moreover given all values of $c(u)$ , it is possible to map to the corresponding $S_{i},S_{ij},\ldots$ and vice versa.

Refer to caption — Figure 1: Connection between Sobol indices and Shapley effects.

2.5 Borgonovo metric

The Borgonovo metric [27] is a method which measures the shift between the unconditional and conditional pdfs, which is illustrated in Figure 2.

Assume that a function $Y$ has $i$ uncertain variables. The shift can then be used to provide a ranking order of importance for all uncertain variables. The main idea is that variables with little impact on the output $Y$ show a small shift and variables with a large impact have a large shift. The definition of a shift, which is noted as $Z$ , between two pdfs can be written as [27]:

\displaystyle Z(X_{i})=\int\lvert f_{Y}(y)-f_{Y|X_{i}}(y)\rvert dy

(31)

where $f_{Y}(y)$ is the pdf of the output $Y$ and $f_{Y|X_{i}}(y)$ is the conditional pdf for a given variable $x_{i}$ . Furthermore, the expected shift between $f_{Y}(y)$ and $f_{Y|X_{i}}(y)$ can be written as:

\displaystyle E_{X_{i}}\left[Z(X_{i})\right]=\int f_{X_{i}}(x_{i})\left(\int\lvert f_{Y}(y)-f_{Y|X_{i}}(y)\rvert dy\right)dx_{i}.

(32)

Finally, the Borgonovo index [27] can be derived from the expected shift, which is:

\displaystyle\delta_{i}=\frac{1}{2}E_{X_{i}}\left[Z(X_{i})\right]

(33)

where $\delta_{i}$ is a moment independent sensitivity measurement for an uncertain variable $x_{i}$ with respect to an output $Y$ .

3 Algebraic Examples

This section will present two algebraic examples including the Ishigami function and a user-defined function.

3.1 Ishigami function

A benchmark problem for global sensitivity analysis is the Ishigami function, especially because of its bi-modal pdf. The function is given as:

\displaystyle Y

\displaystyle=\sin(x_{1})+a\sin(x_{2})^{2}+bx_{3}^{4}\sin(x_{1})

(34)

where $a=7$ and $b=0.1$ are parameters consistent with those used in the literature. $x_{1}$ , $x_{2}$ and $x_{3}$ are uncertain variables, where it is assumed that all variables are uniformly distributed between $-\pi$ and $\pi$ . The pdf of the Ishigami function is shown in Figure 3.

The expected value and the total variance can be analytically evaluated by integrating over all three uncertain variables:

	$\displaystyle E[Y]$	$\displaystyle=\mu=\left(\frac{1}{2\pi}\right)^{3}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\sin(x_{1})+7\sin(x_{2})^{2}+0.1x_{3}^{4}\sin(x_{1})dx_{1}dx_{2}dx_{3}$		(35)
	$\displaystyle Var(Y)$	$\displaystyle=\left(\frac{1}{2\pi}\right)^{3}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\left(\sin(x_{1})+7\sin(x_{2})^{2}+0.1x_{3}^{4}\sin(x_{1})-\mu\right)^{2}dx_{1}dx_{2}dx_{3}.$		(36)

To determine the Sobol indices, three first-order, three second-order and one third-order indices need to be computed. For the first-order Sobol index, the conditional variances can be evaluated as:

$\displaystyle E[Y\|x_{i}]$	$\displaystyle=\mu_{i}=\left(\frac{1}{2\pi}\right)^{2}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\sin(x_{1})+7\sin(x_{2})^{2}+0.1x_{3}^{4}\sin(x_{1})d\mathcal{X}_{\sim i}$	(37)
$\displaystyle Var(E[Y\|x_{i}])$	$\displaystyle=\frac{1}{2\pi}\int^{\pi}_{-\pi}\left(\mu_{i}-\frac{1}{2\pi}\int^{\pi}_{-\pi}\mu_{i}dx_{i}\right)^{2}dx_{i}$	(38)
$\displaystyle\rightarrow\mathbb{V}_{i}$	$\displaystyle=Var(E[Y\|x_{i}])$	(39)

where $d\mathcal{X}_{\sim i}$ means that the integral is calculated for all variables except $x_{i}$ . For the second-order Sobol indices, the conditional variances can be specifically calculated as:

$\displaystyle E[Y\|x_{ij}]$	$\displaystyle=\mu_{ij}=\frac{1}{2\pi}\int^{\pi}_{-\pi}\sin(x_{1})+7\sin(x_{2})^{2}+0.1x_{3}^{4}\sin(x_{1})d\mathcal{X}_{\sim i,j}$	(40)
$\displaystyle Var(E[Y\|x_{ij}])$	$\displaystyle=\left(\frac{1}{2\pi}\right)^{2}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\left(\mu_{ij}-\left(\frac{1}{2\pi}\right)^{2}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\mu_{ij}dx_{i}dx_{j}\right)^{2}dx_{i}dx_{j}$	(41)
$\displaystyle\rightarrow\mathbb{V}_{ij}$	$\displaystyle=Var(E[Y\|x_{ij}])-\mathbb{V}_{i}-\mathbb{V}_{j}.$	(42)

For the third-order Sobol the conditional variances can be specifically calculated as:

$\displaystyle E[Y\|x_{123}]$	$\displaystyle=\mu_{123}=\sin(x_{1})+7\sin(x_{2})^{2}+0.1x_{3}^{4}\sin(x_{1})$	(43)
$\displaystyle Var(E[Y\|x_{123}])$	$\displaystyle=\left(\frac{1}{2\pi}\right)^{2}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\left(\mu_{123}-\left(\frac{1}{2\pi}\right)^{2}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\int^{\pi}_{-\pi}\mu_{123}dx_{1}dx_{2}dx_{3}\right)^{2}dx_{1}dx_{2}dx_{3}$	(44)
$\displaystyle\rightarrow\mathbb{V}_{123}$	$\displaystyle=Var(E[Y\|x_{123}])-\mathbb{V}_{12}-\mathbb{V}_{13}-\mathbb{V}_{23}-\mathbb{V}_{1}-\mathbb{V}_{2}-\mathbb{V}_{3}.$	(45)

The Sobol indices from the analytical evaluation are presented in Table LABEL:table:Ishigami_Sobol_indices_analytical_and_PC to permit their comparison to those determined using the PCE surrogate model.

Table 3: Sobol indices of the Ishigami function.

Index	Analytical	PCE order
		$P=3$	$P=5$	$P=7$	$P=9$
$S_{1}$	0.3139	0.6582	0.3699	0.3166	0.3140
$S_{2}$	0.4424	0.0540	0.3545	0.4377	0.4423
$S_{3}$	0	0	0	0	0
$S_{12}$	0	0	0	0	0
$S_{13}$	0.2437	0.2878	0.2756	0.2456	0.2437
$S_{23}$	0	0	0	0	0
$S_{123}$	0	0	0	0	0
$S_{T,1}$	0.5576	0.9460	0.6455	0.5623	0.5577
$S_{T,2}$	0.4424	0.0540	0.3545	0.4377	0.4423
$S_{T,3}$	0.2437	0.2878	0.2756	0.2456	0.2437

Table LABEL:table:Ishigami_Sobol_indices_analytical_and_PC shows that, with increasing PC degree, the Sobol indices converge to the analytical solution. To determine the Shapley effects, the worth of a coalition can also be calculated analytically with Eqn. (14) and is presented in Table LABEL:table:Ishigami_Shapley_values_analytical_and_PC.

Table 4: Shapley effects of the Ishigami function.

Worth	Analytical	PCE order
		$P=3$	$P=5$	$P=7$	$P=9$
$c(\{1\})$	0.3139	0.6582	0.3699	0.3166	0.3140
$c(\{2\})$	0.4424	0.0540	0.3545	0.4377	0.4423
$c(\{3\})$	0	0	0	0	0
$c(\{1,2\})$	0.7563	0.7122	0.7244	0.7544	0.7563
$c(\{1,3\})$	0.5576	0.9460	0.6455	0.5623	0.5577
$c(\{2,3\})$	0.4424	0.0540	0.3545	0.4377	0.4423
$c(\{1,2,3\})$	1	1	1	1	1
$Sh_{1}$	0.4357	0.8021	0.5077	0.4395	0.4358
$Sh_{2}$	0.4424	0.0540	0.3545	0.4377	0.4423
$Sh_{3}$	0.1218	0.1439	0.1378	0.1228	0.1219

Note that the worth of an empty coalition is $c(\{\varnothing\})=0$ . Similar to the Sobol indices in Table LABEL:table:Ishigami_Sobol_indices_analytical_and_PC, Table LABEL:table:Ishigami_Shapley_values_analytical_and_PC illustrates the convergence of the coalitions worth to their analytical values when the PC degree $P$ is increased and therefore a convergence of the Shapley effects to the analytical solution can also be noted. The remaining questions is: What exactly makes computing Shapley effects different from the Sobol analysis? Sobol and Shapley can be used to rank the importance of uncertain variables in a system. This means that a higher ranked uncertain variable contributes more to the perturbation of the output variable. Taking the route of computing Shapley effects via the Sobol indices from Figure 1, we can try to look at the weighting of those indices. Table LABEL:table:Ishigami_Sobol_vs_Shapley_analytical_weights is specifically tailored for the Ishigami function and provides a deeper understanding of the weighting.

Table 5: Difference in weighting the Sobol indices when calculating total Sobol indices and Shapley effects for the Ishigami function.

Sobol Index	Weight for $S_{T}$	Weight for $Sh$	Analytical
$S_{1}$	1	1	0.3139
$S_{2}$	1	1	0.4424
$S_{3}$	1	1	0
$S_{12}$	1	1/2	0
$S_{13}$	1	1/2	0.2437
$S_{23}$	1	1/2	0
$S_{123}$	1	1/3	0

If we use the analytical values from Table LABEL:table:Ishigami_Sobol_vs_Shapley_analytical_weights we can calculate the Shapley effects from the Sobol indices as follows:


	$\displaystyle Sh_{1}=S_{1}+\frac{1}{2}S_{12}+\frac{1}{2}S_{13}+\frac{1}{3}S_{123}=0.4357$		(46a)
	$\displaystyle Sh_{2}=S_{2}+\frac{1}{2}S_{12}+\frac{1}{2}S_{23}+\frac{1}{3}S_{123}=0.4424$		(46b)
	$\displaystyle Sh_{3}=S_{3}+\frac{1}{2}S_{13}+\frac{1}{2}S_{23}+\frac{1}{3}S_{123}=0.1218.$		(46c)

Conversely, the total Sobol index can be calculated as follows:


	$\displaystyle S_{T,1}=S_{1}+S_{12}+S_{13}+S_{123}=0.5576$		(47a)
	$\displaystyle S_{T,2}=S_{2}+S_{12}+S_{23}+S_{123}=0.4424$		(47b)
	$\displaystyle S_{T,3}=S_{3}+S_{13}+S_{23}+S_{123}=0.2437.$		(47c)

The results from Eqns. (46a)-(46c) compared to Eqns. (47a)-(47c) show that Shapley is weighing the higher order Sobol indices lower compared to Sobol where each index is equally weighted. Evaluating the ranking from a Sobol perspective suggests to order the uncertain variables as $x_{1}$ , $x_{2}$ and $x_{3}$ , while Shapley suggests $x_{2}$ , $x_{1}$ and $x_{3}$ . In fact, we can note if only first-order Sobol indices exist and higher order ones are zero, the total Sobol indices will be identical with the Shapley effects.

As pointed out in [32] and [52], the uncertain variable $x_{2}$ is more important than $x_{1}$ . The clearest way of evaluating the ranking of the uncertain variables is to compute the Borgonovo index as described in subsection 2.5 or other statistical distances [28] because they are based on pdfs which can be computational expensive. $300,000$ MC evaluations to calculate the Borgonovo index reveals the uncertain variables are ranked as $x_{2}$ , $x_{1}$ and $x_{3}$ , which is the same rank order which Shapley suggests. Results of all three methods are listed in Table LABEL:table:Ishigami_Sobol_indices_and_Shapley_values_comparision.

Table 6: Total Sobol indices, Shapley effects and Borgonovo index for the Ishigami function.

Variable	Total Sobol	Shapley	Borgonovo
$x_{1}$	0.5576	0.4357	0.2392
$x_{2}$	0.4424	0.4424	0.4222
$x_{3}$	0.2437	0.1218	0.1957

3.2 User-defined function

We consider a polynomial function to illustrate the evaluation of the Shapley effects where:

\displaystyle Y=1.9x_{1}^{2}+2x_{2}^{2}+1.05x_{3}x_{1}^{3}+0.35x_{4}

(48)

where it is assumed that $x_{1}$ to $x_{4}$ are uniformly distributed and each variable lies between $-1$ and $1$ . Figure 4 illustrates the pdf of the user-defined function.

The Sobol indices are listed in Table LABEL:table:User_defined_Sobol_indices_analytical_and_PC where it is evident that for a PCE order of $P=4$ , we obtain the exact analytical solution. Since the highest order of the user-defined function is $4$ , a PCE of order $4$ can exactly reproduce the original function. Note that all higher order Sobol indices are $0$ except $S_{13}$ .

Table 7: Sobol indices of the user-defined function.

Index	Analytical	PCE order
		$P=1$	$P=2$	$P=3$	$P=4$
$S_{1}$	0.4169	0	0.4215	0.4215	0.4169
$S_{2}$	0.4619	0	0.4670	0.4670	0.4169
$S_{3}$	0	0	0	0	0
$S_{4}$	0.0530	1	0.0536	0.0536	0.0530
$S_{12}$	0	0	0	0	0
$S_{13}$	0.0682	0	0.0579	0	0.0682
$S_{14}$	0	0	0	0	0
$S_{23}$	0	0	0	0	0
$S_{24}$	0	0	0	0	0
$S_{34}$	0	0	0	0	0
$S_{123}$	0	0	0	0	0
$S_{124}$	0	0	0	0	0
$S_{134}$	0	0	0	0	0
$S_{234}$	0	0	0	0	0
$S_{1234}$	0	0	0	0	0
$S_{T,1}$	0.4851	0	0.4794	0.4794	0.4851
$S_{T,2}$	0.4619	0	0.4670	0.4670	0.4619
$S_{T,3}$	0.0682	0	0.0579	0.0579	0.0682
$S_{T,4}$	0.0530	1	0.0536	0.0536	0.0530

Table LABEL:table:User_defined_Shapley_values_analytical_and_PC illustrates the convergence of the coalition worths and Shapley effects when higher order PCE are used. Again, for a PCE order of $P=4$ , we can observe that the PCE approximation is identical to the analytical results.

Table 8: Shapley effects of the user-defined function.

Worth	Analytical	PCE order
		$P=1$	$P=2$	$P=3$	$P=4$
$c(\{1\})$	0.4169	0	0.4215	0.4215	0.4169
$c(\{2\})$	0.4619	0	0.4670	0.4670	0.4619
$c(\{3\})$	0	0	0	0	0
$c(\{4\})$	0.0530	1	0.0536	0.0536	0.0530
$c(\{1,2\})$	0.8788	0	0.8884	0.8884	0.8788
$c(\{1,3\})$	0.4851	0	0.4794	0.4794	0.4851
$c(\{1,4\})$	0.4699	1	0.4751	0.4751	0.4699
$c(\{2,3\})$	0.4619	0	0.4670	0.4670	0.4619
$c(\{2,4\})$	0.5149	1	0.5206	0.5206	0.5149
$c(\{3,4\})$	0.0530	1	0.0536	0.0536	0.0530
$c(\{1,2,3\})$	0.9470	0	0.9464	0.9464	0.9470
$c(\{1,2,4\})$	0.9318	1	0.9421	0.9421	0.9318
$c(\{1,3,4\})$	0.5381	1	0.5330	0.5330	0.5381
$c(\{2,3,4\})$	0.5149	1	0.5206	0.5206	0.5149
$c(\{1,2,3,4\})$	1	1	1	1	1
$Sh_{1}$	0.4510	0	0.4504	0.4504	0.4510
$Sh_{2}$	0.4619	0	0.4670	0.4670	0.4619
$Sh_{3}$	0.0341	0	0.0290	0.0290	0.0341
$Sh_{4}$	0.0530	1	0.0536	0.0536	0.0530

Evaluating the ranking from a Sobol perspective suggests to order the uncertain variables as $x_{1}$ , $x_{2}$ , $x_{3}$ and $x_{4}$ , while Shapley suggests $x_{2}$ , $x_{1}$ , $x_{4}$ and $x_{3}$ . Table LABEL:table:User_defined_Sobol_indices_and_Shapley_values_comparision provides the Borgonovo indices, where the ranking is identical to the ranking based on the Shapley effects.

Table 9: Total Sobol indices, Shapley effects and Borgonovo index for the user-defined function.

Variable	Total Sobol	Shapley	Borgonovo
$x_{1}$	0.4851	0.4510	0.2734
$x_{2}$	0.4619	0.4619	0.3309
$x_{3}$	0.0682	0.0341	0.0178
$x_{4}$	0.0530	0.0530	0.0800

4 Dynamic Systems

Global sensitivity analysis of dynamic models is necessary for control and forecasting problems. Therefore, for illustrative purposes, this section will present two models including a SEIR disease model [15] and the Bergman model, which is used to model Type 1 diabetes.

4.1 SEIR model

The nonlinear SEIR model is presented by the equations:


$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=-\beta\frac{I(t)S(t)}{N}$	(49a)
$\displaystyle\frac{\partial E}{\partial t}$	$\displaystyle=\beta\frac{I(t)S(t)}{N}-\sigma E(t)$	(49b)
$\displaystyle\frac{\partial I}{\partial t}$	$\displaystyle=\sigma E(t)-\gamma I(t)$	(49c)
$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=\gamma I(t)$	(49d)

where $\beta$ , $\sigma$ and $\gamma$ are assumed to be uncertain. $S$ , $E$ , $I$ and $R$ , stand for the susceptible, exposed, infected and recovered group respectively. Variable $N$ is the total population and is set to $1$ for simplicity. The uncertain parameters can be interpreted as follows:

1.

$\gamma=\frac{1}{\tau_{inf}}$ , where $\tau_{inf}$ is the time constant of the infectious state
2.

$\sigma=\frac{1}{\tau_{inc}}$ , where $\tau_{inc}$ is the time constant of the exposed state
3.

$\beta=R_{o}\gamma$ , where $R_{o}$ is the reproduction number

We assume that they are uniformly distributed as follows [14]:


	$\displaystyle\beta\sim U(2.5,5.5)$		(50a)
	$\displaystyle\sigma\sim U(0.5,1.5)$		(50b)
	$\displaystyle\gamma\sim U(0.5,1.5)$		(50c)

which can be further rewritten as:


	$\displaystyle\beta=\beta_{0}+\beta_{1}\zeta_{1}$		(51a)
	$\displaystyle\sigma=\sigma_{0}+\sigma_{1}\zeta_{2}$		(51b)
	$\displaystyle\gamma=\gamma_{0}+\gamma_{1}\zeta_{3}$		(51c)

where the subscript $(.)_{0}$ is the mean and $(.)_{1}$ is used to map the uniform distribution with arbitrary bounds to the interval $[-1,\>1]$ . PCE of the random SEIR variables are:


$\displaystyle S(t,\zeta_{1},\zeta_{2},\zeta_{3})$	$\displaystyle=\sum^{\infty}_{i=0}S_{i}(t)\Phi(\zeta_{1},\zeta_{2},\zeta_{3})=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}S_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})$	(52a)
$\displaystyle E(t,\zeta_{1},\zeta_{2},\zeta_{3})$	$\displaystyle=\sum^{\infty}_{i=0}E_{i}(t)\Phi(\zeta_{1},\zeta_{2},\zeta_{3})=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}E_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})$	(52b)
$\displaystyle I(t,\zeta_{1},\zeta_{2},\zeta_{3})$	$\displaystyle=\sum^{\infty}_{i=0}I_{i}(t)\Phi(\zeta_{1},\zeta_{2},\zeta_{3})=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}I_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})$	(52c)
$\displaystyle R(t,\zeta_{1},\zeta_{2},\zeta_{3})$	$\displaystyle=\sum^{\infty}_{i=0}R_{i}(t)\Phi(\zeta_{1},\zeta_{2},\zeta_{3})=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}R_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3}).$	(52d)

The polynomial $\Phi(\zeta_{1},\zeta_{2},\zeta_{3})$ consists of all tensor product of Legendre polynomial of $\zeta_{1}$ , $\zeta_{2}$ and $\zeta_{3}$ . The ordinary differential equations can now be expressed as:


	$\displaystyle\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\frac{dS_{ijk}(t)}{dt}\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})=-\left(\beta_{0}+\beta_{1}\zeta_{1}\right)...$
	$\displaystyle...\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\sum^{\infty}_{m=0}\sum^{\infty}_{n=0}\sum^{\infty}_{q=0}S_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})\times I_{mnq}(t)\Phi_{m}(\zeta_{1})\Phi_{n}(\zeta_{2})\Phi_{q}(\zeta_{3})$		(53a)
	$\displaystyle\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\frac{dE_{ijk}(t)}{dt}\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})=\left(\beta_{0}+\beta_{1}\zeta_{1}\right)...$
	$\displaystyle...\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\sum^{\infty}_{m=0}\sum^{\infty}_{n=0}\sum^{\infty}_{q=0}S_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})\times I_{mnq}(t)\Phi_{m}(\zeta_{1})\Phi_{n}(\zeta_{2})\Phi_{q}(\zeta_{3})...$
	$\displaystyle...-\left(\sigma_{0}+\sigma_{1}\zeta_{2}\right)\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}E_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})$		(53b)
	$\displaystyle\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\frac{dI_{ijk}(t)}{dt}\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})=\left(\sigma_{0}+\sigma_{1}\zeta_{2}\right)...$
	$\displaystyle...\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}E_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})-\left(\gamma_{0}+\gamma_{1}\zeta_{3}\right)\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}I_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})$		(53c)
	$\displaystyle\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\frac{dR_{ijk}(t)}{dt}\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})=\left(\gamma_{0}+\gamma_{1}\zeta_{3}\right)\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}I_{ijk}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3}).$		(53d)

To approximate the original system shown in Eqns. (49a)-(49d), the Galerkin projection onto the subspace of the orthogonal basis functions formed by the tensor product of the Weiner-Askey basis functions in each random dimension leads to the equations:


	$\displaystyle\frac{dS_{uvw}}{dt}=-\frac{\beta_{0}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}\sum^{P}_{m=0}\sum^{P-m}_{n=0}\sum^{P-m-n}_{q=0}S_{ijk}(t)I_{mnq}(t)\Psi_{2}(1)...$
	$\displaystyle...-\frac{\beta_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}\sum^{P}_{m=0}\sum^{P-m}_{n=0}\sum^{P-m-n}_{q=0}S_{ijk}(t)I_{mnq}(t)\Psi_{2}(\zeta_{1})$		(54a)
	$\displaystyle\frac{dE_{uvw}}{dt}=\frac{\beta_{0}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}\sum^{P}_{m=0}\sum^{P-m}_{n=0}\sum^{P-m-n}_{q=0}S_{ijk}(t)I_{mnq}(t)\Psi_{2}(1)...$
	$\displaystyle...+\frac{\beta_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}\sum^{P}_{m=0}\sum^{P-m}_{n=0}\sum^{P-m-n}_{q=0}S_{ijk}(t)I_{mnq}(t)\Psi_{2}(\zeta_{1})-\sigma_{0}E_{uvw}(t)-\frac{\sigma_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}E_{ijk}(t)\Psi_{1}\left(\zeta_{2}\right)$		(54b)
	$\displaystyle\frac{dI_{uvw}}{dt}=\sigma_{0}E_{uvw}(t)+\frac{\sigma_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}E_{ijk}(t)\Psi_{1}\left(\zeta_{2}\right)-\gamma_{0}I_{uvw}(t)-\frac{\gamma_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}I_{ijk}(t)\Psi_{1}\left(\zeta_{3}\right)$		(54c)
	$\displaystyle\frac{dR_{uvw}}{dt}=\gamma_{0}I_{uvw}(t)+\frac{\gamma_{1}}{\Psi_{3}}\sum^{P}_{i=0}\sum^{P-i}_{j=0}\sum^{P-i-j}_{k=0}I_{ijk}(t)\Psi_{1}\left(\zeta_{3}\right).$		(54d)

$\Psi_{1}$ , $\Psi_{2}$ and $\Psi_{3}$ are calculated as:

	$\displaystyle\Psi_{1}\left(F\left(\zeta_{1},\zeta_{2},\zeta_{3}\right)\right)=\int_{\Omega_{3}}\int_{\Omega_{2}}\int_{\Omega_{1}}F\left(\zeta_{1},\zeta_{2},\zeta_{3}\right)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})\Phi_{u}(\zeta_{1})\Phi_{v}(\zeta_{2})...$
	$\displaystyle...\times\Phi_{w}(\zeta_{3})p_{1}\left(\zeta_{1}\right)p_{2}\left(\zeta_{2}\right)p_{3}\left(\zeta_{3}\right)d\zeta_{1}d\zeta_{2}d\zeta_{3}$		(55)
	$\displaystyle\Psi_{2}\left(F\left(\zeta_{1},\zeta_{2},\zeta_{3}\right)\right)=\int_{\Omega_{3}}\int_{\Omega_{2}}\int_{\Omega_{1}}F\left(\zeta_{1},\zeta_{2},\zeta_{3}\right)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})\Phi_{m}(\zeta_{1})\Phi_{n}(\zeta_{2})\Phi_{q}(\zeta_{3})...$
	$\displaystyle...\times\Phi_{u}(\zeta_{1})\Phi_{v}(\zeta_{2})\Phi_{w}(\zeta_{3})p_{1}\left(\zeta_{1}\right)p_{2}\left(\zeta_{2}\right)p_{3}\left(\zeta_{3}\right)d\zeta_{1}d\zeta_{2}d\zeta_{3}$		(56)
	$\displaystyle\Psi_{3}=\int_{\Omega_{3}}\int_{\Omega_{2}}\int_{\Omega_{1}}\left(\Phi_{u}\left(\zeta_{1}\right)\right)^{2}\left(\Phi_{v}\left(\zeta_{2}\right)\right)^{2}\left(\Phi_{w}\left(\zeta_{3}\right)\right)^{2}p_{1}\left(\zeta_{1}\right)p_{2}\left(\zeta_{2}\right)p_{3}\left(\zeta_{3}\right)d\zeta_{1}d\zeta_{2}d\zeta_{3}$		(57)

and $p_{1}$ , $p_{2}$ and $p_{3}$ are the pdfs of $\zeta_{1}$ , $\zeta_{2}$ and $\zeta_{3}$ respectively. Note that $\Psi_{3}$ is simply a coefficients which is the result of the Galerkin projection of the left-hand side of Eqns. (53a)-(53d) onto the orthogonal basis functions.

4.1.1 Mean and variance

For demonstrative purposes, consider a group of infected people. The mean (expected value) can be calculated as follows:

$\displaystyle\operatorname{\mathbb{E}}\left[I(t,\zeta_{1},\zeta_{2},\zeta_{3})\right]$	$\displaystyle=\operatorname{\mathbb{E}}\left[\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}I_{ijk}(t)\Phi_{1}\left(\zeta_{1}\right)\left(\zeta_{2}\right)\left(\zeta_{3}\right)\right]$
	$\displaystyle=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}I_{ijk}(t)\int_{\Omega_{3}}\int_{\Omega_{2}}\int_{\Omega_{1}}\Phi_{i}\left(\zeta_{1}\right)\Phi_{j}\left(\zeta_{2}\right)\Phi_{k}\left(\zeta_{3}\right)p_{1}\left(\zeta_{1}\right)p_{2}\left(\zeta_{2}\right)p_{3}\left(\zeta_{3}\right)d\zeta_{1}d\zeta_{2}d\zeta_{3}$
	$\displaystyle=I_{000}(t)$	(58)

and the variance as:

	$\displaystyle\mathbb{V}\left[I(t,\zeta_{1},\zeta_{2},\zeta_{3})\right]=\operatorname{\mathbb{E}}\left[I(t,\zeta_{1},\zeta_{2},\zeta_{3})^{2}\right]-\left(\operatorname{\mathbb{E}}\left[I(t,\zeta_{1},\zeta_{2},\zeta_{3})\right]\right)^{2}$
	$\displaystyle=\operatorname{\mathbb{E}}\left[\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\sum^{\infty}_{m=0}\sum^{\infty}_{n=0}\sum^{\infty}_{q=0}I_{ijk}(t)I_{mnq}(t)\Phi_{i}(\zeta_{1})\Phi_{j}(\zeta_{2})\Phi_{k}(\zeta_{3})\Phi_{m}(\zeta_{1})\Phi_{n}(\zeta_{2})\Phi_{q}(\zeta_{3})\right]-I_{000}(t)^{2}$
	$\displaystyle=\sum^{\infty}_{i=0}\sum^{\infty}_{j=0}\sum^{\infty}_{k=0}\left(I_{ijk}\right)^{2}\int_{\Omega_{3}}\int_{\Omega_{2}}\int_{\Omega_{1}}\left(\Phi_{i}\left(\zeta_{1}\right)\right)^{2}\left(\Phi_{j}\left(\zeta_{2}\right)\right)^{2}\left(\Phi_{k}\left(\zeta_{3}\right)\right)^{2}...$
	$\displaystyle....p_{1}\left(\zeta_{1}\right)p_{2}\left(\zeta_{2}\right)p_{3}\left(\zeta_{3}\right)d\zeta_{1}d\zeta_{2}d\zeta_{3}-\left(I_{000}(t)\right)^{2}.$		(59)

Furthermore we can simply calculate the conditional variance as shown in Eqn. (7) and (8) which permit the calculation of the Sobol indices and Shapley effects.

We start by gauging the approximation accuracy in evaluating the GSA metrics using PCE as compared to MC estimates. To compute the sensitivities with a variance-based approach, we only need the expected values and the variances, which is why we compare the MC to the PCE in terms of expected value and total variance. For the MC simulation, $5,000$ samples were chosen and considered sufficient to capture the uncertainties of the system and represent the states behavior. The PCE was computed until a degree of $P=4$ . For the initial conditions it is assumed that $S_{000}=0.99$ , $E_{000}=0$ , $I_{000}=0.01$ and $R_{000}=0$ . All other polynomial chaos states are $0$ . Then we arbitrarily selected a simulation time of $15$ days. Figure 5 illustrates the time variation of the expected value. The red dashed line corresponds to results computed with the MC method for the susceptible, exposed, infected and recovered group. The blue lines illustrate the results from the PCE where the increase in order of the PCE expansion corresponds to the increased depth of the blue color of the solid lines. From the left side of Figure 5, it is evident that the PCE approximations are converging to the results of the MC simulation. On the right side of Figure 5 the absolute error over time between the MC and PCE is shown. The absolute error overall is fairly small; however, on average, a higher PCE degree yields a smaller error over time.

Figure 6 illustrates via a red dashed line, the time-variation of the total variance generated by the MC simulations of the susceptible, exposed, infected and recovered group. The blue solid lines represent the simulation results of the PCE, with the higher PCE degrees being denoted by the increasingly darker blue colors. Compared to the expected value, the total variance shows a higher dependence on the degree of the PCE. The left side of Figure 6 clearly illustrates that a low PC degree yields a poor approximation of the MC results; however, a degree of $P=4$ closely matches the MC results. The right side of Figure 6 illustrates the absolute error between the MC and the PCE for all degrees and states.

In disease forecasting, it is important to know when the infected group will peak, and at what magnitude. Given the PC-coefficients $I_{ijk}$ and basis functions $\Phi(\zeta_{1},\zeta_{2},\zeta_{3})$ , we can easily calculate the statistics of $I(t)$ , by sampling the variables $\zeta_{1}$ , $\zeta_{2}$ and $\zeta_{3}$ which are uniformly distributed over $[-1,1]$ . The PC-coefficients must be calculated once and, by sampling the random variables, different trajectories of $I(t)$ over time can be synthesized. From these simulations, the largest magnitude of $I(t)$ and the peak time can be calculated and plotted in a 2D histogram. Figure 7 illustrates a 2D heat-chart of the peak time and the magnitude of the infected group. It should be noted that the PCE is a lot faster in calculating such heat maps than a sample based approach such as classic MC methods. In Figure 7 the peak time is most likely around time unit $5$ , with a magnitude of approximately $0.18$ , which matches the observation from the expected value in Figure 5.

The remaining question is: Do Shapley effects differ from Sobol indices for dynamical system as they did for the static equations which are presented in subsection 3.1 and 3.2? Table LABEL:table:Ishigami_Sobol_vs_Shapley_analytical_weights indicates that Shapley is assigning less weight on higher order Sobol indices than when calculating the total Sobol indices. Figure 8 illustrates the importance via Sobol and Shapley of the three uncertain variables $\beta$ , $\sigma$ and $\gamma$ over time on the infected group. The solid lines represent the Shapley effects while the dashed lines are the total Sobol indices. In subfigure a) initially the variable $\gamma$ has a high importance for the infected group. With the given initial conditions and looking at the Eqns. (49a)-(49d), it can be noted the initial high importance of $\gamma$ is justified because the infected group reduces its value by $\gamma$ to the recovered group. At the initial time the variables $\beta$ and $\sigma$ have no impact on the infected group and are therefore ranked very low because the exposed group is initialized with $0$ . As time goes on the importance of $\gamma$ drops and the importance of $\beta$ and $\sigma$ rises, while $\beta$ is higher valued than $\sigma$ . This occurs because $\beta$ provides the growth of the exposed group which then goes over to the infected group. Looking back at Figure 5, the expected value peaks at time instant $6$ , and from Eqn. (49c), it is obvious that if $I$ becomes large, the variable $\gamma$ has a large impact. This is captured in Figure 8 a) as well. When the infected group is reduced, the importance of $\gamma$ shrinks and $\beta$ and $\sigma$ become again more important because they provide the growth of $I$ .

Figure 8 illustrates the time-variation of a) the total Sobol indices and Shapley effects are plotted on the same scale, where we note that Shapley effects sum up to $1$ for all times while this is not true for the total Sobol indices. However, Figure 8 a) also illustrates the relative importance of all the uncertain variables for both approaches. To answer the question of why the total Sobol indices and Shapley effects differ, we need to look at the higher order Sobol indices, which are shown in Figure 8 b) and c). Their different weighting results in divergence of the Shapley effects from the total Sobol indices. It can be seen that especially when these higher Sobol indices peak that the total Sobol indices don’t coincide with the Shapley effects (for instance around time $6$ ) and when they are nearly $0$ (around time $5$ ) that the Shapley effects and total Sobol indices are almost the same. However, it is rather difficult to conclude a “ranking” order from Figure 8 a).

Therefore, we introduce Figure 9 where the colors blue, red, and green represent the uncertain variables $\beta$ , $\sigma$ and $\gamma$ respectively. The ranking of each uncertain variables are provided between the initial and final time in an intuitive ranking scheme of $1$ st, $2$ nd and $3$ rd. A significant ranking difference of the uncertain variables can be seen between the times $1$ to $4$ , which is just before the expected value of the infected group peaks. This information might be very helpful for entities that are trying to take certain measures to slow down the growth of infections in a pandemic. The total Sobol indices rank variable $\sigma$ at $2$ nd place until almost time instant $3$ while Shapley suggests that the rank switch between $\sigma$ and $\gamma$ is happening earlier (around time $2.5$ ). After the expected value of the infected group peaks, the total Sobol indices and the Shapley effects are similar.

Owen [36] pointed out that the Shapley effect is “sandwiched” in between the first-order Sobol index and the total order Sobol index. This can be seen in Figure 10 for the uncertain variables $\beta$ , $\sigma$ and $\gamma$ , where the dashed, solid and dotted lines represent the total Sobol indices, Shapley effects, and the first-order Sobol indices respectively. The sandwiching size is particularly large for $\beta$ and $\sigma$ around time instant $6$ , when the higher order Sobol indices peak as indicated in Figure 8 b) and c). In contrast, the width of the sandwich is negligible around time instant $5$ , which is when the higher order Sobol indices are nearly $0$ .

4.2 Type 1 diabetes Bergman model

A popular model to study type 1 Diabetes (T1D) is the minimal Bergman model [21]:


$\displaystyle\dot{G}(t)$	$\displaystyle=-\left(X(t)+p_{1}\right)G(t)+p_{1}G_{b}+D(t)$	(60a)
$\displaystyle\dot{X}(t)$	$\displaystyle=-p_{2}X(t)+p_{3}\left(I(t)-I_{b}\right)$	(60b)
$\displaystyle\dot{I}(t)$	$\displaystyle=-p_{4}\left(I(t)-I_{b}\right)+U(t)$	(60c)

where $G$ is the blood glucose level, $X$ is the intermediate state and $I$ is the insulin concentration. The variable $D(t)$ is the meal disturbance term and is assumed to be the Fisher model [54]:

\displaystyle D(t)=\begin{cases}0&\>\>\>t<t_{m}\\ Be^{-d\left(t-t_{m}\right)}&\>\>\>t\geq t_{m}\end{cases}

(61)

where $B$ stands for the meal quantity and is assumed to be $B=28.98$ , calculated from a meal of $45$ gm Carbohydrates (CHO) and the meal time is given as $t_{m}=15$ min [21]. This simulation considers an insulin bolus at $t=0$ . For the simulation the insulin bolus is represented by the time-varying term $U(t)$ in Eqn. (60c). The input $U(t)$ is given by:

\displaystyle U(t)=\begin{cases}\frac{1000\cdot CHO(\text{in g})}{CR\cdot V_{i}}&\>\>\>0<t<1\\ 0&\>\>\>1<t\end{cases}

(62)

where $CR$ is the insulin-to-carb ratio and $V_{i}$ is the distribution volume-to-insulin. For a $45$ gm meal, $CR=18.477$ and $V_{i}=12$ , and therefore $U(0<t<1)=202.96$ mU/L. The uncertainties lie in the parameters $p_{1}$ , $p_{2}$ , $p_{3}$ and the initial glucose level $G_{0}$ . We assume that they are uniformly distributed as follows [21]:


$\displaystyle p_{1}$	$\displaystyle\sim U(0.0201,0.0373)$	(63a)
$\displaystyle p_{2}$	$\displaystyle\sim U(0.0198,0.0368)$	(63b)
$\displaystyle p_{3}$	$\displaystyle\sim U(3.525e^{-5},6.545e^{-5})$	(63c)
$\displaystyle G_{0}$	$\displaystyle\sim U(83.43,154.93).$	(63d)

When applying the PCE it should be noted that initial blood glucose is calculated by the Galerkin projections. Note that the initial blood glucose levels are divided by the coefficient which results from the Galerkin projection on the state equations. For the MC to PCE comparison, $10,000$ samples for the MC are chosen. The system is further initialized with $X_{0000}=0$ and $I_{0000}=15.3872$ , while the rest of the expanded states of $X$ and $I$ are initialized to $0$ . Figure 11 illustrates the expected value and total variance of the blood glucose level. A higher PCE degree leads to a smaller error between the MC and PCE approach. On the right side of Figure 11, it appears that even with a higher degree of the PCE, a small error remains.

Figure 12 a) shows the total Sobol indices and Shapley effects over time of $p_{1}$ , $p_{2}$ , $p_{3}$ , and $G_{0}$ with respect to the blood glucose level. As expected, the initial blood glucose level is the most important, meanwhile, as time progresses, the importance of $p_{3}$ peaks at approximately $25$ minutes, and then falls continuously. This is due to the insulin input in Eqn. (60c) which makes the insulin level rise, and therefore the uncertain parameters $p_{3}$ in Eqn. (60b) becomes more important. It can be observed that the importance of variable $p_{2}$ increases over time while the importance of variable $p_{1}$ is nominally the least important factor over all time. However, Figure 12 a) also shows that, in the end, the Shapley effects and total Sobol indices diverge due to the higher Sobol indices as shown in Figure 12 b), c), and d) which are relatively small except towards the end of the simulation time. The insets in Figure 12 a) illustrate the divergence; however, this does not change the rank order of the system.

Figure 13 illustrates the time-varying rank order of the uncertain variables with respect to the blood glucose level. The Shapley effects provides exactly the same answer as the total Sobol indices. This is due to the small values of the higher order Sobol indices which means they have a small impact on the calculation of the total Sobol indices and Shapley effects. Hence, both methods provide similar rankings. If we compare the higher order Sobol indices from the SEIR model as shown in Figure 8 with ones in Figure 12 we can see that the magnitude of the second-order Sobol indices in the SEIR model are one order higher than in the Bergman Type 1 Diabetes model. The third-order Sobol indices are even two orders higher.

Figure 14 presents the sandwiching effect of the First and Total Sobol indices on the uncertain variables $p_{1}$ , $p_{2}$ , $p_{3}$ and $G_{0}$ . As shown in Figure 12, the total Sobol indices do not differ too much from the Shapley effects except for $p_{1}$ during the end of the simulation time. Therefore, a wide sandwiching effect can be mostly observed over the latter part of the simulation for $p_{1}$ .

5 Conclusion

This paper proposes to use the Shapley value concept as a variance based global sensitivity metric and compares it to the well established variance based Sobol indices. The Shapley effects are evaluated for the Ishigami benchmark problem in addition to a polynomial model. It is noteworthy that the Shapley effects provide a ranking which is in conflict with the Total Sobol indices for the Ishigami function - a fact which has been illustrated using the Borgonovo delta function and other statistical distances. To address the issue of the computational burden in determining the Shapley effects, a polynomial chaos expansion has been proposed to be used to efficiently evaluate the total and conditional variances which are required to determine both the Sobol indices and Shapley effects. In effect, once a polynomial chaos model has been determined, the Sobol indices and the Shapley effects can be algebraically evaluated. The polynomial chaos approach has been used to illustrate the determination of Shapley effects of uncertain variables for an SEIR model and the Bergman Type 1 diabetes model.

The comparison between the total Sobol indices, Shapley effects and the Borgonovo metric illustrated the potential of ranking uncertain variables using Shapley effects. Apart from diabetes research, we aim to work on controller design techniques where the Shapley value concept could be used to desensitize controllers to uncertainties. In doing so, a mapping from polynomial chaos coefficients to Shapley effects will enable a low computational cost approach for robust controller design.

Acknowledgment

The authors acknowledge the support of this work by the US National Science Foundation through CMMI Award number 2021710.

References

[1] S. Razavi, A. Jakeman, A. Saltelli, C. Prieur, B. Iooss, E. Borgonovo, E. Plischke, S. Lo Piano, T. Iwanaga, W. Becker, S. Tarantola, J. H. Guillaume, J. Jakeman, H. Gupta, N. Melillo, G. Rabitti, V. Chabridon, Q. Duan, X. Sun, S. Smith, R. Sheikholeslami, N. Hosseini, M. Asadzadeh, A. Puy, S. Kucherenko, H. R. Maier, The future of sensitivity analysis: An essential discipline for systems modeling and policy support, Environmental Modelling & Software 137 (2021) 104954. doi:10.1016/j.envsoft.2020.104954.
[2] M. D. Morris, Factorial sampling plans for preliminary computational experiments, Technometrics 33 (2) (1991) 161. doi:10.2307/1269043.
[3] I. M. Sobol, Sensitivity analysis for non-linear mathematical models, Mathematical Modelling and Computational Experiment 1 (1993) 407–414.
URL https://cir.nii.ac.jp/crid/1570009751284338304
[4] I. M. Sobol’, A. Gresham, On an alternative global sensitivity estimators, in: 1995 Proceedings of SAMO, 1995, pp. 40–42.
[5] S. Kucherenko, M. Rodriguez-Fernandez, C. Pantelides, N. Shah, Monte carlo evaluation of derivative-based global sensitivity measures, Reliability Engineering & System Safety 94 (7) (2009) 1135–1148. doi:10.1016/j.ress.2008.05.006.
[6] I. M. Sobol’, S. Kucherenko, Derivative based global sensitivity measures and their link with global sensitivity indices, Mathematics and Computers in Simulation 79 (10) (2009) 3009–3017. doi:10.1016/j.matcom.2009.01.023.
[7] B. Iooss, P. Lemaître, A review on global sensitivity analysis methods (2014).
URL http://arxiv.org/pdf/1404.2405v1
[8] S. Nandi, T. Singh, P. Singla, Derivative based global sensitivity analysis using conjugate unscented transforms, in: 2019 American Control Conference (ACC), IEEE, 2019, pp. 2458–2463. doi:10.23919/ACC.2019.8815110.
[9] T. A. Mara, W. E. Becker, Polynomial chaos expansion for sensitivity analysis of model output with dependent inputs, Reliability Engineering & System Safety 214 (6) (2021) 107795. doi:10.1016/j.ress.2021.107795.
[10] K.-K. K. Kim, D. E. Shen, Z. K. Nagy, R. D. Braatz, Wiener’s polynomial chaos for the analysis and control of nonlinear dynamical systems with probabilistic uncertainties [historical perspectives], IEEE Control Systems 33 (5) (2013) 58–67. doi:10.1109/MCS.2013.2270410.
[11] D. Xiu, G. E. Karniadakis, The wiener–askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing 24 (2) (2002) 619–644. doi:10.1137/S1064827501387826.
[12] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (7) (2008) 964–979. doi:10.1016/j.ress.2007.04.002.
[13] B. Sudret, C. V. Mai, Computing derivative-based global sensitivity measures using polynomial chaos expansions, Reliability Engineering & System Safety 134 (3) (2015) 241–250. doi:10.1016/j.ress.2014.07.009.
[14] D. B. Harman, P. R. Johnston, Applying the stochastic galerkin method to epidemic models with uncertainty in the parameters, Mathematical biosciences 277 (2016) 25–37. doi:10.1016/j.mbs.2016.03.012.
[15] B. C. Jensen, A. P. Engsig-Karup, K. Knudsen, E. Augeraud, M. Banerjee, J.-S. Dhersin, A. d’Onofrio, T. Lipniacki, S. Petrovskii, C. Tran, A. Veber-Delattre, E. Vergu, V. Volpert, Efficient uncertainty quantification and variance-based sensitivity analysis in epidemic modelling using polynomial chaos, Mathematical Modelling of Natural Phenomena 17 (2022) 8. doi:10.1051/mmnp/2022014.
[16] G. Blatman, B. Sudret, Efficient computation of global sensitivity indices using sparse polynomial chaos expansions, Reliability Engineering & System Safety 95 (11) (2010) 1216–1229. doi:10.1016/j.ress.2010.06.015.
[17] N. Lin, X. Xie, R. Schenkendorf, U. Krewer, Efficient global sensitivity analysis of 3d multiphysics model for li-ion batteries, Journal of The Electrochemical Society 165 (7) (2018) A1169–A1183. doi:10.1149/2.1301805jes.
[18] T. Singh, P. Singla, U. Konda, Polynomial chaos based design of robust input shapers, Journal of Dynamic Systems, Measurement, and Control 132 (5) (2010) 76. doi:10.1115/1.4001793.
[19] M. Braband, M. Scherer, H. Voos, Global sensitivity analysis of economic model predictive longitudinal motion control of a battery electric vehicle, Electronics 11 (10) (2022) 1574. doi:10.3390/electronics11101574.
[20] R. Ballester-Ripoll, M. Leonelli, Global sensitivity analysis in probabilistic graphical models (2021).
URL http://arxiv.org/pdf/2110.03749v1
[21] S. Nandi, T. Singh, Global sensitivity analysis on the bergman minimal model, IFAC-PapersOnLine 53 (2) (2020) 16112–16118. doi:10.1016/j.ifacol.2020.12.431.
[22] G. Fan, X. Li, R. Zhang, Global sensitivity analysis on temperature-dependent parameters of a reduced-order electrochemical model and robust state-of-charge estimation at different temperatures, Energy 223 (5) (2021) 120024. doi:10.1016/j.energy.2021.120024.
[23] A. S. Yeardley, P. J. Bugryniec, R. A. Milton, S. F. Brown, A study of the thermal runaway of lithium-ion batteries: A gaussian process based global sensitivity analysis, Journal of Power Sources 456 (2020) 228001. doi:10.1016/j.jpowsour.2020.228001.
[24] X. Lai, Z. Meng, S. Wang, X. Han, L. Zhou, T. Sun, X. Li, X. Wang, Y. Ma, Y. Zheng, Global parametric sensitivity analysis of equivalent circuit model based on sobol’ method for lithium-ion batteries in electric vehicles, Journal of Cleaner Production 294 (3) (2021) 126246. doi:10.1016/j.jclepro.2021.126246.
[25] H. Wu, W. Niu, S. Wang, S. Yan, Sensitivity analysis of control parameters errors and current parameters to motion accuracy of underwater glider using sobol’ method, Applied Ocean Research 110 (1) (2021) 102625. doi:10.1016/j.apor.2021.102625.
[26] S. Yousefian, G. Bourque, R. F. D. Monaghan, Uncertainty quantification of nox emission due to operating conditions and chemical kinetic parameters in a premixed burner, Journal of Engineering for Gas Turbines and Power 140 (12) (2018) 2721. doi:10.1115/1.4040897.
[27] E. Borgonovo, A new uncertainty importance measure, Reliability Engineering & System Safety 92 (6) (2007) 771–784. doi:10.1016/j.ress.2006.04.015.
[28] S. Nandi, T. Singh, Global sensitivity analysis measures based on statistical distances, International Journal for Uncertainty Quantification 11 (6) (2021) 1–30. doi:10.1615/Int.J.UncertaintyQuantification.2021035424.
[29] M.-H. Chun, S.-J. Han, N.-I. Tak, An uncertainty importance measure using a distance metric for the change in a cumulative distribution function, Reliability Engineering & System Safety 70 (3) (2000) 313–321. doi:10.1016/S0951-8320(00)00068-5.
[30] L. S. Shapley, 17. a value for n-person games, in: H. W. Kuhn, A. W. Tucker (Eds.), Contributions to the Theory of Games (AM-28), Volume II, Princeton University Press, 1953, pp. 307–318. doi:10.1515/9781400881970-018.
[31] E. Algaba, V. Fragnelli, J. Sánchez-Soriano, Handbook of the shapley value, 1st Edition, Chapman & Hall/CRC series in operations research, Chapman & Hall/CRC, Boca Raton, 2019. doi:10.1201/9781351241410.
[32] E. Plischke, G. Rabitti, E. Borgonovo, Computing shapley effects for sensitivity analysis, SIAM/ASA Journal on Uncertainty Quantification 9 (4) (2021) 1411–1437. doi:10.1137/19M1304738.
[33] M. I. Radaideh, S. Surani, D. O’Grady, T. Kozlowski, Shapley effect application for variance-based sensitivity analysis of the few-group cross-sections, Annals of Nuclear Energy 129 (8) (2019) 264–279. doi:10.1016/j.anucene.2019.02.002.
[34] E. Song, B. L. Nelson, J. Staum, Shapley effects for global sensitivity analysis: Theory and computation, SIAM/ASA Journal on Uncertainty Quantification 4 (1) (2016) 1060–1083. doi:10.1137/15M1048070.
[35] B. Iooss, C. Prieur, Shapley effects for sensitivity analysis with correlated inputs: Comparisons with sobol’ indices, numerical estimation and applications (2019).
URL http://arxiv.org/pdf/1707.01334v7
[36] A. B. Owen, Sobol’ indices and shapley value, SIAM/ASA Journal on Uncertainty Quantification 2 (1) (2014) 245–251. doi:10.1137/130936233.
[37] B. Rozemberczki, L. Watson, P. Bayer, H.-T. Yang, O. Kiss, S. Nilsson, R. Sarkar, The shapley value in machine learning (2022).
URL http://arxiv.org/pdf/2202.05594v2
[38] S. C. Littlechild, G. F. Thompson, Aircraft landing fees: A game theory approach, The Bell Journal of Economics 8 (1) (1977) 186. doi:10.2307/3003493.
[39] L. Pheulpin, N. Bertrand, V. Bacchi, Uncertainty quantification and global sensitivity analysis with dependent inputs parameters: Application to a basic 2d-hydraulic model, LHB 108 (1) (2022) 221. doi:10.1080/27678490.2021.2015265.
[40] W. Xie, B. Wang, C. Li, D. Xie, J. Auclair, Interpretable biomanufacturing process risk and sensitivity analyses for quality–by–design and stability control, Naval Research Logistics (NRL) 69 (3) (2022) 461–483. doi:10.1002/nav.22019.
[41] H. Wang, Z. Yan, X. Xu, K. He, Evaluating influence of variable renewable energy generation on islanded microgrid power flow, IEEE Access 6 (2018) 71339–71349. doi:10.1109/ACCESS.2018.2881189.
[42] Y. Chen, L. Juan, X. Lv, L. Shi, Bioinformatics research on drug sensitivity prediction, Frontiers in pharmacology 12 (2021) 799712. doi:10.3389/fphar.2021.799712.
[43] Y. Li, N. Friedman, P. Teatini, A. Benczur, S. Ye, L. Zhu, C. Zoccarato, Sensitivity analysis of factors controlling earth fissures due to excessive groundwater pumping, Stochastic Environmental Research and Risk Assessment 36 (1) (2022) 201. doi:10.1007/s00477-022-02237-8.
[44] L. Schröder, N. K. Dimitrov, J. Aasted Sørensen, Uncertainty propagation and sensitivity analysis of an artificial neural network used as wind turbine load surrogate model, Journal of Physics: Conference Series 1618 (4) (2020) 042040. doi:10.1088/1742-6596/1618/4/042040.
[45] J. A. Carta, S. Díaz, A. Castañeda, A global sensitivity analysis method applied to wind farm power output estimation models, Applied Energy 280 (2020) 115968. doi:10.1016/j.apenergy.2020.115968.
[46] B. Broto, F. Bachoc, M. Depecker, J.-M. Martinez, Sensitivity indices for independent groups of variables, Mathematics and Computers in Simulation 163 (6) (2019) 19–31. doi:10.1016/j.matcom.2019.02.008.
[47] A. B. Owen, C. Prieur, On shapley value for measuring importance of dependent inputs (2017).
URL http://arxiv.org/pdf/1610.02080v3
[48] M. Il Idrissi, V. Chabridon, B. Iooss, Developments and applications of shapley effects to reliability-oriented sensitivity analysis with correlated inputs, Environmental Modelling & Software 143 (2021) 105115. doi:10.1016/j.envsoft.2021.105115.
[49] T. Goda, A simple algorithm for global sensitivity analysis with shapley effects, Reliability Engineering & System Safety 213 (2021) 107702. doi:10.1016/j.ress.2021.107702.
[50] M. B. Heredia, C. Prieur, N. Eckert, Global sensitivity analysis with aggregated shapley effects, application to avalanche hazard assessment, Reliability Engineering & System Safety 222 (4) (2022) 108420. doi:10.1016/j.ress.2022.108420.
[51] E. Gayrard, C. Chauvière, H. Djellout, P. Bonnet, D.-P. Zappa, Global sensitivity analysis for stochastic processes with independent increments, Probabilistic Engineering Mechanics 62 (1) (2020) 103098. doi:10.1016/j.probengmech.2020.103098.
[52] N. Benoumechiara, K. Elie-Dit-Cosaque, B. Bouchard, J.-F. Chassagneux, F. Delarue, E. Gobet, J. Lelong, Shapley effects for sensitivity analysis with dependent inputs: Bootstrap and kriging-based algorithms, ESAIM: Proceedings and Surveys 65 (2019) 266–293. doi:10.1051/proc/201965266.
[53] S. S. Kucherenko, et al., Global sensitivity indices for nonlinear mathematical models, review, Wilmott Mag 1 (2005) 56–61.
[54] M. E. Fisher, A semiclosed-loop algorithm for the control of blood glucose levels in diabetics, IEEE transactions on bio-medical engineering 38 (1) (1991) 57–61. doi:10.1109/10.68209.