How to organize plot of sbiosobol results

5 vues (au cours des 30 derniers jours)
Jesse Chao
Jesse Chao le 10 Nov 2021
Hello,
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.
Jesse

Réponse acceptée

Florian Augustin
Florian Augustin le 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);
Best,
Florian
  6 commentaires
Jesse Chao
Jesse Chao le 15 Nov 2021
Hello Florian,
Thank you very much. You make my day with all the suggestions. I appreciate it a lot. They all work perfectly.
I still have a few questions, when I further move on to access my data and plot.
First, I would like to plot my raw data on top of my sobolResults. I found that I couldn't use the hold on and subplot() to access each subplot generate by plot(sobolResults). By executing the code, it overwrites the sobolResults plot. However, I could add my raw data by clicking on each subplot and adding them manually. Is there any code allow me to plot on top of the plot(sobolResults) instead of clicking on the plot?
Second, is there any way that I could have the actual value from the calculation of the classifier by the sbiompgs? Because I have my classifier scripted as abs(trapz(time,PlasmaConc)- AUC_plasma_mean)< 0.5*AUC_plasma_mean. The AUC_plasma_mean is from NCA results which is the mean value of the plasma AUC (correspond to the data 1-5). For some reason, I get all samples are rejected which is quite bizarre to me. Because when I plot my raw data on top of the sobolResults, I expect to have some accepted samples. (the plot is shown below) Though I just want to check with the results of the classifier, if there is anything wrong. Or if you have any idea why this is the case please also let me know. Thank you very much.
By the way, if you think we should discuss through email, please also let me know. I am kind of worry that I drag the conversation in this question too long.
Thank you very much.
Have a nice day.
Best,
Jesse
Florian Augustin
Florian Augustin le 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)", ...
"0.5*AUC_plasma_mean"]);
classifierData = selectbyname(simData, ["lhs", "rhs"]);
% Plot lhs and rhs
sbioplot(classifierData);
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.
Best,
Florian

Connectez-vous pour commenter.

Plus de réponses (1)

Jesse Chao
Jesse Chao le 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,...
'Bounds',modelparambounds,'SamplingMethod',samplingmethod,'NumberSamples',numbersamples,...
'UseParallel',true,'Accelerate',true,'ShowWaitBar',true);
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.
Best,
Jesse
  1 commentaire
Florian Augustin
Florian Augustin le 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.
-Florian

Connectez-vous pour commenter.

Communautés

Plus de réponses dans  SimBiology Community

Catégories

En savoir plus sur Perform Sensitivity Analysis dans Help Center et File Exchange

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by