# Global sensitivity analysis - Morris Analysis

**URL:** <https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825>\
**Category:** New to Julia\
**Created:** [September 18, 2020, 7:21am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825 "2020-09-18T07:21:51Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [September 18, 2020, 7:21am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/1 "2020-09-18T07:21:52Z")

</div>

I am new to Julia and I am looking for a way to do Morris analysis on an system of ODEs. I would like to get a analysis as time series. Additionally, I would like to get the tested parameters to use these as a virtual patient cohort. Unfortunately, I did not manage to run the example successfully…

I tested the Lotka-Volterra Global Sensitivities Analysis example ([Global Sensitivity Analysis · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/analysis/global_sensitivity/)) for the Morris analysis and did not manage to get comparable results for means and variances.

```julia
using DiffEqSensitivity, Statistics, OrdinaryDiffEq, Plots

function f(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2] #prey
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2] #predator
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,u0,tspan,p)
t = collect(range(0, stop=10, length=200))

f1 = function (p)
  prob1 = remake(prob;p=p)
  sol = solve(prob1,Tsit5();saveat=t)
  [mean(sol[1,:]), maximum(sol[2,:])]
end

m = gsa(f1,Morris(total_num_trajectory=1000,num_trajectory=150),[[1,5],[1,5],[1,5],[1,5]])

```

Since this is a copy of the given example I do not unterstand why my output looks like that:

```julia
m.means

2×4 Array{Float64,2}:
 0.0186182 -0.0224811 0.345511 -0.392407
 0.716902 -1.5517 -0.0182938 0.0223239

```

Even though I think that getting a 2x4 output for a system of 2 equations and 4 parameters do make sense.  
Additionally, I am wondering about negative outputs…

Thanks,  
S

---

<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:** [September 18, 2020, 3:34pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/2 "2020-09-18T15:34:49Z")

</div>

@Vaibhavdixit02

---

<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:** [September 20, 2020, 4:30pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/3 "2020-09-20T16:30:58Z")

</div>

> a 2x4 output for a system of 2 equations and 4 parameters do make sense.

Yes, this looks correct, I think the docs example might not have been updated the last time. If you want to create a PR for it feel free, or else I’ll fix it.

> Additionally, I am wondering about negative outputs…

Right, yeah you shouldn’t get negative values but I think it might just be fixed by increasing the `num_trajectory` since 150 is quite low. Try something like `total_num_trajectory=5000,num_trajectory=1000`

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [September 21, 2020, 7:19am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/4 "2020-09-21T07:19:31Z")

</div>

Thank you for your reply!

I increased the values even further (`total_num_trajectory=100000,num_trajectory=10000`) but I still get negative values in the second and the last column of `m.means`

---

<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:** [September 21, 2020, 1:07pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/5 "2020-09-21T13:07:14Z")

</div>

Does the value remain consistent on repeated runs? Btw it might be better to move this to slack or GitHub

---

<div class="post-metadata">

**Author:** ![platawiec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/platawiec/32/31914_2.png) [@platawiec](https://discourse.julialang.org/u/platawiec)\
**Post date:** [September 21, 2020, 1:49pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/6 "2020-09-21T13:49:37Z")

</div>

You may be looking for `m.means_star`, which provides the absolute value of the mean effect. The mean effect, given by `m.mean`, can be negative.

I think `?Morris` should have a doc string.

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [September 21, 2020, 2:29pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/8 "2020-09-21T14:29:38Z")

</div>

No, it does not, but negative values appear in the same positions - at least in all the cases I have looked at.

However, I should have had looked at `m.means_star`

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [September 21, 2020, 2:42pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/9 "2020-09-21T14:42:45Z")

</div>

`m.means_star` … Thanks!

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [September 22, 2020, 6:35am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/10 "2020-09-22T06:35:14Z")

</div>

One last questions… Does anyone has a suggestion how to get the trajectory of the sensitivity (as it is shown in one of the older documentations)?

---

<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:** [September 28, 2020, 7:47am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/11 "2020-09-28T07:47:09Z")

</div>

Do you mean the sensitivity at all time points of the solution?

In that case instead of returning the `[mean(sol[1,:]), maximum(sol[2,:])]` from `f1` just return the `sol`

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [October 2, 2020, 12:30pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/12 "2020-10-02T12:30:37Z")

</div>

Thank you! I figures this one out.

One last thing - I am sorry for all the questions - I just want to make sure that I understand the output correct. When I consider 4 equations and 7 parameters for m.means\_star for example I will get a output like this:

```julia
m.means_star
1×7 Array{Array{Float64,2},2}:
 [0.0 0.00187342 … 4.13169e-7 4.203e-7; 0.0 6.49071 … 0.00874392 0.00818493; 0.0 0.863008 … 0.951201 0.951201; 0.0 1.52598 … 0.108425 0.1084] … [0.0 2645.61 … 274.499 274.465; 0.0 7823.24 … 810.615 809.859; 0.0 15176.7 … 1.1333e7 1.1331e7; 0.0 5.26665e5 … 3.57822e8 3.57762e8]

```

where each of the 7 array represent one parameter (in the same order as I have initialized them).

```julia
m.means_star[1]
4×10980 Array{Float64,2}:
 0.0 0.00187342 6.56241e-6 5.7061e-6 3.85373e-6 3.80246e-6 3.84851e-6 … 4.13261e-7 4.12201e-7 4.10419e-7 4.24247e-7 4.13169e-7 4.203e-7
 0.0 6.49071 0.0839967 0.0836014 0.0822764 0.0856696 0.0819229 0.00880887 0.00884013 0.00894719 0.00791318 0.00874392 0.00818493
 0.0 0.863008 0.863646 0.863671 0.863698 0.863724 0.86375 0.951198 0.951199 0.951199 0.9512 0.951201 0.951201
 0.0 1.52598 1.52665 1.52627 1.52589 1.52552 1.52514 0.108524 0.108499 0.108474 0.10845 0.108425 0.1084

```

And each of this arrays has the sensitivity at each time point and the order of the time arrays correspond to the order of the equations in my system.  
Right? Because my output still looks not as I would have expected it to be.

Thank you!

---

<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 3, 2020, 4:30am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/13 "2020-10-03T04:30:55Z")

</div>

@SoFi I think your understanding of the output seems to be correct. I can not comment of the correctness without knowing more about what seems off to you, it might be much more clearer if I could see the code you are running

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [October 5, 2020, 11:47am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/14 "2020-10-05T11:47:56Z")

</div>

That’s my model. And one thong that seams off for me it that I could have expected from my gut feeling and previous analysis of the system that the p[1] should have the biggest effect on equation 2 but the Julia results indicate something way different.

```julia
VKid = 0.386
VBlood = 4.7
VOrg1 = 1.71
VOrg2 = 0.2

p = [10000.0,2500.0,1.0,1000.0,200.0,0.001,0.1]
u0 = [800000000.0/VBlood,0.0,0.0,0.0]

function CisPlatin(du,u,p,t)
    du[1] = - p[2]*u[1] + p[3]*u[2]*(VKid/VBlood) - p[4]*u[1] + p[6]*u[3]*(VOrg1/VBlood) - p[5]*u[1] + p[7]*u[4]*(VOrg2/VBlood)  
    du[2] = + p[2]*u[1]*(VBlood/VKid) - p[3]*u[2] - p[1]*u[2] 
    du[3] = - p[6]*u[3] + p[4]*u[1] *(VBlood/VOrg1)
    du[4] = - p[7]*u[4] + p[5]*u[1] *(VBlood/VOrg2)
end

```

```julia
include("CisModel.jl")

tspan = (0.0,30.0)
t = collect(range(0, stop=30, length=30*366))
prob = ODEProblem(CisPlatin,u0,tspan,p)

f1 = function (p)
    prob1 = remake(prob;p=p)
    sol = solve(prob1,Tsit5();saveat=t)
    sol
end

a = 0.9
b = 1.1

m = gsa(f1,Morris(total_num_trajectory=100000,num_trajectory=10000),[[10000.0*a,10000.0*b],[2500.0*a,2500.0*b],[1.0*a,1.0*b],[1000.0*a,1000.0*b],[200.0*a,200.0*b],[0.001*a,0.001*b],[0.1*a,0.1*b]])

```

---

<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, 7:30pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/15 "2020-10-05T19:30:20Z")

</div>

I was hoping to reply with something more concrete but I couldn’t get your `gsa` call to finish for me today, I’ll try it again tomorrow.

Meanwhile, did you try to confirm your intuition/morris result with a different GSA method like Sobol? I would suspect it might be happening due to interaction effects which are not easily captured intuitively.

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [October 6, 2020, 1:06pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/16 "2020-10-06T13:06:08Z")

</div>

Yes, previously I worked in R and performed Morris and Sobol analyses for the same system. I assume that the overall result should be the same regardless if I perform it in Julia or R.  
I am going to perform the Sobol analysis in Julia and then compare all my results again.

---

<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, 3:49pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/17 "2020-10-06T15:49:09Z")

</div>

We do test against values from Morris so I would be very surprised if we see a difference. It might be worth checking the function definition to ensure we are comparing the same things.

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [October 7, 2020, 9:18am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/18 "2020-10-07T09:18:01Z")

</div>

Which function definition you mean?  
I checked that it is the same ODE system by comparing the simulation results from R and Julia. They were identical.

I also performed a Sobol analysis in Julia and the results are comparable to the results I got before - but very different from the Morris analysis in Julia.

I am not trying to imply that there is something wrong with the method. I just can not make sense of the results. And I am sorry for that.

---

<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 7, 2020, 9:42am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/19 "2020-10-07T09:42:35Z")

</div>

> I am not trying to imply that there is something wrong with the method. I just can not make sense of the results. And I am sorry for that.

No worries, I am just trying to get a better understanding of the situation as well.

> [@SoFi](#):
>
> Which function definition you mean?  
> I checked that it is the same ODE system by comparing the simulation results from R and Julia. They were identical.

I mean the function we are running the GSA on. Were you were running Morris on the entire time series of the solution in R as well? Are the timepoints you are trying to compare (at which the solutions are being recorded) the same in both R and Julia?

> [@SoFi](#):
>
> I also performed a Sobol analysis in Julia and the results are comparable to the results I got before - but very different from the Morris analysis in Julia.

That’s a little reassuring, was Sobol run for the entire ODE solution as well or on some reduction of the solution?

---

<div class="post-metadata">

**Author:** ![SoFi](https://avatars.discourse-cdn.com/v4/letter/s/5f8ce5/32.png) [@SoFi](https://discourse.julialang.org/u/SoFi)\
**Post date:** [October 7, 2020, 11:25am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/20 "2020-10-07T11:25:31Z")

</div>

The time array were of the same length (from 0 to 30 by 30/366). And I compared the time curves. I was not exacting to get identical results but I thought that the results should be comparable in a qualitative way like which parameters do have a big effect which do not.

Sobol was running for the entire ODE solution, because I want to study if the effect of the parameter gets less important over time.

Would it helpful if I share the plots of my results?

---

<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 9, 2020, 8:28pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825/21 "2020-10-09T20:28:20Z")

</div>

Yeah it should help to see both Sobol and Morris results’ plots to see what’s happening

[Next page](https://discourse.julialang.org/t/global-sensitivity-analysis-morris-analysis/46825.md?page=2)