How to organize plot of sbiosobol results

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

Answer 1

[Florian Augustin]

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

[Jesse Chao]

Hello Florian,

Thank you very mcuh. It works perfectly.

Beside the plot setting in the GlobalSensitivityAnalysisApp, there is also a check box for "Reuse Sobol model simulation for MPGSA" in the user interface. How could I execute this command by code? I only know that for the MPGSA is by sbiompgsa()

Could you please provide me with further instructions? Thank you very much.

Best,

Jesse

[Florian Augustin]

Hi Jesse,

you can provide the parameter samples and model simulations from the results object returned from sbiosobol as input arguments to sbiompgsa. To get those samples and simulation results you can use the getSimulationResults function.

For example, you can run sbiompgsa using the two independent sets of parameter samples and associated model simulations from sbiosobol as follows:

% Run Sobol analysis:

% sobolResults = sbiosobol(...);

% Get parameter samples and model simulations:

[samplesTable1,simdata1] = getSimulationResults(sobolResults,1);

[samplesTable2,simdata2] = getSimulationResults(sobolResults,2);

% you can also get all simulation results for use in the MPGSA;

% the documentation for getSimulationResults explains what those samples are.

% Stack the samples tables and simulation results:

samplesTable = [samplesTable1; samplesTable2];

simData = [simdata1; simdata2];

% Example classifiers:

% classfier = "scalarModelResponse <= mean(scalarModelResponse)"; % or

% classfier = "mean(timeDependentModelResponse) <= mean(timeDependentModelResponse, 'all')";

% Run MPGSA:

mpgsaResults = sbiompgsa(simData, samplesTable, classifier);

% Plot results (again assuming at least 6 sens. inputs / parameters):

plot(mpgsaResults, "Parameters", 1:3);

plot(mpgsaResults, "Parameters", 4:6);

I hope this helps.

-Florian

[Jesse Chao]

Hello Florian,

Thank you very much. I appreciate all your precious suggestions. They help a lot.

However, when I try to stake my parameter samples and model simulations as you suggested, my computer freeze due to "System runs out of application memory". The sobolResults I have has the size of 4.89 GB which is quite a large data and please find my computer specifications below.

However, I never have this problem when I was using the GlobalSensitivityAnalysisApp. I used the following code to stake my data.

samplesTable_mpgsa  = [];
simdata_mpgsa       = [];
validRuns_mpgsa     = [];

parfor i = 1:length(sobolResults.ParameterSamples.Properties.VariableNames)+2

    [samplesTable_ans,simdata_ans,validRuns_ans] = getSimulationResults(sobolResults,i);
    samplesTable_mpgsa = [samplesTable_mpgsa; samplesTable_ans];
    simdata_mpgsa      = [simdata_mpgsa     ; simdata_ans     ];
    validRuns_mpgsa    = [validRuns_mpgsa   ; validRuns_ans   ];

    samplesTable_ans = []; % free memory
    simdata_ans      = []; % free memory
    validRuns_ans    = []; % free memory

end

Could you please give me some suggestions to make my code more efficient like the way that the GlobalSensitivityAnalysisApp handle? Thank you very much.

Best,

Jesse

[Florian Augustin]

Hi Jesse,

I see you are using a parfor loop to concatenate the samples table and simulation data. For the parallel execution, the data used within this loop is copied to each worker (even if the parallel pool is local). This is probably too much memory.

I assume you tried parfor to speed up the concatenation of the big data. Instead of using parfor, can you try using cells to avoid the iterative increase of the arrays within the loop?

numSetsOfSamples = length(sobolResults.ParameterSamples.Properties.VariableNames) + 2;

% Avoid iterative dynamic increase of array size using cell

samplesTable_mpgsa = cell(numSetsOfSamples, 1);

simdata_mpgsa = cell(numSetsOfSamples, 1);

validRuns_mpgsa = cell(numSetsOfSamples, 1);

for i = 1:numSetsOfSamples

[samplesTable_mpgsa{i}, simdata_mpgsa{i}, validRuns_mpgsa{i}] = getSimulationResults(sobolResults,i);

end

% Concatenate samples tables and simulation results

samplesTable_mpgsa = vertcat(samplesTable_mpgsa{:});

simdata_mpgsa = vertcat(simdata_mpgsa{:});

validRuns_mpgsa = vertcat(validRuns_mpgsa{:});

I just checked and the Global Sensitivity Analysis app only uses the first set of independent samples for the MPGSA, but there is no reason why you couldn't use at least the first two sets of samples.

Unrelated to your question, but if you ever run into problems saving your GSA results to a mat file, you can specify the mat-file version 7.3 during save; this saves mat-files in the HDF5 format that allows bigger file sizes than the default mat-files.

Best,

Florian

[Jesse Chao]

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]

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

Answer 2

[Jesse Chao]

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

[Florian Augustin]

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