--- a
+++ b/LV_evaluation_code/functions/BlandAltman.m
@@ -0,0 +1,432 @@
+% BlandAltman - draws a Blant-Altman and correlation graph for two
+% datasets.
+%
+% BlandAltman(data1, data2)
+% BlandAltman(data1, data2,label) - Names of data sets. Formats can be
+%   - {'Name1'}
+%   - {'Name1, 'Name2'}
+%   - {'Name1, 'Name2', 'Units'}
+% BlandAltman(data1, data2,label,tit,gnames)
+% BlandAltman(data1, data2,label,tit,gnames,corrinfo) - specifies what
+% information to display on the correlation chart as a cell of string in 
+% order of top to bottom. The following codes are available:
+%  - 'eq'     slope and intercept equation
+%  - 'int%'   intercept as % of mean values
+%  - 'r'      pearson r-value
+%  - 'r2'     pearson r-value squared
+%  - 'rho'    Spearman rho value
+%  - 'SSE'    sum of squared error
+%  - 'n'      number of data points used
+%  * if not specified or empty, default is: {'eq';'r2';'SSE';'n'}
+% BlandAltman(data1, data2,label,tit,gnames,corrinfo,BAinfo) - specifies what
+% information to display on the Bland-Altman chart as for corrinfo, but
+% with the following codes:
+%  - 'RPC'    reproducibility coefficient (1.96*SD)
+%  - 'RPC(%)' reproducibility coefficient and % of mean values
+%  - 'CV'     coefficient of variation (SD of mean values in %)
+%  * if not specified or empty, default is: {'RPC(%)';'CV'}
+% BlandAltman(data1, data2,label,tit,gnames,corrinfo,BAinfo,limits) -
+% specifies the axes limits:
+%  - scalar - lower limit (eg. 0)
+%  - [min max] - specifies minimum and maximum
+%  - 'tight' - minimum and maximum of data.
+%  - 'auto' - plot default. {default option}
+% BlandAltman(data1, data2,label,tit,gnames,corrinfo,BAinfo,limits,colors) -
+% specify the order of group colors. (eg. 'brg' for blue, red, green) or
+% RGB columns.
+% BlandAltman(data1, data2,label,tit,gnames,corrinfo,BAinfo,limits,colors,symbols)
+% - specify the order of symbols. (eg. 'sod' for squares, circles, dots).
+% Alternatively can be set to 'Num' to display the subject number.
+% BlandAltman(fig, ...) - specify a figure handle in which to
+% display
+% figure in which the Bland-Altman and correlation will be displyed 
+% BlandAltman(ah, ...) - specify an axes which will be replaced by the 
+% Bland-Altman and correlation exes.
+% cr = BlandAltman(...) - return the coefficient of reproducibility
+% (1.96*sd)
+% [cr fig] = BlandAltman(...) - also return the figure handles
+% [cr fig sstruct] = BlandAltman(...) - also return the structure of
+% statistics for the analysis
+
+% by Ran Klein 2010
+% 2013-02-20  RK  added unit lableling in SSE and RPC labels.
+% 2013-02-20  RK  switched to equal scaling on BA yaxis.
+% 2014-01-15  RK  added colors and symbols input arguments.
+
+function [cr, fig, sstruct] = BlandAltman(varargin)
+
+% Added feature to control y-axis scaling on BA figure. Follow code below
+% for details. May become a parameter with demand.
+fixylim = true;
+
+if isscalar(varargin{1}) && isequal(size(varargin{1}),[1 1]) && ishandle(varargin{1})
+	shift = 1;
+	fig = varargin{1};
+else
+	shift = 0;
+	fig = [];
+end
+data1 = varargin{shift+1};
+data2 = varargin{shift+2};
+if nargin>=shift+3
+	label = varargin{shift+3};
+else
+	label = '';
+end
+if nargin>=shift+4
+	tit = varargin{shift+4};
+else
+	tit = '';
+end
+if nargin>=shift+5
+	gnames = varargin{shift+5};
+else
+	gnames = '';
+end
+if nargin>=shift+6 && ~isempty(varargin{shift+6})
+	if ischar(varargin{shift+6})
+		corrinfo = varargin(shift+6);
+	else
+		corrinfo = varargin{shift+6};
+	end
+else
+	corrinfo = {'eq';'r';'SSE';'n'};
+end
+if nargin>=shift+7 && ~isempty(varargin{shift+7})
+	if ischar(varargin{shift+7})
+		BAinfo = varargin(shift+7);
+	else
+		BAinfo = varargin{shift+7};
+	end
+else
+	BAinfo = {'RPC(%)';'CV'};
+end
+
+if nargin>=shift+8
+	axesLimits = varargin{shift+8};
+else
+	axesLimits = 'auto';
+end
+	
+if nargin>=shift+9 && ~isempty(varargin{shift+9})
+	colors = varargin{shift+9};
+else
+	colors = 'rbgmcky';
+end
+if ~ischar(colors)
+	if size(colors,2)~=3
+		if size(colors,1)==3
+			colors = colors';
+		else
+			error('Colors must be specified in either character codes or RGB');
+		end
+	end
+elseif size(colors,1)==1
+	colors = colors';
+end
+
+if nargin>=shift+10 && ~isempty(varargin{shift+10})
+	symb = varargin{shift+10};
+else
+	symb = 'sodp^v';
+end
+
+markersize = 4;
+xdatamode = 1; % use the (1) mean of data1 and data2 or (2) data1 for x position on bland altman
+
+units = '';
+if iscell(label)
+	if length(label)==1
+		labelx = label{1};
+		labely = label{1};
+		labelm = label{1};
+		labeld = ['\Delta ' label{1}];
+	elseif length(label)==2
+		labelx = label{1};
+		labely = label{2};
+		labelm = ['Mean ' label{1} ' & ' label{2}];
+		labeld = [label{2} ' - ' label{1}];
+	else % units also provided
+		units = label{3};
+		labelx = [label{1} ' (' units ')'];
+		labely = [label{2} ' (' units ')'];
+		labelm = ['Mean ' label{1} ' & ' label{2} ' (' units ')'];
+		labeld = [label{2} ' - ' label{1} ' (' units ')'];
+	end	
+else
+	labelx = label;
+	labely = label;
+	labelm = label;
+	labeld = ['\Delta ' label];
+end
+if isempty(units)
+	unitsstr = '';
+else
+	unitsstr = [' ' units];
+end
+
+s = size(data1);
+if ~isequal(s,size(data2));
+	error('data1 and data2 must have the same size');
+end
+
+switch length(s)
+	case 1
+		s = [s 1 1];
+	case 2
+		s = [s 1];
+	case 3
+	otherwise
+		error('Data have too many dimension');
+end
+n = s(1); % number of elements in each group
+groups = numel(data1)/n;
+if size(colors,1)<s(3)
+	error('More groups than colors specified. Use the colors input variable to specify colors for each group.');
+end
+if ~strcmpi(symb,'Num') && length(symb)<s(2)
+	error('More subgroups than symbolss specified. Use the symbols input variable to specify symbols for each subgroup, or use the ''Num'' option.');
+end
+
+data1 = reshape(data1, [numel(data1),1]);
+data2 = reshape(data2, [numel(data2),1]);
+mask = isfinite(data1) & isnumeric(data1) & isfinite(data2) & isnumeric(data2);
+
+if isempty(fig)
+	fig = figure;
+	set(fig,'units','centimeters','position',[3 3 20 10],'color','w');
+	cah = subplot(121);
+	dah = subplot(122);
+elseif strcmpi(get(fig,'type'),'figure')
+	cah = subplot(121);
+	dah = subplot(122);
+elseif strcmpi(get(fig,'type'),'axes')
+	ah = fig;
+	pos = get(ah,'position');
+	fig = get(ah,'parent');
+	delete(ah);
+	cah = axes('parent',fig,'position',[pos(1) pos(2) pos(3)/2 pos(4)]);
+	dah = axes('parent',fig,'position',[pos(1)+pos(3)/2 pos(2) pos(3)/2 pos(4)]);
+else
+	error('What in tarnations is the handle that was passed to Bland-Altman????')
+end
+set(cah,'tag','Correlation Plot');
+set(dah,'tag','Bland Altman Plot');
+
+%% Correlation
+hold(cah,'on');
+for groupi=1:groups
+	if strcmpi(symb,'Num')
+		for i=1:n
+			text(data1((groupi-1)*n+i),data2((groupi-1)*n+i),num2str(i),'parent',cah,'fontsize',markersize,'color',colors(floor((groupi-1)/s(2))+1,:), 'HorizontalAlignment','Center', 'VerticalAlignment','Middle');
+		end
+	else
+		if s(3)==1
+			marker = symb(1);
+			color = colors(groupi,:);
+		else
+			marker = symb(rem(groupi-1,s(2))+1);
+			color = colors(floor((groupi-1)/s(2))+1,:);
+		end
+		ph=plot(cah,data1((groupi-1)*n+(1:n)),data2((groupi-1)*n+(1:n)),marker,'color',color);
+		set(ph,'markersize',markersize);
+	end
+end
+% Linear regression
+[polyCoefs, S] = polyfit(data1(mask),data2(mask),1);
+r = corrcoef(data1(mask),data2(mask)); r=r(1,2);
+rho = corr(data1(mask),data2(mask),'type','Spearman');
+N = sum(mask);
+SSE = sqrt(sum((polyval(polyCoefs,data1(mask))-data2(mask)).^2)/(N-2));
+
+if ischar(axesLimits)
+	if strcmpi(axesLimits,'Auto')
+		% Workaround - Add invisible minimum and maximum point to fix Auto axes limits (text
+		% does not count for axis('auto')
+		if strcmpi(symb,'Num')
+			mindata = min( min(data1(mask)), min(data2(mask)) );
+			maxdata = max( max(data1(mask)), max(data2(mask)) );
+			ph = plot(cah, [mindata maxdata], [mindata maxdata], '.', 'Visible','on');
+		end
+		axesLimits = axis(cah); 
+		axesLimits(1) = min(axesLimits(1),axesLimits(3));
+		axesLimits(2) = max(axesLimits(2),axesLimits(4));
+		if strcmpi(symb,'Num')
+			delete(ph);
+		end
+	elseif strcmpi(axesLimits,'Tight')
+		axesLimits(1) = min( min(data1(mask)), min(data2(mask)) );
+		axesLimits(2) = max( max(data1(mask)), max(data2(mask)) );
+	else
+		error(['Unknown axes limit option (' axesLimits ') detected.']);
+	end
+else
+	if length(axesLimits)==1
+		a = axis(cah);
+		axesLimits(2) = max(a(2),a(4));
+	else
+		% Do nothing
+	end
+end
+axesLimits(3) = axesLimits(1);
+axesLimits(4) = axesLimits(2);
+
+axis(cah,axesLimits); axis(cah,'square');
+plot(cah,axesLimits(1:2), polyval(polyCoefs,axesLimits(1:2)),'-k');
+h = plot(cah,axesLimits(1:2),axesLimits(1:2),':'); set(h,'color',[0.6 0.6 0.6]);
+if 0 % Add 95% CI lines
+	xfit = axesLimits(1):(axesLimits(2)-axesLimits(1))/100:axesLimits(2);
+	[yfit, delta] = polyconf(polyCoefs,xfit,S);
+	h = [plot(cah,xfit,yfit+delta);...
+		plot(cah,xfit,yfit-delta)];
+	set(h,'color',[0.6 0.6 0.6],'linestyle','-');
+end
+corrtext = {};
+for i=1:length(corrinfo)
+	switch lower(corrinfo{i})
+		case 'eq'
+			if polyCoefs(2)>0
+				corrtext = [corrtext; ['y=' num2str(polyCoefs(1),3) 'x+' num2str(polyCoefs(2),3)]];
+			else
+				corrtext = [corrtext; ['y=' num2str(polyCoefs(1),3) 'x' num2str(polyCoefs(2),3)]];
+			end
+		case 'int%', corrtext = [corrtext; ['intercept=' num2str(polyCoefs(2)/mean(data1+data2)*2*100,3) '%']];
+		case 'r2', corrtext = [corrtext; ['r^2=' num2str(r^2,4)]];
+		case 'r', corrtext = [corrtext; ['r=' num2str(r,4)]];
+		case 'rho', corrtext = [corrtext; ['rho=' num2str(rho,4)]];
+		case 'sse', corrtext = [corrtext; ['SSE=' num2str(SSE,2) unitsstr]];
+		case 'n', corrtext = [corrtext; ['n=' num2str(N)]];
+	end
+end
+text(axesLimits(1)+0.01*(axesLimits(2)-axesLimits(1)),axesLimits(1)+0.9*(axesLimits(2)-axesLimits(1)),corrtext,'parent',cah);
+xlabel(cah,labelx); ylabel(cah,labely);
+
+%% Differences
+set(dah,'units','normalized');
+hold(dah,'on');
+for groupi=1:groups
+	d1 = data1((groupi-1)*n+(1:n));
+	d2 = data2((groupi-1)*n+(1:n));
+	dif = d2-d1;
+	if strcmpi(symb,'Num')
+		for i=1:n
+			if xdatamode==1
+				text(mean([d1(i),d2(i)]), dif(i), num2str(i), 'parent',dah,'fontsize',markersize,'color',colors(floor((groupi-1)/s(2))+1,:));
+			else
+				text(d1(i), dif(i), num2str(i), 'parent',dah,'fontsize',markersize,'color',colors(floor((groupi-1)/s(2))+1,:));
+			end
+		end
+	else
+		if s(3)==1
+			marker = symb(1);
+			color = colors(groupi,:);
+		else
+			marker = symb(rem(groupi-1,s(2))+1);
+			color = colors(floor((groupi-1)/s(2))+1,:);
+		end
+		
+		if xdatamode==1
+			ph = plot(dah,mean([d1,d2],2),dif,marker,'color',color);
+		else
+			ph = plot(dah,d1,dif,marker,'color',color);
+		end
+		set(ph,'markersize',markersize);
+	end
+end
+axis(dah,'square')
+xlabel(dah,labelm); ylabel(dah,labeld);
+
+% add std-dev lines
+st = std(data2(mask)-data1(mask));
+mn = mean(data2(mask)-data1(mask));
+[h, p] = ttest(data2(mask)-data1(mask),0);
+cr = 1.96*st;
+
+% fix limits to +/- data limit
+if fixylim 
+	a = [axesLimits(1:2) [-1 1]*abs(axesLimits(2)-axesLimits(1))/2];
+	axis(dah, a);
+end
+
+plot(a(1:2),mn+[0 0],'k')
+% plot(a(1:2),mn+st*[1 1],'k')
+% plot(a(1:2),mn-st*[1 1],'k')
+plot(a(1:2),mn+cr*[1 1],':k')
+plot(a(1:2),mn-cr*[1 1],':k')
+a = axis(dah);
+if fixylim
+	fontsize = 6;
+	text(a(2),mn+cr, [num2str(mn+cr,2) ' (+1.96SD)'],'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+	text(a(2),mn,[num2str(mn,2) ' [p=' num2str(p,2) ']'],'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+	text(a(2),mn-cr, [num2str(mn-cr,2) ' (-1.96SD)'],'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+else
+	fontsize = 8;
+	text(a(2),mn+cr,{'+1.96SD',num2str(mn+cr,2)},'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+	text(a(2),mn,{num2str(mn,2),['p=' num2str(p,2)]},'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+	text(a(2),mn-cr,{num2str(mn-cr,2),'-1.96SD'},'HorizontalAlignment','left','VerticalAlignment','middle','fontsize',fontsize);
+end
+BAtext = {};
+
+for i=1:length(BAinfo)
+	switch lower(BAinfo{i})
+		case 'rpc', BAtext = [BAtext; ['{\bfRPC: ' num2str(cr,2) unitsstr '}']];
+		case 'rpc(%)', BAtext = [BAtext; ['{\bfRPC: ' num2str(cr,2) unitsstr '} (' num2str(100*cr/mean((data2(mask)+data1(mask))/2),2) '%)']];
+		case 'cv', BAtext = [BAtext; ['CV: ' num2str(100*st/mean((data1(mask)+data2(mask))/2),2) '%']];
+		case 'p', BAtext = [BAtext; ['p-value: ' num2str(p,4) ]]; 
+		case 'ks' % Kolmogorov-Smirnov test that difference-data is Gaussian
+			ddata = data1(:)-data2(:);
+			[h, p] = kstest((ddata-mean(ddata))/std(ddata));
+			BAtext = [BAtext; ['KS p-value: ' num2str(p)]];
+		case 'kurtosis' % Kolmogorov-Smirnov test that difference-data is Gaussian
+			ddata = data1(:)-data2(:);
+			BAtext = [BAtext; ['kurtosis: ' num2str(kurtosis(ddata))]];
+	end
+end
+text(a(2),a(4),BAtext,'interpreter','tex','HorizontalAlignment','right','VerticalAlignment','top');
+
+if ~isempty(tit)
+	h = suptitle(tit);
+	set(h,'interpreter','tex');
+end
+
+% Add legend
+if ~strcmpi(symb,'Num') && ~isempty(gnames)
+	lh = legend('show');
+	if iscell(gnames)
+		if length(gnames)==2 
+			if iscell(gnames{1}) 
+				temp = cell(1,groups);
+				for groupi=1:length(gnames{1})
+					for j=1:length(gnames{2})
+						temp{groupi+(j-1)*length(gnames{1})} = [gnames{1}{groupi} '-' gnames{2}{j}];
+					end
+				end	
+				gnames = temp;
+			elseif iscell(gnames{2})
+				gnames = strcat(gnames{1}, '-', gnames{2});
+			end
+		end
+	end
+	cpos = get(cah,'Position');
+	dpos = get(dah,'Position');
+	set(cah,'Position',cpos+[0 0.07 0 0]);
+	set(dah,'Position',dpos+[0 0.07 0 0]);
+	set(lh,'string',gnames,'orientation','horizontal');
+	drawnow;
+	set(lh,'units','normalized');
+	pos = get(lh,'position'); pos = min(pos(3),0.9);
+	set(lh,'position',[(1-pos)/2 0.02 pos 0.05]);
+else
+% 	set(lh,'visible','off');
+end
+
+if nargout>2
+	sstruct = struct('N',N,...
+		'CR', cr,...
+		'r',r,...
+		'r2',r^2,...
+		'SSE',SSE,...
+		'rho',rho,...
+		'Slope',polyCoefs(1),...
+		'Intercept',polyCoefs(2));
+end
\ No newline at end of file