How can I plot only one fitted variable after using sbiofit?

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!

Accepted Answer

Jeremy Huard
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:
predci_plot.png
Would this work for you?
Best,
Jérémy
  3 Comments
Jim Bosley
Jim Bosley on 30 Aug 2019
Jérémy
Thanks - I think your method (plot the data with biotrells, and then plot just the drug) is what I wanted. So I could say
figure('Name','Drug');
t = sbiotrellis(gmidata, 'ID', 'Time', 'DrugObs', 'Marker', '+','LineStyle','none');
plot(t, simdat,'', '', 'Drug');
figure('Name','Prot')
t = sbiotrellis(gmidata, 'ID', 'Time', 'ProtObs', 'Marker', '+','LineStyle','none');
plot(t, simdat,'', '', 'Prot');
Of course, I'd already written a script to take the plot from plot(fitcon), make two copies, and to delete one variable from one new copy, and the other varaible from the other copy. Ridiculously complex compared to the above!
That said, I think that the plot commands re-run the simulations, right? So my crude script may be faster?
I will emphasize one thing. If you plot all varaibles on the same trellis, difference in scale can obscure a poor fit in the variable with the smaller value. In my case, I had a gross error in the fitting (used the wrong model variable), but the fit looked fine until I separated them. Strongly recommend a flag of some sort like
plot(simdata,'PlotGrouping','together'); % or
plot(simdata,'PlotGrouping','separate');
Jeremy Huard
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

Sign in to comment.

More Answers (1)

Priya Moorthy
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
  2 Comments
Jim Bosley
Jim Bosley on 3 Sep 2019
Thanks, Priya,
I've heard that one goal of a current programming language is that there is one (and only one!) way of doing things. Can't say that about Matlab!
So I set up a script that takes a two-variable biotrellis plot that contains both observed and predicted values for each variable, and it makes two additional plots, sp variables.
It's crude, but it's easy for me to use. I click on the two-variable plot that plot(fitcon) creates, and run the script and I get two biotrellis plots, one with each variable.
I'm copying child objects of the sbiotrellis plots. Another way to do this would be to get the variable names and data object names from the original sbiotrellis plot and then run the commands you very kindly outlined above.
Feel free to criticize the following:
fh = gcf; % This is the current figure - a simdata plot with two variables
fhc = fh.Children; % all children of the figure
fhp = fh.Position; % get position and size
fh2 = figure('Name','1st Variable'); % Create figure for var 1
copyobj(fhc,fh2); % Copy plot(fitcon) to new fig
fh2c = fh2.Children; % variable contains new fig children
fh2.Position = fhp; % Make new fig same size as old
fh3 = figure('Name','2nd Variable'); % Create figure for var 2
copyobj(fhc,fh3); % Copy plot(fitcon) to new fig
fh3c = fh3.Children; % variable contains new fig children
fh3.Position = fhp; % Make new fig same size as old
nfhc = numel(fh2c); % Number of children of the figs
Ymax = 1.0e-5; % Preset
Ymin = 0; % Preset
for i = 1:nfhc
if strcmp(class(fh2c(i)),'matlab.graphics.axis.Axes') % this child == axes?
fh2cic = fh2c(i).Children; % Get children of figure 2 child
if ~isempty(fh2cic) % If no children, subplot is blank
fh3cic = fh3c(i).Children; % Get children of figure 3 child
fh2cic(1).Visible = 'off'; % Turn of obs and pred for one var
fh2cic(3).Visible = 'off';
fh3cic(2).Visible = 'off'; % Turn off obs and pred for other var
fh3cic(4).Visible = 'off';
Ymax = max([fh3cic(1).YData fh3cic(3).YData Ymax]); % Get Max
Ymin = min([fh3cic(1).YData fh3cic(3).YData, Ymin]);% Get Min
end
elseif strcmp(class(fh2c(i)), 'matlab.graphics.illustration.Legend')
lgdtxt = fh2c(i).String; % Get legend text cell array
fh2c(i).String = lgdtxt([1 3]); % Get rid of one var on one plot
fh3c(i).String = lgdtxt([2 4]); % Get rid of other var on other plot
end
end
nYmax = Ymax +0.05 * (Ymax-Ymin); % Kludge a range together for 2nd figure
if Ymin <0
nYmin = Ymin - 0.05*(Ymax-Ymin);
else
nYmin = 0;
end
% Apply range for all 2nd figure axes
for i = 1:nfhc
if strcmp(class(fh3c(i)),'matlab.graphics.axis.Axes')
fh3cic = fh3c(i).Children;
if ~isempty(fh3cic)
fh3c(i).YLim = [Ymin Ymax];
end
end
end
Priya Moorthy
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

Sign in to comment.

Communities

More Answers in the  SimBiology Community

Products

Community Treasure Hunt

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

Start Hunting!