# How to interpret 2nd order Sobol global sensitivity analysis result?

**URL:** <https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804>\
**Category:** Statistics\
**Created:** [October 5, 2020, 5:00pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804 "2020-10-05T17:00:43Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Janssena](https://avatars.discourse-cdn.com/v4/letter/j/dfb087/32.png) [@Janssena](https://discourse.julialang.org/u/Janssena)\
**Post date:** [October 5, 2020, 5:00pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/1 "2020-10-05T17:00:43Z")

</div>

Hey All,

I used DiffeqSensitivity’s global sensitivity analysis using the Sobol method in order to explain the influence of specific input parameters of a model I’m building. I’m already getting some exciting results from the first order indices, and am also interested in the second order effects.  
I am however unsure how to interpret the second order results.  
My model has 35 input variables and 4 output variables. The resulting second order indices matrix has a size of 595 x 4.  
How do I know which indexes in the second order indices matrix correspond to the interaction between two input variables?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 5, 2020, 5:12pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/2 "2020-10-05T17:12:04Z")

</div>

@Vaibhavdixit02 might want to comment here.

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [October 5, 2020, 6:22pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/3 "2020-10-05T18:22:18Z")

</div>

@Janssena the second order indices proceed as `1-2, 1-3, 1-4...1-n, 2-3, 2-4,...2-n...` i.e. for each parameter starting from the first we take the combination of that parameter with all subsequent parameters. This leads to 35 \choose 2 total combinations which is 595. I hope this clarifies it?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 5, 2020, 9:43pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/4 "2020-10-05T21:43:48Z")

</div>

@Vaibhavdixit02 food for thought: reshape and display it as a matrix?

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [October 6, 2020, 4:30am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/5 "2020-10-06T04:30:57Z")

</div>

Hmm, not sure what the reshape should be? Currently it is stored as rows for dependent variables with the sensitivity for that variable due to the parameter (or combination of parameters in second order) as the columns indexed as per the order they are passed in. I think the op might have the dimensions for their output mentioned the other way round (I checked to ensure I remembered correctly)

---

<div class="post-metadata">

**Author:** ![Janssena](https://avatars.discourse-cdn.com/v4/letter/j/dfb087/32.png) [@Janssena](https://discourse.julialang.org/u/Janssena)\
**Post date:** [October 6, 2020, 5:58am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/6 "2020-10-06T05:58:44Z")

</div>

That is exaclty what I meant yeah! Shape indeed is 4×595. I think it would be most clear if you have a input x input matrix per output variable. That makes it easy to look at all interactions for each input variable and can immediately be visualized using something like a heatmap.

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [October 6, 2020, 4:02pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/7 "2020-10-06T16:02:52Z")

</div>

If you want to take a look at using heat map with the `gsa` output there is an example [https://github.com/PumasAI/PumasTutorials.jl/blob/f97b04da047197643456c61b04deabc0b797892d/tutorials/pkpd/hcvgsa.jmd#L203-L212](https://github.com/PumasAI/PumasTutorials.jl/blob/f97b04da047197643456c61b04deabc0b797892d/tutorials/pkpd/hcvgsa.jmd#L203-L212), it differs a bit since the dispatch for `gsa` in Pumas returns a separate struct but it should hopefully give you enough hint to work it out. You can ping me if you get stuck.

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [November 17, 2021, 10:18am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/8 "2021-11-17T10:18:03Z")

</div>

I was wondering about the second order indices as well.

Using the MWE from the docs

```julia
using GlobalSensitivity, QuasiMonteCarlo

function ishi(X)
    A= 7
    B= 0.1
    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
end

n = 600000
lb = -ones(4)*π
ub = ones(4)*π
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,sampler)

res1 = gsa(ishi,Sobol(order=[0,1,2]),A,B)

```

The sum over all 2nd order indices considerably exceeds 1. Are the 2nd order indices somehow cumulative?

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [November 17, 2021, 11:07am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/9 "2021-11-17T11:07:16Z")

</div>

Second order indices wouldn’t sum to 1, that’s only true for the first order. I am not sure what you mean by cumulative? The definition of second order should make it clear why you can’t expect them to sum to 1, it’s the sensitivities wrt each pair of the parameters, x1&x2, x1&x3…x1&xn, x2&x3, x2&x4 …xn-1&xn. I hope that makes sense 😅

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [November 17, 2021, 11:41am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/10 "2021-11-17T11:41:50Z")

</div>

If the `gsa` function returns sobol indices, then in the case of orthogonal inputs, the sum of all indices should be equal to 1.

See also the ishigami example [here](https://salib.readthedocs.io/en/latest/basics.html).

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [November 17, 2021, 11:45am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/11 "2021-11-17T11:45:12Z")

</div>

```julia
using GlobalSensitivityAnalysis
using Distributions, DataStructures
params = [Symbol("p$i") for i =1:4]
pbounds = [(l,u) for (l,u) in zip(lb,ub)]
data = SobolData(
    params = OrderedDict([p => Uniform(b[1],b[2]) for (b,p) in zip(pbounds,params)]...),
    N = 10000
)
samples = GlobalSensitivityAnalysis.sample(data)
function ishi_batch(X)
    A= 7
    B= 0.1
    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
end
y = ishi_batch(samples')
sob = analyze(data, y)

```

gives  
 ![grafik](https://global.discourse-cdn.com/julialang/original/3X/e/1/e1701093d48379eb16806bb950faa2124d123856.png)

Negative values do not make sense, but this might be avoided by increasing the number of samples, but the direction is correct, compare also with the results from [here](https://salib.readthedocs.io/en/latest/basics.html).

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [November 18, 2021, 1:51pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/12 "2021-11-18T13:51:33Z")

</div>

I don’t see the sum to 1 in the link you shared?

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [November 18, 2021, 1:53pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/13 "2021-11-18T13:53:57Z")

</div>

There seems to be a negative value in the link as well, yes it can be avoided by increasing the number of samples, when it’s so close to 0 it can sometimes turn up -ve, for all measures you can consider it effectively 0 in that case.

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [November 18, 2021, 1:56pm UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/14 "2021-11-18T13:56:25Z")

</div>

If you some up all entries in the first order vector and the second order matrix, you end up around 1:

![grafik](https://global.discourse-cdn.com/julialang/original/3X/e/1/e1701093d48379eb16806bb950faa2124d123856.png)

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [December 15, 2021, 8:49am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/15 "2021-12-15T08:49:20Z")

</div>

Is there an update on this? Shall I open an issue on github linking to this discussion?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 15, 2021, 9:11am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/16 "2021-12-15T09:11:03Z")

</div>

yes please

---

<div class="post-metadata">

**Author:** ![moesphere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moesphere/32/3348_2.png) [@moesphere](https://discourse.julialang.org/u/moesphere)\
**Post date:** [December 15, 2021, 9:20am UTC](https://discourse.julialang.org/t/how-to-interpret-2nd-order-sobol-global-sensitivity-analysis-result/47804/17 "2021-12-15T09:20:06Z")

</div>

Issue opened [here](https://github.com/SciML/GlobalSensitivity.jl/issues/46).