function varargout=sigstar(groups,stats,nosort)
% SIGSTAR Add significance stars to bar charts, boxplots, line charts, etc,
%
% H = SIGSTAR(GROUPS,STATS,NSORT)
%
% Purpose
% Add stars and lines highlighting significant differences between pairs of groups.
% The user specifies the groups and associated p-values. The function handles much of
% the placement and drawing of the highlighting. Stars are drawn according to:
% * represents p<=0.05
% ** represents p<=1E-2
% *** represents p<=1E-3
%
%
% Inputs
% GROUPS - a cell array defining the pairs of groups to compare. Groups defined
% either as pairs of scalars indicating locations along the X axis or as
% strings corresponding to X-tick labels. Groups can be a mixture of both
% definition types.
% STATS - a vector of p-values the same length as GROUPS. If empty or missing it's
% assumed to be a vector of 0.05s the same length as GROUPS. Nans are treated
% as indicating non-significance.
% NSORT - optional, 0 by default. If 1, then significance markers are plotted in
% the order found in GROUPS. If 0, then they're sorted by the length of the
% bar.
%
% Outputs
% H - optionally return handles for significance highlights. Each row is a different
% highlight bar. The first column is the line. The second column is the text (stars).
%
%
% Examples
% 1.
% bar([5,2,1.5])
% sigstar({[1,2], [1,3]})
%
% 2.
% bar([5,2,1.5])
% sigstar({[2,3],[1,2], [1,3]},[nan,0.05,0.05])
%
% 3. **DOESN'T WORK IN 2014b**
% R=randn(30,2);
% R(:,1)=R(:,1)+3;
% boxplot(R)
% set(gca,'XTick',1:2,'XTickLabel',{'A','B'})
% H=sigstar({{'A','B'}},0.01);
% ylim([-3,6.5])
% set(H,'color','r')
%
% 4. Note the difference in the order with which we define the groups in the
% following two cases.
% x=[1,2,3,2,1];
% subplot(1,2,1)
% bar(x)
% sigstar({[1,2], [2,3], [4,5]})
% subplot(1,2,2)
% bar(x)
% sigstar({[2,3],[1,2], [4,5]})
%
% ALSO SEE: demo_sigstar
%
% KNOWN ISSUES:
% 1. Algorithm for identifying whether significance bar will overlap with
% existing plot elements may not work in some cases (see line 277)
% 2. Bars may not look good on exported graphics with small page sizes.
% Simply increasing the width and height of the graph with the
% PaperPosition property of the current figure should fix things.
%
% Rob Campbell - CSHL 2013
%Input argument error checking
%If the user entered just one group pair and forgot to wrap it in a cell array
%then we'll go easy on them and wrap it here rather then generate an error
if ~iscell(groups) & length(groups)==2
groups={groups};
end
if nargin<2
stats=repmat(0.05,1,length(groups));
end
if isempty(stats)
stats=repmat(0.05,1,length(groups));
end
if nargin<3
nosort=0;
end
%Check the inputs are of the right sort
if ~iscell(groups)
error('GROUPS must be a cell array')
end
if ~isvector(stats)
error('STATS must be a vector')
end
if length(stats)~=length(groups)
error('GROUPS and STATS must be the same length')
end
%Each member of the cell array GROUPS may be one of three things:
%1. A pair of indices.
%2. A pair of strings (in cell array) referring to X-Tick labels
%3. A cell array containing one index and one string
%
% For our function to run, we will need to convert all of these into pairs of
% indices. Here we loop through GROUPS and do this.
xlocs=nan(length(groups),2); %matrix that will store the indices
xtl=get(gca,'XTickLabel');
for ii=1:length(groups)
grp=groups{ii};
if isnumeric(grp)
xlocs(ii,:)=grp; %Just store the indices if they're the right format already
elseif iscell(grp) %Handle string pairs or string/index pairs
if isstr(grp{1})
a=strmatch(grp{1},xtl);
elseif isnumeric(grp{1})
a=grp{1};
end
if isstr(grp{2})
b=strmatch(grp{2},xtl);
elseif isnumeric(grp{2})
b=grp{2};
end
xlocs(ii,:)=[a,b];
end
%Ensure that the first column is always smaller number than the second
xlocs(ii,:)=sort(xlocs(ii,:));
end
%If there are any NaNs we have messed up.
if any(isnan(xlocs(:)))
error('Some groups were not found')
end
%Optionally sort sig bars from shortest to longest so we plot the shorter ones first
%in the loop below. Usually this will result in the neatest plot. If we waned to
%optimise the order the sig bars are plotted to produce the neatest plot, then this
%is where we'd do it. Not really worth the effort, though, as few plots are complicated
%enough to need this and the user can define the order very easily at the command line.
if ~nosort
[~,ind]=sort(xlocs(:,2)-xlocs(:,1),'ascend');
xlocs=xlocs(ind,:);groups=groups(ind);
stats=stats(ind);
end
%-----------------------------------------------------
%Add the sig bar lines and asterisks
holdstate=ishold;
hold on
H=ones(length(groups),2); %The handles will be stored here
y=ylim;
yd=myRange(y)*0.05; %separate sig bars vertically by 5%
for ii=1:length(groups)
thisY=findMinY(xlocs(ii,:))+yd;
H(ii,:)=makeBar(xlocs(ii,:),thisY,stats(ii));
end
%-----------------------------------------------------
%Now we can add the little downward ticks on the ends of each line. We are
%being extra cautious and leaving this it to the end just in case the y limits
%of the graph have changed as we add the highlights. The ticks are set as a
%proportion of the y axis range and we want them all to be the same the same
%for all bars.
yd=myRange(ylim)*0.01; %Ticks are 1% of the y axis range
for ii=1:length(groups)
y=get(H(ii,1),'YData');
y(1)=y(1)-yd;
y(4)=y(4)-yd;
set(H(ii,1),'YData',y)
end
%Be neat and return hold state to whatever it was before we started
if ~holdstate
hold off
elseif holdstate
hold on
end
%Optionally return the handles to the plotted significance bars (first column of H)
%and asterisks (second column of H).
if nargout>0
varargout{1}=H;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Internal functions
function H=makeBar(x,y,p)
%makeBar produces the bar and defines how many asterisks we get for a
%given p-value
if p<=1E-3
stars='***';
elseif p<=1E-2
stars='**';
elseif p<=0.05
stars='*';
elseif isnan(p)
stars='n.s.';
else
stars='';
end
x=repmat(x,2,1);
y=repmat(y,4,1);
H(1)=plot(x(:),y,'-k','LineWidth',1.5);
%Increase offset between line and text if we will print "n.s."
%instead of a star.
if ~isnan(p)
offset=0.005;
else
offset=0.02;
end
H(2)=text(mean(x(:)),mean(y)+myRange(ylim)*offset,stars,...
'HorizontalAlignment','Center',...
'BackGroundColor','none');
function Y=findMinY(x)
%Find the minimum y value needed to clear all the plotted data present
%over a given range of x values. This allows the significance bar to
%be plotted in the best location.
%
%First look for patch objects (bars in a bar-chart, most likely)
if verLessThan('matlab','8.4.0')
p=findobj(gca,'Type','Patch');
xd=get(p,'XData');
else
p=findobj(gca,'Type','bar')
xd=p.XData;
end
if iscell(xd) & verLessThan('matlab','8.4.0')
xd=groupedBarFix(xd,'x');
end
xd(xd<x(1))=0;
xd(xd>x(2))=0;
overlapping=any(xd,1); %These x locations overlap
%Find the corresponding y values
if verLessThan('matlab','8.4.0')
yd=get(p,'YData');
else
yd=p.YData;
end
if iscell(yd) & verLessThan('matlab','8.4.0')
yd=groupedBarFix(yd,'y');
end
yd=yd(:,overlapping);
%So we must have a value of at least Y in order to miss all the
%plotted bar data:
Y=max(yd(:));
%Now let's check if any other plot elements (such as significance bars we've
%already added) exist over this range of x values.
%
% NOTE! This code doesn't identify all cases where there is overlap.
%For example, if you have a significance bar going from 3 to 7 along the x
%axis and you then try to add a new one from 4 to 5, then it won't see the
%existing one as overlapping. However, I've yet to find this a problem in
%practice so I'll leave things be. Can easily be fixed if leads to bugs.
p=findobj(gca,'Type','Line');
tmpY=nan(1,length(p));
for ii=1:length(p)
xd=get(p(ii),'XData');
xd(xd<x(1))=0;
xd(xd>x(2))=0;
overlapping=xd>0; %These x locations overlap
if ~any(overlapping)
continue
end
clear xd
yd=get(p(ii),'YData');
yd=yd(overlapping);
tmpY(ii)=max(yd);
end
Y=max([Y,tmpY]);
%The patch coords of grouped error bars aren't a matrix but a cell array. This needs to be
%converted to a suitable matrix in order for the code to work. This function does that.
function out=groupedBarFix(in,xy)
out=ones([size(in{1}),length(in)]);
for ii=1:length(in)
out(:,:,ii)=in{ii};
end
switch xy
case 'x'
out=mean(out,3);
case 'y'
out=max(out,[],3);
end
%replacement for stats toolbox range function
function rng=myRange(x)
rng = max(x) - min(x);