Global sensitivity analysis - Morris Analysis

New to Julia
1

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) for the Morris analysis and did not manage to get comparable results for means and variances.

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:

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

2

@Vaibhavdixit02

3

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

4

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

5

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

6

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.

8

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

1 Like
9

m.means_star … Thanks!

10

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)?

11

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

2 Likes
12

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:

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).

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!

13

@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

14

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.

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

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]])
15

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.

16

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.

17

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.

18

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.

19

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.

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?

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

20

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?

21

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