Color inside Confidence Interval after a fitlm

Hi

I'd like to color the area in betwen the Confidence Interval after a fitlm, like shown in the attached pic

I tried patch, area, fill, but I couln't make any of them work.

I think the issue is that the x values are not continuous (it's not 1:19 for example)

I have 2 continuous variables

% x = rand(1,19);  y = x+[1:19];  
d = load('d.txt') ;
x = d(:,1) ; 
y = d(:,2) ; 
mdl = fitlm(x,y)

mdl = 
Linear regression model:
    y ~ 1 + x1

Estimated Coefficients:
                    Estimate        SE         tStat        pValue  
                   ___________    ______    ___________    _________

    (Intercept)    -0.00082259    1.1472    -0.00071705      0.99944
    x1                 -16.629    4.6874        -3.5477    0.0024745


Number of observations: 19, Error degrees of freedom: 17
Root Mean Squared Error: 5
R-squared: 0.425,  Adjusted R-Squared: 0.392
F-statistic vs. constant model: 12.6, p-value = 0.00247

mdl.plot

thanks a lot

Answer 1

[dpb]

From the point of generating the linear model incorporating both forms of either only 3 or 4 lines.

hL=lm.plot;                               % returns array of line objects
hCI=findobj(hL,'-regexp','DisplayName','.bounds');  % find the labelled bounds one
x=hCI.XData;                              % retrieve x,y
y=hCI.YData;
if numel(hL)==3                           % newer version with one line for bounds
  x=x(isfinite(x));                       % keep only the finite data
  y=y(isfinite(y));
  nHalf=numel(x)/2;                       % it's all one line object in two halves...
  x=[x(1:nHalf) flip(x(nHalf+1:end))];    % reflect the second half for fill() to
  y=[y(1:nHalf) flip(y(nHalf+1:end))];    % have an outline shape to fill in
elseif numel(hL)==4                       % older version with two lines for bounds
  hCI=findobj(hL,'DisplayName','');       % find the unlabelled bounds one
  x=[x flip(hCI.XData)];                  % reflect and add the second line for fill()
  y=[y flip(hCI.YData)];
else
  error('Unexpected number of lines from linear fit plot method.  Aborting.')
end
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String(end)=[];

[dpb]

Kept nagging at me...an alternate way of splitting the two halves w/o explicitly calculating the half length for newer version is to reshape the two lines data vector into two columns. (Nota bene have to use two columns to split in half or will be taking alternate values in pairss.)

data=readtable('d.txt');
data.Properties.VariableNames={'X','Y'};
lm=fitlm(data);
hL=lm.plot;                               % returns array of line objects
hCI=findobj(hL,'-regexp','DisplayName','.bounds');  % find the labelled bounds one
x=hCI.XData;                              % retrieve x,y
y=hCI.YData;
switch numel(hL)
  case 3                                  % newer version with one line for bounds
    x=reshape(x(isfinite(x)),[],2);       % keep only the finite data and break into half
    y=reshape(y(isfinite(y)),[],2);       % reshape into two columns instead of one long row
    x=[x(:,1); flip(x(:,2))];             % reflect the second half for fill() to
    y=[y(:,1); flip(y(:,2))];             % have an outline shape to fill in
  case 4                                  % older version with two lines for bounds
    hCI=findobj(hL,'DisplayName','');     % find the unlabelled bounds one
    x=[x flip(hCI.XData)];                % reflect and add the second line for fill()
    y=[y flip(hCI.YData)];
  otherwise
    error('Unexpected number of lines from linear fit plot method.  Aborting.')
end
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String(end)=[];

Answer 2

[dpb]

You could sort your data first, but the line boundaries would be jagged and not follow the plot object lines closely. The easiest thing is to get the data from the plot that fitlm already drew for you.

data=readtable('d.txt');
data.Properties.VariableNames={'X','Y'};
lm=fitlm(data);
lm.plot
hAx=gca;
hCI=hAx.Children(1);  % discover this by looking; it isn't documented which line is which
x=hCI.XData(isfinite(hCI.XData));       % keep only the finite data
y=hCI.YData(isfinite(hCI.YData));
nHalf=numel(x)/2;                       % it's all one line object in two halves...
x=[x(1:nHalf) flip(x(nHalf+1:end))];    % reflect the second half for fill() to
y=[y(1:nHalf) flip(y(nHalf+1:end))];    % have an outline shape to fill in
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String=hLg.String(1:end-1);

The "trick" is to know (or find out by exploring) that the confidence intervals are one line object with a NaN between the two sections so they won't be joined on the plot and that fill won't fill in anything that has any NaN elements. So, first remove that element and then flip the second half of the two arrays to make an area instead of two lines. fill will recognize and close the starting end by joinng the first and last x values.

You can adjust color and opacity to suit tastes.....

ADDENDUM

A better route is to save the lines array from the call to the plot() method. Interestingly, they're in the opposite order of the axes children array. One can be a little more generic to avoid that issue

figure
hL=lm.plot;                             % returns array of line objects
hCI=findobj(hL,'-regexp','DisplayName','.bounds');  % find the bounds one
x=hCI.XData(isfinite(hCI.XData));       % keep only the finite data
y=hCI.YData(isfinite(hCI.YData));
nHalf=numel(x)/2;                       % it's all one line object in two halves...
x=[x(1:nHalf) flip(x(nHalf+1:end))];    % reflect the second half for fill() to
y=[y(1:nHalf) flip(y(nHalf+1:end))];    % have an outline shape to fill in
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String(end)=[];

[Pat]

Hi

thanks a lot!!

it works fine but only with Matlab 2025

I also have Matlab 2023 and 2021 and it does not work!! I'm usually using 2021, but I'll switch!!!

best

[dpb]

I had just done it interactively here on the forum...in R2021b,

>> hL=lm.plot
hL = 
  4×1 Line array:

  Line    (Data)
  Line    (Fit)
  Line    (Confidence bounds)
  Line
>> 

The two bounds are separate lines, not one with the NaN in the middle. The second is unlabeled so there aren't duplicated legend entries.

The fix is simple enough, use the data from the two line handles but since are separate lines don't divide the data arrays in half, just flip one of the two.

...
hL=lm.plot;                             % returns array of line objects
hCI=findobj(hL,'-regexp','DisplayName','.bounds');  % find the labelled bounds one
xL=hCI.XData;                           % it's the lower bound line
yL=hCI.YData;
hCI=findobj(hL,'DisplayName','');       % find the unlabelled bounds one
xU=hCI.XData;                           % it's the upper bound line
yU=hCI.YData;
x=[xL flip(xU)];                        % reflect the second for fill() to
y=[yL flip(yU)];                        % have an outline shape to fill in
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String(end)=[];

If you use multiple versions, probably best to write a little function to branch depending on the version being used. Or, you could count the number of lines returned from the plot() method and avoid having to discover which release contained the change.

ADDENDUM

A slightly more efficient way would be to use x,y variables in place in lieu of so many temporaries since there are no modifications needed in this case it's easy to write the result vector directly.

...
hL=lm.plot;                             % returns array of line objects
hCI=findobj(hL,'-regexp','DisplayName','.bounds');  % find the labelled bounds one
x=hCI.XData;                            % start the boundaries vectors
y=hCI.YData;
hCI=findobj(hL,'DisplayName','');       % find the unlabelled bounds one
x=[x flip(hCI.XData)];                  % reflect the second for fill() to
y=[y flip(hCI.YData)];
hold on
hF=fill(x,y,0.95*[1 1 1],'EdgeColor','none','FaceAlpha',0.5);
hLg=findobj(gcf,'Type','Legend');
hLg.String(end)=[];