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

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!

Answer 1

[Jeremy Huard]

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]

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')

[Jim Bosley]

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]

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

Answer 2

[Priya Moorthy]

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

[Jim Bosley]

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]

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