GitHub

Sensitivity Analysis

In this tutorial, we will showcase how to perform global sensitivity analysis of a model. To get started, let's first pull in our modified second order ODE Lorenz equation model from before:

using EasyModelAnalysis

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(D(x)) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

@named sys = ODESystem(eqs)
sys = structural_simplify(sys)

u0 = [D(x) => 2.0,
    x => 1.0,
    y => 0.0,
    z => 0.0]

p = [σ => 28.0,
    ρ => 10.0,
    β => 8 / 3]

tspan = (0.0, 100.0)
prob = ODEProblem(sys, u0, tspan, p, jac = true)
sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 1495-element Vector{Float64}:
   0.0
   0.00014131524176735154
   0.0015544676594408668
   0.012956381792766277
   0.038266700044138006
   0.07852291796751046
   0.12210864065288088
   0.17155375443120538
   0.2311290630460472
   0.30252383872495503
   ⋮
  99.50500396864174
  99.57665667367455
  99.64283640514518
  99.72166966002771
  99.79967562312255
  99.86801143860899
  99.93304884850775
  99.99013919862291
 100.0
u: 1495-element Vector{Vector{Float64}}:
 [2.0, 0.0, 0.0, 1.0]
 [1.9960454103704428, 0.0014132521265988688, 9.986094350108765e-8, 1.0002823510089436]
 [1.956746201469907, 0.01555657123298723, 1.2096178746316906e-5, 1.0030752466310615]
 [1.6563910957343149, 0.13029732534611171, 0.0008464888325952064, 1.023645017393648]
 [1.1011942420041494, 0.3874130558030291, 0.007445311135321615, 1.0582098436702096]
 [0.5557298463522196, 0.7959541011040255, 0.031202921104640697, 1.0901348104540605]
 [0.44950102412404563, 1.229007863529698, 0.07390400544865926, 1.1102092717261387]
 [0.9266726317320947, 1.7062381668292075, 0.14159577334379553, 1.1416871365805596]
 [2.2781005968110417, 2.2773793271385414, 0.25114094267068454, 1.233177108242911]
 [4.867419855792264, 3.014693593518013, 0.44047994481421116, 1.4824570579833207]
 ⋮
 [80.19625687470715, 0.22319498068561278, 6.377998459496999, -10.92101238125288]
 [93.91126047193579, -1.9435112337378833, 5.7144265707275865, -4.580907294954086]
 [92.47053186268388, -2.2128845219633724, 4.966589670343976, 1.653117805222435]
 [78.34378276619968, 0.2736929875090918, 3.6720087645277166, 8.448758815091098]
 [59.61988722700268, 5.207887350906977, 5.364596386122552, 13.834326658988568]
 [42.06608771509131, 6.848721336825014, 10.923029356998175, 17.327760613533872]
 [17.746361864729497, 2.561462421673927, 14.710142641007621, 19.336065831253446]
 [-13.425204706776736, -1.9987153748274389, 12.635789734327103, 19.49627693453924]
 [-19.39944052157736, -2.411052627709995, 11.88863957641005, 19.334495160672596]

On this model, we wish to perform sensitivity analyses. Global sensitivity analysis requires the specification of two things:

  1. Global sensitivity of what? What is the value to be measured from the model that you want to assess the sensitivity of?
  2. Over what space of parameter values? This is a box of potential values for the parameters of the model.

The function get_sensitivity(prob, t, obs, bounds) first takes (t,obs), where t is the time point to take measurements from the model and y is the desired observable to measure. bounds is specified as an array of pairs, which maps parameter symbols to arrays specifying the lower and upper bounds of its parameter range.

Thus for example, let's calculate the sensitivity of y(100) over the parameters (ρ, β) where $\rho \in [0,20]$ and $\beta \in [0,100]$:

pbounds = [ρ => [0.0, 20.0], β => [0.0, 100.0]]
sensres = get_sensitivity(prob, 100.0, y, pbounds; samples = 2000)
Dict{Symbol, Float64} with 5 entries:
  :ρ_first_order    => 0.107475
  :β_first_order    => 0.0326214
  :ρ_total_order    => 0.988595
  :β_total_order    => 0.879228
  :ρ_β_second_order => 0.825519

The output shows values of first_order, second_order and total_order sensitivities. These are quantities that define the nonlinear effect of the variables on the output. The first order values can be thought of as the independent interaction effects, similar to the values of a linear regression for $R^2$, i.e. it's the variance explained by independent effects of a given parameter. The first indices are the scaled variance terms, for example:

ρ_first_order = V[y(100) changing only ρ] / V[y(100)]

where V denotes the variance. I.e., ρ_first_order is the percentage of variance explained by only changing ρ. Being normalized, if the model is linear then ρ_first_order + β_first_order == 1, and thus its total summation tells us the degree of nonlinearity. Our simulation here has the sum of first indices as <0.2, an indication of a high degree of nonlinear interaction between the measured output and the parameters.

The second order indices then say how much can be attributed to changes of combinations of variables. I.e.:

ρ_β_second_order = V[y(100) changing only ρ and β] / V[y(100)]

which thus gives the percentage of variance explained by the nonlinearities of ρ and β combined. These sensitivity functions only output up to the second indices, since there is a combinatorial explosion in the number of terms that need to be computed for models with more parameters. However, in this tutorial there are only two variables, and thus all variance should be explained by just these two parameters. This means that ρ_first_order + β_first_order + ρ_β_second_order should be approximately equal to 1, as all variance should be explained by the linearity or second order interactions between the two variables. Let's check this in action:

sensres[:ρ_first_order] + sensres[:β_first_order] + sensres[:ρ_β_second_order]
0.9656149538655133

This is not exactly equal to 1 due to the numerical error in the integral approximations, but you can see theory in action! (Also, this is a good test for correctness of the implementation).

Now, if you had more than two variables, is there a good way to get a sense of the “sensitivity due to ρ”? Using the first indices is a bad approximation to this value, since it's only the measurement of the independent or linear sensitivity of the output due to ρ. Instead, what one would want to do, is to get the sum of the first order index ρ_first_order, plus all second order effects which include ρ, plus all third order effects that include ρ, plus … all the way to the Nth order for N variables. Surprisingly, this summation can be computed without requiring the computation of all the higher order indices. This is known as the “total order index” of a variable. In a two parameter model, we can see this in action:

sensres[:ρ_first_order] + sensres[:ρ_β_second_order], sensres[:ρ_total_order]
(0.932993563491674, 0.9885950888540156)
sensres[:β_first_order] + sensres[:ρ_β_second_order], sensres[:β_total_order]
(0.8581400072304508, 0.879227712171329)

Thus, the total indices are a good measurement of the relative size of the total effect of each parameter on the solution of the model.

In summary:

  • First order indices showcase the amount of linearity and the direct linear attributions to each variable
  • The second order indices show the linear correlations in the outputs
  • The total indices measure the total effect a given variable has on the variance of the output

and notably, all values are normalized relative quantities.

Thus we can finally use the create_sensitivity_plot function to visualize the field of sensitivity results:

create_sensitivity_plot(prob, 100.0, y, pbounds; samples = 2000)
Example block output

which shows the relative sizes of the values in plots for the first, second, and total index values.