|
| 1 | +# [Global Sensitivity Analysis](@id global_sensitivity_analysis) |
| 2 | +*Global sensitivity analysis* (GSA) is used to study the sensitivity of a function with respect of its input [^1]. Within the context of chemical reaction network modelling it is primarily used for two purposes: |
| 3 | +- [When fitting a model's parameters to data](@ref TBA_AT_LATER_STAGE), it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters does, and does not, affects the model's fit to the data. This can be used to identify parameters that are less relevant for the observed data. |
| 4 | +- [When measuring some system behaviour or property](@ref TBA_AT_LATER_STAGE), it can help determine which parameters influence that property. E.g. for a model of a biofuel producing circuit in a synthetic organism, GSA could determine which system parameters has the largest impact on the total rate of biofuel production. |
| 5 | + |
| 6 | +GSA can be carried out using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package. This tutorial contain a brief introduction of how to use it for GSA on Catalyst models, with [GlobalSensitivity providing a more complete documentation](https://docs.sciml.ai/GlobalSensitivity/stable/). |
| 7 | + |
| 8 | +#### Global vs local sensitivity |
| 9 | +A related concept to global sensitivity is *local sensitivity*. This, rather than measuring a function's sensitivity (with regards to its inputs) across its entire (or large part of its) domain, measures it at a specific point. This is equivalent to computing the function's gradients at a specific point in phase space, which is an important routine for most gradient-based optimisation methods (typically carried out through [*automatic differentiation*](https://en.wikipedia.org/wiki/Automatic_differentiation)). For most Catalyst related functionalities, local sensitivities are computed using the [SciMLSensitivity.jl](https://github.com/SciML/SciMLSensitivity.jl) package. While certain GSA methods can utilise local sensitivities, this is not necessarily the case. |
| 10 | + |
| 11 | +While local sensitivities are primarily used as a subroutine of other methodologies (such as optimisation schemes), it also has direct uses. E.g., in the context of fitting parameters to data, local sensitivity analysis can be used to, at the parameter set of the optimal fit, [determine the cost function's sensitivity to the system parameters](@ref TBA_AT_LATER_STAGE). |
| 12 | + |
| 13 | +## Basic example |
| 14 | +We will consider a simple [SEIR model of an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). This is an expansion of the classic SIR model with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection others. |
| 15 | +```@example gsa_1 |
| 16 | +using Catalyst |
| 17 | +seir_model = @reaction_network begin |
| 18 | + 10^β, S + I --> E + I |
| 19 | + 10^a, E --> I |
| 20 | + 10^γ, I --> R |
| 21 | +end |
| 22 | +``` |
| 23 | +We will study the peak number of infected cases's ($max(I(t))$) sensitivity to the system's three parameters. We create a function which simulates the system from a given initial condition and measures this property: |
| 24 | +```@example gsa_1 |
| 25 | +using OrdinaryDiffEq |
| 26 | +u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0] |
| 27 | +function peak_cases(p) |
| 28 | + oprob = ODEProblem(seir_model, u0, (0.0, 10000.0), p) |
| 29 | + sol = solve(oprob, Tsit5(); maxiters=100000, verbose=false) |
| 30 | + SciMLBase.successful_retcode(sol) || return Inf |
| 31 | + return maximum(sol[:I]) |
| 32 | +end |
| 33 | +nothing # hide |
| 34 | +``` |
| 35 | +Now, GSA can be applied to our `peak_cases` function using GlobalSensitivity's `gsa` function. It takes 3 mandatory inputs: |
| 36 | +- The function for which we wish to carry out GSA. |
| 37 | +- A method with which we wish to carry out GSA. |
| 38 | +- A domain on which we carry out GSA. This is defined by a vector, which contains one two-valued Tuple for each parameter. These Tuples contain a lower and a upper bound for their respective parameter's value. |
| 39 | + |
| 40 | +E.g., here we carry out GSA using [Morris's method](https://en.wikipedia.org/wiki/Morris_method): |
| 41 | +```@example gsa_1 |
| 42 | +using GlobalSensitivity |
| 43 | +global_sens = gsa(peak_cases, Morris(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]) |
| 44 | +nothing # hide |
| 45 | +``` |
| 46 | +on the domain $10^β ∈ (-3.0,-1.0), 10^a ∈ (-2.0,0.0), 10^γ ∈ (-2.0,0.0)$. The output of `gsa` varies depending on which GSA approach is used. GlobalSensitivity implements a range of methods for GSA. Bellow, we will describe the most common ones, as well as how to apply them and interpret their outputs. |
| 47 | + |
| 48 | +!!! note |
| 49 | + We should make a couple of notes about the example above: |
| 50 | + - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref TBA_AT_LATER_STAGE), this is advantageous in the context of inverse problems such as this one. |
| 51 | + - For simplicity, we create a new `ODEProblem` in each evaluation of the `peak_cases` function. For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. As [previously described](@ref TBA_AT_LATER_STAGE), creating a single `ODEProblem` initially, and then using `remake` to modify it in each evaluations of `peak_cases` will increase performance. |
| 52 | + - Again, as [previously described in other inverse problem tutorials](@ref TBA_AT_LATER_STAGE), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters=100000` and `verbose=false` arguments to `solve`. |
| 53 | + - As we have encountered in [a few other cases](@ref TBA_AT_LATER_STAGE), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, in its third argument, we have to ensure that the i'th Tuple corresponds to the parameter bounds of the i'th parameter in the `parameters(seir_model)` vector. |
| 54 | + |
| 55 | + |
| 56 | +## Sobol's method based global sensitivity analysis |
| 57 | +The most common method for GSA is [Sobol's method](https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis). This can be carried out using: |
| 58 | +```@example gsa_1 |
| 59 | +global_sens = gsa(peak_cases, Sobol(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]; samples=500) |
| 60 | +nothing # hide |
| 61 | +``` |
| 62 | +Note: when `Sobol()` is used as the method, the `samples` argument must also be used. |
| 63 | + |
| 64 | +Sobol's method computes so called *Sobol indices*, each measuring some combination of input's effect on the output. Here, when `Sobol()` is used, the *first order*, *second order*, and *total order* Sobol indices are computed. These can be accessed through the following fields: |
| 65 | +- `global_sens.S1`: A vector where the i'th element is the output's sensitivity to variations in the i'th input. |
| 66 | +- `global_sens.S2`: A matrix where element i-j contain the output's sensitivity to simultaneous variations in the i'th and j'th inputs. |
| 67 | +- `global_sens.ST`: A vector where the i'th element is the output's sensitivity to any simultaneous variation of any combination of inputs that contain the i'th input. While only the first and second order (and the total) Sobol indices are computed, the total order index compounds the information contained in Sobol indices across all orders. |
| 68 | + |
| 69 | +We can plot the first order Sobol indices to analyse their content: |
| 70 | +```@example gsa_1 |
| 71 | +using Plots |
| 72 | +bar(["β", "a", "γ"], global_sens.S1; group=["β", "a", "γ"], fillrange=1e-3) |
| 73 | +``` |
| 74 | +Here, we see that $β$ have a relatively low effect on the peak in infected cases, as compared to $a$ and $γ$. Plotting the total order indices suggests the same: |
| 75 | +```@example gsa_1 |
| 76 | +bar(["β", "a", "γ"], global_sens.ST; group=["β", "a", "γ"], fillrange=1e-3) |
| 77 | +``` |
| 78 | + |
| 79 | +GlobalSensitivity implements several version of Sobol's method, and also provides several options. These are described [here](https://docs.sciml.ai/GlobalSensitivity/stable/methods/sobol/). Specifically, it is often recommended to, due to its quick computation time, use the related extended Fourier amplitude sensitivity test (EFAST) version. We can run this using: |
| 80 | +```@example gsa_1 |
| 81 | +global_sens = gsa(peak_cases, eFAST(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]; samples=500) |
| 82 | +nothing # hide |
| 83 | +``` |
| 84 | +It should be noted that when EFAST is used, only the first, and total, order Sobol indices are computed (and not the second order ones). |
| 85 | + |
| 86 | +## Morris's method based global sensitivity analysis |
| 87 | +An alternative to using Sobol's method is to use [Morris's method](https://en.wikipedia.org/wiki/Morris_method). The syntax is similar to previously (however, the `samples` argument is no longer required): |
| 88 | +```@example gsa_1 |
| 89 | +global_sens = gsa(peak_cases, Morris(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]) |
| 90 | +nothing # hide |
| 91 | +``` |
| 92 | + |
| 93 | +Morris's method computes, for parameter samples across parameter space, their *elementary effect* on the output. Next, the output's sensitivity with respect to each parameter is assessed through various statistics on these elementary effects. In practise, the following two fields are considered: |
| 94 | +- `global_sens.means_star` (called $μ*$): Measures each parameter's influence on the output. A large $μ*$ indicates a parameter to which the output is sensitive. |
| 95 | +- `global_sens.variances`: Measures the variance of each parameter's influence on the output. A large variance suggests that a parameter's influence on the output is highly dependent on other parameter values. |
| 96 | + |
| 97 | +We can check these values for our example: |
| 98 | +```@example gsa_1 |
| 99 | +mean_star_plot = bar(["β" "a" "γ"], global_sens.means_star; labels=["β" "a" "γ"], title="μ*") |
| 100 | +variances_plot = bar(["β" "a" "γ"], global_sens.variances; labels=["β" "a" "γ"], title="σ²") |
| 101 | +plot(mean_star_plot, variances_plot) |
| 102 | +``` |
| 103 | +Like previously, we note that the peak number of infected cases is more sensitive to $a$ and $γ$ than to $β$. |
| 104 | + |
| 105 | +!!! note |
| 106 | + The syntax for plotting the output using Sobol's and Morris's methods are slightly different. The reason is that `global_sens.means_star` and `global_sens.variances` (for Morris's method) are 1x3 Matrices, while for Sobol's method, `global_sens.S1` and `global_sens.ST` are length-3 vectors. |
| 107 | + |
| 108 | +Generally, Morris's method is computationally less intensive, and have easier to interpret output, as compared to Sobol's method. However, if computational resources are available, Sobol's method is more comprehensive. |
| 109 | + |
| 110 | + |
| 111 | +## Other global sensitivity analysis methods |
| 112 | +GlobalSensitivity also implements additional methods for GSA, more details on these can be found in the [package's documentation](https://docs.sciml.ai/GlobalSensitivity/stable/). |
| 113 | + |
| 114 | +## Global sensitivity analysis for non-scalar outputs |
| 115 | +Previously, we have demonstrated GSA on functions with scalar outputs. However, it is also possible to apply it to functions with vector outputs. Let us consider our previous function, but where it provides both the peak number of exposed *and* infected individuals: |
| 116 | +```@example gsa_1 |
| 117 | +function peak_cases_2(p) |
| 118 | + oprob = ODEProblem(seir_model, u0, (0.0, 10000.0), p) |
| 119 | + sol = solve(oprob, Tsit5(); maxiters=100000, verbose=false) |
| 120 | + SciMLBase.successful_retcode(sol) || return Inf |
| 121 | + return [maximum(sol[:E]),maximum(sol[:I])] |
| 122 | +end |
| 123 | +nothing # hide |
| 124 | +``` |
| 125 | + |
| 126 | +We can apply `gsa` to this function as previously: |
| 127 | +```@example gsa_1 |
| 128 | +global_sens = gsa(peak_cases_2, Morris(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]) |
| 129 | +nothing # hide |
| 130 | +``` |
| 131 | +however, each output field is now a multi-row matrix, containing one row for each of the outputs. E.g., we have |
| 132 | +```@example gsa_1 |
| 133 | +global_sens.means_star |
| 134 | +``` |
| 135 | +Here, the function's sensitivity is evaluated with respect to each output independently. Hence, GSA on `peak_cases_2` is equivalent to first carrying out GSA on a function returning the peak number of exposed individuals, and then on one returning the peak number of infected individuals. |
| 136 | + |
| 137 | +--- |
| 138 | +## [Citations](@id global_sensitivity_analysis_citations) |
| 139 | +If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following paper to support the authors of the GlobalSensitivity.jl package: |
| 140 | +``` |
| 141 | +@article{dixit2022globalsensitivity, |
| 142 | + title={GlobalSensitivity. jl: Performant and Parallel Global Sensitivity Analysis with Julia}, |
| 143 | + author={Dixit, Vaibhav Kumar and Rackauckas, Christopher}, |
| 144 | + journal={Journal of Open Source Software}, |
| 145 | + volume={7}, |
| 146 | + number={76}, |
| 147 | + pages={4561}, |
| 148 | + year={2022} |
| 149 | +} |
| 150 | +``` |
| 151 | + |
| 152 | +--- |
| 153 | +## References |
| 154 | +[^1]: [Saltelli, A et al. *Global Sensitivity Analysis. The Primer*, Wiley (2008).](http://www.andreasaltelli.eu/file/repository/A_Saltelli_Marco_Ratto_Terry_Andres_Francesca_Campolongo_Jessica_Cariboni_Debora_Gatelli_Michaela_Saisana_Stefano_Tarantola_Global_Sensitivity_Analysis_The_Primer_Wiley_Interscience_2008_.pdf) |
0 commit comments