%Hesterberg - Example 18.12
%Permutation Test of Treatment vs. Control Effects on Reading Performance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read reading performance scores %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
fid=fopen([pwd,'\data\ta18_004.txt']);
fgetl(fid);  %skip header line
AoA=textscan(fid,'%d%s'); %read Numeric and Text columns from data file
reading=AoA{1};   %extract reading score
group=AoA{2};     %extract group label (Treatment/Control)
fclose(fid);

%segregate treatment and control groups
treatment = double(reading(1:21));
control = double(reading(22:44));

%seed the random number generator
rand('state',sum(100*clock));

% traditional 1-tailed t-test
[h,p,ci,stats]=ttest2(treatment,control,0.01,'right');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Resampling to generate Permutation Distribution %
% of Difference between Means                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
samplesize1=21;
samplesize2=23;
reading=randomize_matrix(reading); %shuffle pseudo-population first
n_resamples=1000;
for i=1:n_resamples
    
    clear s1; clear s2;
    
    %generate permutation samples
    [s1,s2]=randperm2(reading,samplesize1);
    %compute/save difference between permutation sample means
    remdiff(i)=mean(s1)-mean(s2);
end
%
%compute and display summary statistics
remean=mean(remdiff);  %bootstrap distribution mean
restderr=std(remdiff); %bootstrap distribution standard error
remdiff=remdiff';
%percentile cutoffs
p1=prctile(remdiff,1);
p25=prctile(remdiff,2.5);
p5=prctile(remdiff,5);
p95=prctile(remdiff,95);
p975=prctile(remdiff,97.5);
p99=prctile(remdiff,99);
display('Permutation Distribution Summary statistics:');
display(['Mean:        ',num2str(remean)]);
display(['Std. Error:  ',num2str(restderr)]);
display(['1st  %ile:   ',num2str(p1)]);
display(['2.5  %ile:   ',num2str(p25)]);
display(['5th  %ile:   ',num2str(p5)]);
display(['95th %ile:   ',num2str(p95)]);
display(['97.5 %ile:   ',num2str(p975)]);
display(['99th %ile:   ',num2str(p99)]);
display('');
delta=mean(treatment)-mean(control);
display(['Difference between Observed Means = ',num2str(delta)]);
display(['t-test:  ',num2str(stats.tstat)])
display(['p-level: ',num2str(p)])
display(['df:      ',num2str(stats.df)])
display(['sd:      ',num2str(stats.sd)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot bootstrap distribution with parametric overlay %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear n; clear hours;
hours=[min(remdiff):1:max(remdiff)]; %data bins
n=histc(remdiff,hours);
figure(1); clf;
subplot(2,1,1);
set(gcf,'PaperPositionMode','auto');
set(gcf,'Position',[650 100 600 800]);
bar(hours,n,'k');
hold on;
axis([-20 20 0 max(n)+2]);
xlabel('delta Permutation Means');
ylabel('Frequency');
title_str(1)={['Permutation Distribution of Difference between Means']};
title_str(2)={'Hesterberg Example 18.12'};
title_str(3)={['(Number of Resamples = ',num2str(n_resamples),')']};
title(title_str);
%annotate randomization distribution
f11str(1)={['Mean       : ',num2str(remean)]};
f11str(2)={['Std. Error:  ',num2str(restderr)]};
f11str(3)={['1st  %ile:   ',num2str(p1)]};
f11str(4)={['5th  %ile:   ',num2str(p5)]};
f11str(5)={['95th %ile:   ',num2str(p95)]};
f11str(6)={['99th %ile:   ',num2str(p99)]};
text(-18,max(n)/2,f11str);
%generate normal-quantile plot to assess normality
subplot(2,1,2);
qqplot(remdiff);
hold on;
ylabel('delta Permutation Means');
title('');
hold off;
