3 views (last 30 days)

I'm fitting a model to a dataset. I have two variables I'm fitting simultaneously: the therapy, and the quantity of an endogenous protein. Call the model variables Drug and Prot. The corresponding columns in the dataset woudl be DrugObs and ProtObs. Drug and Prot are of different magnitude. I fit the model:

[fitcon, simdat] = sbiofit( m1, gmidata, resmap, estpars, doses, 'UseParallel',true);

If I use plot(fitcon) I get a trellis plot of data and mode for each group (that is, each subject). This is great! However, since the magnitudes are different I'd like to plot the PRED and OBS of each variable in a different plot.

Also, it's useful to look at these things in both log and linear form.

So is there any easy way to say

plot(fitcon,'PREDVarToPlot','Drug', 'OBSVarToPlot', 'DrugObs'); % and

plot(fitcon,'PREDVarToPlot','Prot', 'OBSVarToPlot', 'ProtObs'); %

and also

plot(fitcon,'PREDVarToPlot','Drug', 'OBSVarToPlot', 'DrugObs',YScale,'log'); % and

plot(fitcon,'PREDVarToPlot','Prot', 'OBSVarToPlot', 'ProtObs',YScale,'log'); %

I can do this by going into the plot but its tedious:

fh = gcf;

fhc = fh.Children

nfhc = numel(fhc)

for i = 1:nfhc

test the specific figure child to see if its an axes or a legend

if axes get axes children

fhcc = fhc(i).Childen % These are lines

fhcc(1).Visible = 'off' % This is one way. Any way to delete the line completely?

fhcc(3).Visible = 'off' % fhcc([1 3]) are two of the four plotted variables corresponding to etiher drug, or protein. Could be fhcc([2 4])

elseif legend

castr = fhc(i).String

castr = castr([2 4])

end

Another approach would be to allow the second variable to be plotted using a different (right hand) yscale. I've not tried this in script, but as I understand it I'd have to loop through current axes and create a new Child axes with YAxisLocation = 'right'. Not sure then how I would transfer one set of variables (Drug and DrugObs or Prot and ProtObs) to that new axis. Again, it would be ideal if there was an ability to say:

plot(fitcon, 'LeftAxis', 'Drug','RightAxis','Prot');

Are there any easy ways to do this, or must I create a script function? Thanks!

Jeremy Huard
on 30 Aug 2019

Hi Jim,

thanks for this feedback. It would be a good thing to have.

Have you had a look at sbiopredictionci? It computes the compute confidence intervals for model predictions. It also has a plot method that will create a grid with one column per group and one row per response. This is basically what you are asking for, except that you will also have a plot of the confidence intervals.

It would look like this:

Would this work for you?

Best,

Jérémy

Jeremy Huard
on 30 Aug 2019

By the way, if you prefer without the confidence intervals you can achieve it with a couple of lines of code:

t = sbiotrellis(gmidata, 'ID', 'Time', 'DrugObs',...

'Marker', '+','LineStyle','none');

plot(t, simdat,'', '', 'Drug');

% the following lines are there to have both obs and sim data

% in the same color

gobj = findobj(t.hFig, 'Type','Line');

set(gobj, 'Color','black')

Jeremy Huard
on 3 Sep 2019

Hi Jim,

Thanks for your suggestion.

Just to clarify one thing: the plots commands do not rerun the simulations. They simply add the simulation data contained in simdat to the trellis.

This is in essence what plot(fitcon) does. fitcon (the OptimResults object) saves the simulation data in its property fitted. So, the sbiotrellis + plot commands are more efficient because your code does what they do twice plus some other things to delete the unnecessary lines.

Best,

Jérémy

Priya Moorthy
on 3 Sep 2019

Edited: Priya Moorthy
on 4 Sep 2019

Hi, Jim,

Currently, SimBiology plot functions do not support using multiple y-axes. You can pass a custom plot function to sbiotrellis that will allow you to use different axis scales. You will need a helper function that allows you to use plotting functions like @semilogy with simData objects.

If your original call is:

[fitcon, simdat] = sbiofit( m1, gmidata, resmap, estpars, doses, 'UseParallel',true );

...then the the plotting function calls are as follows, where groupVariableName and independentVariableName are the group and time column names from your experimental data:

To plot linearly:

h = sbiotrellis(expData, groupVariableName, independentVariableName, ‘DrugObs’, ':o');

plot(h, simData, [], 'time', ‘Drug’);

To plot on a log scale:

h = sbiotrellis(gmidata, @semilogy, groupVariableName, independentVariableName, ‘DrugObs’, ':o');

plot(h, simdat, @(sd, xcol, ycol)plotSimDataHelperFcn(sd, xcol, ycol, @semilogy), 'time', 'Drug');

The helper function would look something like this:

function plotSimDataHelperFcn(sd, xcol, ycol, fcnhandle)

% SD is a single SimData object

% XCOL and YCOL are strings or char vectors that refer to states in the

% simData object

% FCNHANDLE must have signature FCNHANDLE(X,Y)

if strcmpi(xcol, 'time')

x = sd.Time;

else

x = sd.selectbyname(xcol).Data;

end

if strcmpi(ycol, 'time')

y = sd.Time;

else

y = sd.selectbyname(ycol).Data;

end

fcnhandle(x,y);

There are some features coming up in a future release to make it easier to plot select variables from simData objects and experimental objects from the desktop interface (you can try it out in the R2019b prerelease, although there will be some changes to the features in the final release).

Best,

Priya

Priya Moorthy
on 4 Sep 2019

Jim,

Your copyobj method certainly works, although it is relying on the order of the lines to decide which ones to turn off and on. You could also delete lines instead of just changing their visibility. My solution is virtually the same as Jérémy's; I just added the function handle argument to enable changing the y-axis scale.

As Jérémy said, using the plot command with SimData objects will not re-run the simulation; however using the plot command with the fit results object may rerun simulations to get the predicted data. So, using sbiotrellis with the SimData objects returned from the fitting function is probably the more efficient solution.

Best,

Priya

More Answers in the SimBiology Community

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

Start Hunting!
## 0 Comments

Sign in to comment.