How to organize plot of sbiosobol results

6 views (last 30 days)
I used to use the GlobalSensitivityAnalysisApp from file exhange. However, I would like to use the script to carry out the sbiosobol in order to be able to define more details for my global sensitivity analysis (sobolResults = sbiosobol(...)). At the very bottom of the GlobalSensitivityAnalysisApp user interface, there is a plot setting for defining how many parameters and obserables per plot. How could I do the same plot setting with code for the result from the sbiosobol?
Please give me any suggestion or advice. Thank you very much.
Have a nice day.

Accepted Answer

Florian Augustin
Florian Augustin on 10 Nov 2021
Hello Jesse,
You can use the name-value arguments "Parameters" and "Observables" of the plot and bar functions to specify which sensitivity inputs (parameters) and sensitivity outputs (observables) are plotted. You can specify the values as names of input parameters / observables, or as numeric indices. The plot and bar function respect the order of the specified parameters / observables / indices, this way you can also sort the Sobol indices.
For example, assuming you have 6 sensitivity inputs / parameters and want to split the plot into two figures, you can run:
% sobolResults = sbiosobol(...);
plot(sobolResults, "Parameters", 1:3);
plot(sobolResults, "Parameters", 4:6);
Florian Augustin
Florian Augustin on 15 Nov 2021
Hi Jesse,
You should be able to plot on top of any GSA plot by turning on "hold" for the axes you want to plot into. You'd first have to get the axes from the plot, hold the existing data of the axes, and then plot on top of it:
% sobolResults = sbiosobol(...);
hFig = plotData(sobolResults); % get figure handle
hAxs = findobj(hFig, "Type", "Axes") % find axes to draw on
selectedAxes = hAxs(1); % assuming axes 1 is the one you want to add content to
hold(selectedAxes, "on"); % hold existing content in axes
plot(selectedAxes, ...); % plot into axes
hold(selectedAxes, "off");
The function sbiompgsa only evaluates the whole classifier and does not store individual evaluations of the left- and right-hand sides. Understanding what is going on during the evaluation of classifiers can be tricky; let me share my workflow how I ususally "debug" the evaluation of classifiers. I typically run a couple of model simulations to see the general ballpark of the ranges of the left- and right-hand sides of the classifier's inequality.
% sobolResults = sbiosobol(...);
% Get SimFunction to reproduce model simulations
simFun = sobolResults.SimulationInfo.SimFunction;
% Get time steps
simTimes = sobolResults.Time;
% Get samples for parameter inputs
paramSamples = sobolResults.ParameterSamples{:,:};
% This is equivalent to writing:
% paramSamples = getSimulationResults(sobolResults, 1);
% If simulating all paramSamples takes too long, you can just select a subset, e.g.
% every 10th sample:
% paramSamples = sobolResults.ParameterSamples{1:10:end, :};
% Run model simulations:
simData = simFun(paramSamples, [], [], simTimes);
% Compute left- and right-hand side of inequality of classifier:
simData = addobservable(simData, ["lhs", "rhs"], ...
["abs(trapz(time, PlasmaConc) - AUC_plasma_mean)", ...
classifierData = selectbyname(simData, ["lhs", "rhs"]);
% Plot lhs and rhs
SimBiology's SimFunctions and how to use them is documented here. The above sample code should give you a spaghetti plot of values of the left- and right-hand side of the classifiers for a range of parameter samples. I assume in your case AUC_plasma_mean is a constant scalar parameter, so you should be see from the evaluation of all lhs evaluations are above, below, or distributed around the rhs. This should match the evaluation of the classifier, i.e. the property SupportHypothesis, on the results object returned from sbiompgsa.
No worries about having a public conversation on this tread. You raise very interesting questions other SimBiology users may also be interested in. However, feel free to send me an email via my community profile if you have model specific questions. I know many users are unable to share model details publicly.
I hope this helps.

Sign in to comment.

More Answers (1)

Jesse Chao
Jesse Chao on 18 Nov 2021
Dear Florian,
Thank you very much for your expertise and patient. The plotting works perfectly. I appreciate your example code.
For the mpgsa, I try to use the debug method you suggested. However, I found something strange which is that the plot doesn't show any results after using the simFun method you suggested. Thus, I further checked with the simData.Data after the simFun(). It showed all the observable column values are zero and the observable we added to debug the classifiers are NaN. On the other hand, I also went back to check the sobolResults.SimulationInfo.SimData.Data and they look normal.
Thus, I am not sure where could be the problem.
It seems that is the simFun = sobolResults.SimulationInfo.SimFunction didn't work properly. I am not sure is it because I have a lot of parameters and variants that need to be applied to my model. When I ran sbiosobol I code it like this.
sobolResults = sbiosobol(model,modelParamNames,OutputName,'Variants',variants,'Doses',doses,'OutputTimes',outputtimes,...
Should I commit all the variants and parameters to the model first and do the sbiosobol() after the commit()? Do you think this could be the reason why simData = simFun(paramSamples, [], [], simTimes); gave me all zero in the simData.Data? or do you have any other thoughts on this?
Please let me know what's your thought and suggestion. Thank you very much.
  1 Comment
Florian Augustin
Florian Augustin on 19 Nov 2021
Hi Jesse,
The SimFunction you get from sobolResults already knows about the model and variants you specified in the call to sbiosobol. But I think the problem is that the dosing information is missing. My example assumed that there are no doses. To apply doses, you need to supply them as follows:
% Assuming doses is a row vector of SimBiology RepeatDose / ScheduleDose objects.
% If you are not sure, you can turn your doses into a row vector as follows:
% doses = doses(:)';
simData = simFun(paramSamples, [], getTable(doses), simTimes);
The reason you need to supply the doses is because they describe "runtime" information that you can change between model simulations. All the different ways to apply doses in SimFunctions is documented here (look for the input argument "u").
The first step would be to get the model simulation working, then we can look at the observables for the left- and right-hand sides of the classifier. Feel free to contact me via my profile page.
I hope that helped.

Sign in to comment.


More Answers in the  SimBiology Community


Find more on Scan Parameter Ranges in Help Center and File Exchange




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!