function pairedPlot(data,varname,xlabels,titles)

% plot_prepost(data,varname,xlabels,titles)
%
% INPUTS
% data     : Nx2 matrix
% varname  : variable name string
% xlabels  : {'Before','After'} optional
% titles   : {'Left Panel','Right Panel'} optional
%
% DEFAULTS
% xlabels = {'1','2'}
% titles  = {'',''}
%
% Sample call
% 

%% ---------- Defaults ----------
if nargin < 2 || isempty(varname)
    varname = '';
end

if nargin < 3 || isempty(xlabels)
    xlabels = {'1','2'};
end

if nargin < 4 || isempty(titles)
    titles = {'',''};
end

%% ---------- Data cleaning ----------
if size(data,2) ~= 2
    error('Data must be Nx2 matrix')
end

data = data(all(~isnan(data),2),:);

before = data(:,1);
after  = data(:,2);

n = length(before);

%% ---------- Means + 95% CI ----------
means = [mean(before) mean(after)];

sem_before = std(before)/sqrt(n);
sem_after  = std(after)/sqrt(n);

tcrit = tinv(0.975,n-1);

cis = [tcrit*sem_before tcrit*sem_after];

%% ---------- Stats ----------
[~,p,~,stats] = ttest(before,after);

tval = stats.tstat;
df   = stats.df;

diff_scores = before-after;

cohens_d = mean(diff_scores)/std(diff_scores);

%% Percent change
pct_change = (mean(after)-mean(before))/mean(before)*100;

%% Bayes Factor
BF10 = jzs_bf_ttest(tval,n);

if BF10 > 100
    bftext = 'BF_{10} > 100';
else
    bftext = sprintf('BF_{10} = %.2f',BF10);
end

%% Format p
if p < 0.001
    ptxt = 'p < 0.001';
else
    ptxt = sprintf('p = %.3f',p);
end

stat_text = sprintf(['t(%d) = %.2f\n' ...
                     '%s\n' ...
                     'Cohen''s d = %.2f\n' ...
                     '%s\n' ...
                     'Percent Change = %.1f%%'],...
                     df,tval,ptxt,cohens_d,bftext,pct_change);

%% ---------- Axis limits ----------
ymin = min([before;after]);
ymax = max([before;after]);

if ymin >= 0 && ymax <= 100
    yl = [0 100];
else
    pad = (ymax-ymin)*0.15;
    yl = [ymin-pad ymax+pad];
end

%% ---------- Figure ----------
figure
set(gcf,'color','w','position',[100 100 1150 500])

font_size = 24;
stats_font = 18;

blue  = [0.12 0.35 0.95];
green = [0.05 0.65 0.25];

%% =====================================================
%% LEFT PANEL
%% =====================================================

subplot(1,2,1)
hold on

b = bar(1:2,means,0.65,'FaceColor','flat','LineWidth',2.8);
b.CData(1,:) = blue;
b.CData(2,:) = green;

errorbar(1:2,means,cis,'k','LineStyle','none','LineWidth',2.8,'CapSize',14)

xlim([0.4 2.6])
ylim(yl)

set(gca,'XTick',[1 2],...
        'XTickLabel',xlabels,...
        'FontSize',font_size,...
        'LineWidth',2.5)

ylabel('Score (mean ± 95% CI)','FontSize',font_size)

title(titles{1},'FontSize',font_size,'FontWeight','bold')

box off

%% Stats block
text(0.55,yl(2)-0.09*(yl(2)-yl(1)),stat_text,...
    'FontSize',stats_font,...
    'VerticalAlignment','top')

%% Significance bracket
y_bracket = max(means+cis)+0.07*(yl(2)-yl(1));
h = 0.02*(yl(2)-yl(1));

plot([1 1 2 2],[y_bracket-h y_bracket y_bracket y_bracket-h],...
    'k','LineWidth',2)

if p < 0.001
    sig='***';
elseif p < 0.01
    sig='**';
elseif p < 0.05
    sig='*';
else
    sig='n.s.';
end

text(1.5,y_bracket+0.012*(yl(2)-yl(1)),sig,...
    'HorizontalAlignment','center',...
    'FontSize',stats_font,...
    'FontWeight','bold')

%% =====================================================
%% RIGHT PANEL
%% =====================================================

subplot(1,2,2)
hold on

xvals = [ones(n,1); 2*ones(n,1)];
yvals = [before; after];

colors = [repmat(blue,n,1); repmat(green,n,1)];

s = swarmchart(xvals,yvals,34,colors,'filled');

s.XJitter = 'density';
s.XJitterWidth = 0.24;
s.MarkerFaceAlpha = 0.9;
s.MarkerEdgeAlpha = 0.9;

plot(1,means(1),'ko','MarkerFaceColor','k','MarkerSize',9)
plot(2,means(2),'ko','MarkerFaceColor','k','MarkerSize',9)

errorbar(1:2,means,cis,'k','LineStyle','none','LineWidth',2.8,'CapSize',14)

xlim([0.4 2.6])
ylim(yl)

set(gca,'XTick',[1 2],...
        'XTickLabel',xlabels,...
        'FontSize',font_size,...
        'LineWidth',2.5)

ylabel('Score (mean ± 95% CI)','FontSize',font_size)

title(titles{2},'FontSize',font_size,'FontWeight','bold')

box off

end


function BF10 = jzs_bf_ttest(t,n,r)

if nargin < 3
    r = 0.707;
end

v = n - 1;

integrand = @(g) ...
    (1 + n.*g.*r.^2).^(-1/2) .* ...
    (1 + t.^2 ./ ((1 + n.*g.*r.^2).*v)).^(-(v+1)/2) .* ...
    (2*pi).^(-1/2) .* g.^(-3/2) .* exp(-1./(2.*g));

num = integral(integrand,0,Inf,'ArrayValued',true);

den = (1 + t.^2./v).^(-(v+1)/2);

BF10 = num ./ den;

end