Tutorials › Quantitative Systems Pharmacology (QSP) › Shapley Effects for Correlated QSP Parameters

Sensitivity Analysis and Model Reduction · Tutorial 667 of 799

Shapley Effects for Correlated QSP Parameters

Learn to calculate Shapley effects under a joint parameter distribution and interpret how dependence shares predictive influence across QSP parameters.

Advanced 12 min read

What You'll Learn

  • Define coalition values using conditional expectations under the joint parameter distribution.
  • Calculate Shapley effects by averaging marginal variance contributions over parameter orderings.
  • Interpret the allocation when parameters are correlated, without assuming independent inputs.
  • Work through exact two-parameter examples and a multi-parameter numerical calculation.
  • Identify estimation challenges, diagnostic checks, and limits on interpretation.

Why Shapley Effects Are Useful for Dependent Parameters

Global sensitivity analysis asks how uncertainty in model inputs relates to uncertainty in an output. In QSP, calibrated parameters may be dependent: for example, data may constrain a ratio of two rates more strongly than either rate separately. If a sensitivity measure assumes independent inputs, it can misrepresent how the model’s output uncertainty is distributed across those parameters. Shapley effects provide an allocation based on the joint input distribution, so dependence is part of the calculation rather than something removed by assumption.

The allocation answers a predictive question: how much of the output variance is assigned to each parameter when its marginal contribution is averaged over all possible orders in which parameters could be added to a set? It does not establish that a parameter causes the output to change, nor does it uniquely identify a biological mechanism.

Definition: Let \(\Theta=(\Theta_1,\ldots,\Theta_d)\) have a specified joint distribution, and let \(Y=f(\Theta)\) be a model output with finite, positive variance. For a subset \(S\) of the parameter indices, define its coalition value by $$ v(S)=\frac{\operatorname{Var}\!\left(\operatorname{E}[Y\mid \Theta_S]\right)}{\operatorname{Var}(Y)}. $$ Here \(\Theta_S\) denotes the parameters in \(S\). Set \(v(\varnothing)=0\). The Shapley effect for parameter \(i\) is $$ \phi_i=\sum_{S\subseteq N\setminus\{i\}} \frac{|S|!(d-|S|-1)!}{d!}\left[v(S\cup\{i\})-v(S)\right], $$ where \(N=\{1,\ldots,d\}\) is the full set of parameters. The weights average the contribution of adding \(i\) across all possible parameter orderings.

The conditional expectation \(\operatorname{E}[Y\mid\Theta_S]\) uses the distribution of the parameters not in \(S\), conditional on the values in \(S\). Thus, when parameters are dependent, learning \(\Theta_S\) also provides information about excluded parameters. The coalition value measures the fraction of output variance explained by that information, not merely by direct variation in the included parameters.

How the Allocation Works

For each coalition \(S\), the difference \(v(S\cup\{i\})-v(S)\) is the extra fraction of output variance explained by adding parameter \(i\) to \(S\). This increment can be large when \(i\) is added early and smaller when other correlated parameters are already known. Shapley effects average these increments rather than choosing one arbitrary order.

The main result is the efficiency property: the effects sum to one. To see why, consider any ordering of all \(d\) parameters. Each parameter contributes the increase in \(v\) from the preceding set to the set after it is added. Summing those increases telescopes from \(v(\varnothing)=0\) to \(v(N)=1\). Averaging this identity over all \(d!\) orderings gives the Shapley formula, because each parameter’s marginal contribution after a particular subset \(S\) occurs in exactly \(|S|!(d-|S|-1)!\) orderings. Therefore, $$ \sum_{i=1}^d \phi_i=v(N)-v(\varnothing)=1. $$ The equality \(v(N)=1\) follows because conditioning on all inputs determines \(Y=f(\Theta)\), so \(\operatorname{E}[Y\mid\Theta]=Y\).

This is a variance allocation, not a decomposition into independent effects. A parameter can receive an effect partly because it is informative about another parameter. If two parameters are nearly interchangeable as predictors under the joint distribution, averaging over orderings shares their predictive contribution rather than assigning it all to whichever parameter happens to be entered first.

Worked Example: Two Correlated Normal Parameters

Suppose standardized parameters \(\Theta_1\) and \(\Theta_2\) are jointly normal, each with mean 0 and variance 1, and correlation 0.6. Consider the hypothetical normalized output \(Y=\Theta_1+\Theta_2\). The finite-variance and joint-normal assumptions hold by construction, and the correlation is part of the specified input distribution.

The output variance is $$ \operatorname{Var}(Y)=1+1+2(0.6)=3.2. $$ For jointly normal variables, \(\operatorname{E}[\Theta_2\mid\Theta_1]=0.6\Theta_1\), so $$ \operatorname{E}[Y\mid\Theta_1]=1.6\Theta_1,\qquad v(\{1\})=\frac{1.6^2}{3.2}=0.8. $$ By symmetry, \(v(\{2\})=0.8\), and knowing both inputs gives \(v(\{1,2\})=1\). With two parameters, each parameter’s effect is the average of its contribution when entered first and when entered second: $$ \phi_1=\tfrac12[0.8-0]+\tfrac12[1-0.8]=0.5,\qquad \phi_2=\tfrac12[0.8-0]+\tfrac12[1-0.8]=0.5. $$

Each parameter is assigned 50% of the output variance. This equal allocation does not mean the parameters are independent or that each has an isolated causal effect. It reflects their symmetric roles in this model and the shared predictive information created by their correlation.

A Discrete Example of Shared Predictive Information

The conditional calculation also works for non-normal inputs. The next example makes dependence visible through a simple joint probability table.

Worked Example: A Correlated Parameter as a Proxy

Let two binary inputs have the hypothetical joint distribution \(P(\Theta_1=0,\Theta_2=0)=0.4\), \(P(1,0)=0.3\), \(P(0,1)=0.1\), and \(P(1,1)=0.2\). Let \(Y=\Theta_1\). These probabilities define the joint distribution, and \(Y\) has variance \(0.5(1-0.5)=0.25\).

Knowing \(\Theta_1\) determines the output, so \(v(\{1\})=1\). To calculate the value of \(\Theta_2\), first find the conditional means: $$ \operatorname{E}[Y\mid\Theta_2=1]=\frac{0.2}{0.3}=\frac23,\qquad \operatorname{E}[Y\mid\Theta_2=0]=\frac{0.3}{0.7}=\frac37. $$ Since \(P(\Theta_2=1)=0.3\), the variance of these conditional means is $$ 0.3\left(\frac23-\frac12\right)^2+ 0.7\left(\frac37-\frac12\right)^2 =0.01190. $$ Thus \(v(\{2\})=0.01190/0.25=0.04762\). Knowing both inputs gives \(v(\{1,2\})=1\).

The two Shapley effects are therefore $$ \phi_1=\tfrac12[(1-0)+(1-0.04762)]=0.9762,\qquad \phi_2=\tfrac12[(0.04762-0)+(1-1)]=0.02381. $$ They sum to one. Although \(\Theta_2\) does not appear in the output equation, its association with \(\Theta_1\) makes it a weak predictor and earns it a small share of explained variance. This is predictive allocation, not evidence that \(\Theta_2\) changes \(Y\).

Computing Effects for Several Parameters

For \(d\) parameters, exact calculation requires coalition values for subsets of the inputs. There are \(2^d\) subsets, so computing every value can become expensive. A practical approximation samples parameter orderings, estimates the increase in \(v\) at each step, and averages those increases. To estimate a coalition value, one needs the conditional expectation of the model output under the joint distribution. Depending on the available model evaluations and the input distribution, this may be estimated by conditional simulation, a fitted conditional-mean model, or another suitable method. The estimator must respect the joint distribution; independently resampling excluded parameters would generally change the conditional expectation being estimated.

Worked Example: Three-Parameter Allocation from Coalition Values

Suppose a hypothetical QSP analysis estimates the following normalized coalition values under its joint parameter distribution. The output has positive finite variance, the empty coalition has value 0, and all three parameters together have value 1. The values are rounded estimates from conditional calculations.

CoalitionEstimated value
Empty set0.00
\(\{1\}\), \(\{2\}\), \(\{3\}\)0.50, 0.30, 0.10
\(\{1,2\}\), \(\{1,3\}\), \(\{2,3\}\)0.70, 0.65, 0.50
\(\{1,2,3\}\)1.00

For three parameters, the weights on subsets of sizes 0, 1, and 2 are respectively \(1/3\), \(1/6\), and \(1/3\). Parameter 1 therefore has effect $$ \phi_1=\tfrac13(0.50)+\tfrac16[(0.70-0.30)+(0.65-0.10)] +\tfrac13(1.00-0.50)=0.4917. $$ For parameters 2 and 3, the same calculation gives $$ \phi_2=\tfrac13(0.30)+\tfrac16[(0.70-0.50)+(0.50-0.10)] +\tfrac13(1.00-0.65)=0.3167, $$ $$ \phi_3=\tfrac13(0.10)+\tfrac16[(0.65-0.50)+(0.50-0.30)] +\tfrac13(1.00-0.70)=0.1917. $$ The effects sum to 1.0001 because the coalition values and results are rounded; before rounding, the efficiency identity gives a sum of exactly 1. Under this input distribution and output definition, parameter 1 receives the largest share of explained variance. The ranking is specific to the stated distribution and model output.

Assumptions, Diagnostics, and Alternatives

A Shapley analysis requires a clearly defined output, a joint distribution for the parameters, and a finite, nonzero output variance. In a calibrated QSP model, the joint distribution might represent posterior uncertainty or another explicitly chosen uncertainty distribution. The resulting effects describe sensitivity under that distribution; changing the distribution can change the allocation. The output must also be defined precisely, such as a concentration at a specified time or a summary of a simulated trajectory.

For a deterministic model and a fixed output definition, the conditional quantities are mathematically well-defined under the joint distribution. In computation, however, estimated coalition values can be noisy or biased. Check that estimates are stable across independent simulation batches or larger samples, that the values stay within the theoretical range from 0 to 1, and that they are approximately nondecreasing as parameters are added. Exact coalition values are nondecreasing because conditional expectation explains at least as much variance when more information is available; finite-sample estimates may violate this pattern. Also check that estimated effects are near nonnegative and sum to one before rounding. Material violations can signal inadequate conditional modeling, Monte Carlo error, or a mismatch between sampling and the intended joint distribution.

The earlier material on first- and total-order Sobol indices concerns a decomposition particularly suited to independent inputs. Applying standard independent-input Sobol interpretations to dependent parameters is not a substitute for conditioning on the joint distribution: shared information can make the usual separation into individual and interaction contributions inappropriate. Shapley effects instead provide a symmetric allocation of explained variance, but they do not remove the dependence or reveal a unique “true” division of shared information. Regression coefficients answer a different question, often depending on the chosen model form and scale; they are not generally variance allocations.

The main computational limitation is cost. Exact Shapley calculation requires many coalition values, and conditional expectations may be difficult to estimate in regions with little probability or poor overlap. Sampling orderings reduces the number of coalitions evaluated, but introduces additional Monte Carlo error. Report the input distribution, output definition, conditional-expectation method, number of orderings or model evaluations, and uncertainty or stability checks. A Shapley ranking should not be treated as a universal biological ranking or a causal claim.

Common pitfall: “The parameters are correlated, so I can calculate their independent-input Sobol indices and interpret them as separate shares.” This reasoning ignores the information one parameter provides about another. Use conditional expectations under the joint distribution for Shapley effects, and interpret each effect as an order-averaged allocation of predictive variance, not as an isolated or causal contribution.

Check Your Understanding

Answer each question using the definitions and calculations above. Brief solutions follow the questions.

  1. What does \(v(S)\) measure, and why does it depend on the joint distribution when inputs are correlated?
  2. For two parameters with \(v(\{1\})=0.6\), \(v(\{2\})=0.4\), and \(v(\{1,2\})=1\), calculate both Shapley effects.
  3. Why can a parameter receive a nonzero Shapley effect even if it does not appear explicitly in the output equation?
  4. Name two checks that can reveal problems in an estimated Shapley analysis.

Solutions. (1) \(v(S)\) is the fraction of output variance explained by the conditional mean given the parameters in \(S\); dependence affects the conditional distributions of excluded parameters. (2) \(\phi_1=\tfrac12(0.6)+\tfrac12(1-0.4)=0.6\), and \(\phi_2=\tfrac12(0.4)+\tfrac12(1-0.6)=0.4\); they sum to one. (3) It may predict an included parameter through dependence, providing information about the output. (4) Examples include checking stability across simulation batches, checking approximate monotonicity of coalition values, and checking that effects are approximately nonnegative and sum to one.

Clinical Biostats