You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
*Global sensitivity analysis* (GSA) is used to study the sensitivity of a function's outputs 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 petab_parameter_fitting), 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 behaviour_optimisation), it can help determine which parameters influence that property. E.g. for a model of a biofuelproducing circuit in a synthetic organism, GSA could determine which system parameters has the largest impact on the total rate of biofuel production.
2
+
*Global sensitivity analysis* (GSA) is used to study the sensitivity of a function's outputs with respect to 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 petab_parameter_fitting), it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
4
+
-[When measuring some system behaviour or property](@ref behaviour_optimisation), 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 have the largest impact on the total rate of biofuel production.
5
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/).
6
+
GSA can be carried out using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package. This tutorial contains 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
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 Catalystrelated 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.
8
+
### [Global vs local sensitivity](@id global_sensitivity_analysis_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
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.
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 ref).
12
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 modelwith an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection others.
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](@ref ref) with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection to others.
15
15
```@example gsa_1
16
16
using Catalyst
17
17
seir_model = @reaction_network begin
@@ -23,10 +23,15 @@ end
23
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:
sol = solve(oprob, Tsit5(); maxiters=100000, verbose=false)
32
+
ps = [:β => p[1], :a => p[2], :γ => p[3]]
33
+
oprob = remake(oprob_base; p = ps)
34
+
sol = solve(oprob; maxiters = 100000, verbose = false)
30
35
SciMLBase.successful_retcode(sol) || return Inf
31
36
return maximum(sol[:I])
32
37
end
@@ -35,62 +40,62 @@ nothing # hide
35
40
Now, GSA can be applied to our `peak_cases` function using GlobalSensitivity's `gsa` function. It takes 3 mandatory inputs:
36
41
- The function for which we wish to carry out GSA.
37
42
- 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.
43
+
- 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 an upper bound for their respective parameter's value.
39
44
40
45
E.g., here we carry out GSA using [Morris's method](https://en.wikipedia.org/wiki/Morris_method):
on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0.0)$ (which corresponds to $β ∈ (0.001,0.1)$, $a ∈ (0.01,1.0)$, $γ ∈ (0.01,1.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.
51
+
on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0.0)$ (which corresponds to $β ∈ (0.001,0.1)$, $a ∈ (0.01,1.0)$, $γ ∈ (0.01,1.0)$). The output of `gsa` varies depending on which GSA approach is used. GlobalSensitivity implements a range of methods for GSA. Below, we will describe the most common ones, as well as how to apply them and interpret their outputs.
47
52
48
53
!!! note
49
54
We should make a couple of notes about the example above:
50
55
- Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_logarithmic_scale), 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. 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 optimization_parameter_fitting_basics), 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 optimization_parameter_fitting_basics), 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.
56
+
- For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base`ODEProblem`, and then apply the [`remake`](@ref ref) function to it in each evaluation of `peak_cases`to generate a problem which is solved for that specific parameter set.
57
+
- Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), 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`.
58
+
- As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases`we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.
54
59
55
60
56
-
## Sobol's methodbased global sensitivity analysis
61
+
## [Sobol's method-based global sensitivity analysis](@id global_sensitivity_analysis_sobol)
57
62
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:
Note: when `Sobol()` is used as the method, the `samples` argument must also be used.
63
68
64
-
Sobol's method computes socalled *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:
69
+
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
70
-`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 secondorder (and the total) Sobol indices are computed, the total order index compounds the information contained in Sobol indices across all orders.
71
+
-`global_sens.S2`: A matrix where element i-j contains the output's sensitivity to simultaneous variations in the i'th and j'th inputs.
72
+
-`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
73
69
-
We can plot the firstorder Sobol indices to analyse their content:
74
+
We can plot the first-order Sobol indices to analyse their content:
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:
79
+
Here, we see that $β$ has a relatively low effect on the peak in infected cases, as compared to $a$ and $γ$. Plotting the total order indices suggests the same:
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:
84
+
GlobalSensitivity implements several versions 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:
It should be noted that when EFAST is used, only the first, and total, order Sobol indices are computed (and not the secondorder ones).
89
+
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
90
86
-
## Morris's methodbased global sensitivity analysis
91
+
## [Morris's method-based global sensitivity analysis](@id global_sensitivity_analysis_morris)
87
92
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):
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:
98
+
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 practice, the following two fields are considered:
94
99
-`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
100
-`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
101
@@ -100,25 +105,26 @@ mean_star_plot = bar(["β" "a" "γ"], global_sens.means_star; labels=["β" "a" "
100
105
variances_plot = bar(["β" "a" "γ"], global_sens.variances; labels=["β" "a" "γ"], title="σ²")
101
106
plot(mean_star_plot, variances_plot)
102
107
```
103
-
Like previously, we note that the peak number of infected cases is more sensitive to $a$ and $γ$ than to $β$.
108
+
As previously, we note that the peak number of infected cases is more sensitive to $a$ and $γ$ than to $β$.
104
109
105
110
!!! 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.
111
+
The syntax for plotting the output using Sobol's and Morris's methods is 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
112
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.
113
+
Generally, Morris's method is computationally less intensive, and has easier to interpret output, as compared to Sobol's method. However, if computational resources are available, Sobol's method is more comprehensive.
109
114
110
115
111
-
## Other global sensitivity analysis methods
116
+
## [Other global sensitivity analysis methods](@id global_sensitivity_analysis_other_methods)
112
117
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
118
114
-
## Global sensitivity analysis for non-scalar outputs
119
+
## [Global sensitivity analysis for non-scalar outputs](@id global_sensitivity_analysis_nonscalars)
115
120
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:
0 commit comments