
disp('Multiple Comparison procedures for One-Way layouts');
disp('All-pairs comparisons with multiple-comparison correction');
disp(' ');
disp('Numbers of glucocorticoid receptor sites per Leukocyte cell');
disp('  for controls and 4 types of leukemia/anemia');
disp('Hollander & Wolfe (Table 6.4, p201)');

%% Data in 5 rows
xx1 = [ 3500 3500 3500 4000 4000 4000 4300 4500 4500 4900 5200 6000 6750 8000];
xx2 = [ 5710  6110  8060  8080 11400 ];
xx3 = [ 2930  3330  3580  3880  4280  5120 ];
xx4 = [ 6320  6860 11400 14000 ];
xx5 = [ 3230  3880  7640  7890  8280 16200 18250 29900 ];

n1=length(xx1);  n2=length(xx2);  n3=length(xx3);
   n4=length(xx4);  n5=length(xx5);
num = [ n1  n2  n3  n4  n5 ];
kk = length(num);  % Number of treatment groups
nsum = sum(num);   % Total number of observations
nmax = max(num);   % Maximum number of observations per group

disp(' ');
fprintf('Grp #1 (n=%2d) Normal subjects\n',n1);
fprintf('Grp #2 (n=%2d) Hairy-cell anemia\n',n2);
fprintf('Grp #3 (n=%2d) Chronic lymphatic leukemia\n',n3);
fprintf('Grp #4 (n=%2d) Chronic myelocytic leukemia\n',n4);
fprintf('Grp #5 (n=%2d) Acute leukemia\n',n5);

%% Put data in an nmax by kk array for display
xxa=zeros(nmax,kk);
for i=1:n1;  xxa(i,1)=xx1(i);  end;
for i=1:n2;  xxa(i,2)=xx2(i);  end;
for i=1:n3;  xxa(i,3)=xx3(i);  end;
for i=1:n4;  xxa(i,4)=xx4(i);  end;
for i=1:n5;  xxa(i,5)=xx5(i);  end;
disp('  ');
disp('Data:');
fprintf('    Grp#1    Grp#2    Grp#3    Grp#4    Grp#5\n');
for i=1:nmax;
  for j=1:kk;
     if i <= num(j)  fprintf('    %5d', xxa(i,j));
     else fprintf('         ');  end;  end;
  fprintf('\n');  end;

disp(' ');
fprintf('%d treatment groups with %d total observations:\n',kk,nsum);

% Put data in an nsum by 2 array
zz = [ xx1'    ones(n1,1);
       xx2'  2*ones(n2,1);
       xx3'  3*ones(n3,1);
       xx4'  4*ones(n4,1);
       xx5'  5*ones(n5,1) ];  %'

% Add midranks for all treatment groups together
[midr tiegr] = getmidranks(zz);
zz = [ [1:nsum]'  zz  midr  tiegr];  %'

disp('Ord, Values, Group, Midranks, Tiegroups:');
fprintf(' %3d.  %5d   %d    %4.1f   %2g\n', zz');  %'

disp(' ');
disp('Sample means and rank means for treatment groups:');

num2=zeros(1,kk);  DataSum=zeros(1,kk);  RankSum=zeros(1,kk);
vals=zz(:,2);  group=zz(:,3);  ranks=zz(:,4);

% Find numbers, data sum, and rank sum for each treatment group
for g=1:kk;
    num2(g)=sum(group==g);
    DataSum(g)=sum(vals.*(group==g));
    RankSum(g)=sum(ranks.*(group==g));  end;
% DataSum, RankSum, and num2 are all 1xkk row vectors
DataAv=DataSum./num2;
RankAv=RankSum./num2;

ddsum = [ [1:kk]'  num2'  DataAv'  RankAv' ];

fprintf('  #%d (n=%2d):  DataAv=%9.2f   RankAv=%6.2f\n',ddsum');  %'


% Compute H statistic
Ravc = RankAv - (nsum+1)/2;
hh = (12/(nsum*(nsum+1)))*sum(num.*Ravc.*Ravc);

fprintf('\nKruskal-Wallis H = %g (with no tie correction)\n', hh);

% Observed Hscore:
Hscore=sum(RankSum.*RankSum./num);
fprintf ('  For permutations:  Hscore = %g\n', Hscore);

% Adjust for tie correction
Tiesum=0;
for i=1:nsum;  tt=tiegr(i);
  if tt>0  ntg=i;  Tiesum = Tiesum + (tt-1)*tt*(tt+1);
  else break;  end;  end;

Tiecorr = Tiesum/((nsum-1)*nsum*(nsum+1));
fprintf('  Tiesum=%d    TieCorr=%g    (nsum=%d)\n', Tiesum, Tiecorr, nsum);

hhp = hh/(1-Tiecorr);

fprintf('Kruskal-Wallis HP = %g (with tie correction)\n', hhp);
disp('Large-sample chi-square approximation for H:');
df=kk-1;
Pval = pvchisq(hhp,df);
%Pval = pvtdist(df,hhp);
fprintf('    P= %7.5f  (df=%d)\n', Pval,df);

% Carry out permutation test for large-sample Kruskal-Wallis approximation:
% Note that H = A_n * HScore + C_n   for  A_N=12/(nsum*(nsum+1)), C_N=-3(nsum+1)
%   for  Hscore = Sum RankSum^2/num
% Thus simulating  Hscore  is the same as simulating  H (hh)

fprintf('\nPermutation test for large-sample Kruskal-Wallis approximation:\n');

disp(' ');
disp('Initializing the random-number generator at');
Startseed=12345678;
disp(Startseed);
rand('state',Startseed);

disp('Number of permutations:');
nsims=10000;
disp(nsims);

% Do nsims permutations for a test of Hscore
Nge=0;  group=zz(:,3);  midr=zz(:,4);
for ns=1:nsims;
  %% Permute copy of midranks:
  midrperm=midr(randperm(nsum));
  %% Store rank sums for permuted midranks in rrsum
  rrsum=zeros(1,kk);
  for i=1:nsum;  rrsum(group(i))=rrsum(group(i))+midrperm(i);  end;
%  While the following code works, it is slower
%  for g=1:kk;  rrsum(g)=sum(midrperm.*(group==g));  end;
  %% Calculate the hscore for the permuted midranks
  Hrand = sum(rrsum.*rrsum./num);
  if Hrand >= Hscore  Nge=Nge+1;  end;
  end;

disp('Number of simulations with values >= Hscore and total number:');
fprintf('    Nge=%d    nsims=%d\n', Nge, nsims);

phat=Nge/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% CI for true P-value bracketing estimate of true pvalue:');
disp('(Since H and Hscore >= 0, P-values are inherently two-sided.)');
fprintf('    %7.5f    %7.5f    %7.5f\n', phat-pwid, phat, phat+pwid);

disp(' ');
disp('ABSOLUTE-VALUE DIFFERENCES BETWEEN ALL PAIRS:');
npairs=kk*(kk-1)/2;
disp('The following are MULTIPLE-COMPARISON-UNCORRECTED P-values:');
disp('For each of the  k*(k-1)/2  pairwise comparisons among the k treatment');
disp('  groups, consider the large-sample approximation of the');
disp('  Wilcoxon rank-sum test, including tie corrections:');
fprintf('Results for k(k-1)/2 = %d pairs:\n',npairs);

% Collect data for a 6-column table:
table=zeros(npairs,6);  ord=0;
zsobs=zeros(kk,kk);  wbuff=zeros(2*nmax+3,1);  % Buffer for two groups
maxzstar=0;  % Observed maximum Zstar statistic over pairs of treatments
% Copy data to rows of a matrix of fixed width
mdat=zeros(kk,nmax);
for j=1:n1;  mdat(1,j)=xx1(j);  end;
for j=1:n2;  mdat(2,j)=xx2(j);  end;
for j=1:n3;  mdat(3,j)=xx3(j);  end;
for j=1:n4;  mdat(4,j)=xx4(j);  end;
for j=1:n5;  mdat(5,j)=xx5(j);  end;
% Find pairwise Wilcoxon rank-sum P-values
for i=1:kk;  for j=i+1:kk;
  ni=num(i);  nj=num(j);  nb=ni+nj;  r=0;
  ord=ord+1;  table(ord,1)=i;  table(ord,2)=j;
  for u=1:ni;  r=r+1;  wbuff(r)=mdat(i,u);  end;
  for u=1:nj;  r=r+1;  wbuff(r)=mdat(j,u);  end;
  [midrij tiegrij] = getmidranks(wbuff(1:nb));
  % Wilcoxon score of first sample
  wx=0;  for u=1:ni;  wx=wx+midrij(u);  end;
  table(ord,3)=wx;
  wmean=ni*(nb+1)/2;  wvar=ni*nj*(nb+1)/12;  table(ord,4)=wmean;
  % Include the tiesum 
  tsum=0;
  for g=1:nb;
    if tiegrij(g)>0  tt=tiegrij(g);  tsum=tsum+(tt-1)*tt*(tt+1);
    else break;  end;  end;
  tcorr=tsum/((nb-1)*nb*(nb+1));
  zstar=abs(wx-wmean)/sqrt(wvar*(1-tcorr));
  zsobs(i,j)=zstar;
  table(ord,5)=zstar;  table(ord,6)=normpval(zstar);
  if zstar>maxzstar maxzstar=abs(zstar);  end;
  end;  end;

% Display the table
fprintf ('  #%d vs. #%d:  W=%5.1f  E(W)=%5.1f  Z=%6.3f   P=%6.4f (2-sided, UNCORR)\n', table');  %'
fprintf ('Maximum observed absolute Z-value:  %.4f\n', maxzstar);
ng1=0;  ng2=0;
for i=1:kk*(kk-1)/2;
  if table(i,6)<0.01 ng1=ng1+1;
  elseif table(i,6)<0.05  ng2=ng2+1;  end;  end;
if ng1>0
  fprintf ('Note that %d pair(s) have P<0.01 and %d more have P<0.05.\n',ng1,ng2);
else
  fprintf ('Note that none have P<0.01 and %d have P<0.05.\n',ng2);
end;

disp(' ');
disp('WHICH OF THESE ARE SIGNIFICANT, allowing for MULTIPLE COMPARISONS?');
disp('Using the Steel-Dwass-Crichlow-Fligner multiple-comparison');
disp('  procedure described in Section 6.5 of Hollander & Wolfe.');
disp('Let  ZObs(i,j)  be the pairwise Wilcoxon rank-sum score between all');
disp('  treatment-group pairs, normalized to have mean zero and variance one.');
disp('After a permutation of the midranks, let mstar = max_ab |Z(a,b)|');
disp('  where Z(a,b) is the corresponding normalized Wilcoxon score between');
disp('  pairs (a,b) for the permuted data. Finally, let P(i,j) be the');
disp('  proportion of permutations for which  mstar >= |ZObs(i,j)|.');
disp('Then P(i,j) is the probability that |Z(a,b)| >= |ZObs(i,j)| for at');
disp('  least one pair (a,b), which gives a multiple-comparison-corrected');
disp('  two-sided P-value for each ZObs(i,j).');

disp(' ');
disp('Simulating P-values for MULTIPLE-COMPARISON-CORRECTED COMPARISONS');
nsims=10000;
fprintf('Using %d permutations of the data\n',nsims);
disp('Pairwise MULTIPLE-COMPARISON-CORRECTED P-values:');

vals=zz(:,2);  Nge=zeros(kk,kk);
for ns=1:nsims;
    % Permute the values
    valperm=vals(randperm(nsum));
    % Redistribute to mdat array
    % This does not require that the rows of zz are sorted by groups
    % This results in a random permutation of the values in the data table
    r=0;  for i=1:kk;  for j=1:num(i);
            r=r+1;  mdat(i,j)=valperm(r);  end;  end;
    % Find mstar = max_ij zstar for this permutation
    mstar=0;
    for i=1:kk;  for j=i+1:kk;
      ni=num(i);  nj=num(j);  nb=ni+nj;  r=0;
      for u=1:ni;  r=r+1;  wbuff(r)=mdat(i,u);  end;
      for u=1:nj;  r=r+1;  wbuff(r)=mdat(j,u);  end;
      [midrij tiegrij] = getmidranks(wbuff(1:nb));
      % Wilcoxon score of first sample
      wx=0;  for u=1:ni;  wx=wx+midrij(u);  end;
      wmean=ni*(nb+1)/2;  wvar=ni*nj*(nb+1)/12;
      % Include the tiesum 
      tsum=0;
      for g=1:nb;
        if tiegrij(g)>0  tt=tiegrij(g);  tsum=tsum+(tt-1)*tt*(tt+1);
        else break;  end;  end;
      tcorr=tsum/((nb-1)*nb*(nb+1));
      zstar=abs(wx-wmean)/sqrt(wvar*(1-tcorr));
      if zstar>mstar mstar=zstar;  end;
      end;  end;
    %% Which pairs are farther apart than mstar?
    for i=1:kk;  for j=i+1:kk;
      if mstar>=zsobs(i,j)  Nge(i,j)=Nge(i,j)+1;  end;  end;  end;
    end;


ord=kk*(kk-1)/2;  table=zeros(ord,6);  r=0;
for i=1:kk;  for j=i+1:kk;
  r=r+1;  table(r,1)=i;  table(r,2)=j;  table(r,3)=zsobs(i,j);
  table(r,4)=Nge(i,j);  table(r,5)=nsims;  table(r,6)=Nge(i,j)/nsims;
  end;  end;

fprintf('  #%d vs #%d:  Z=%6.3f  P=%5d/%5d = %7.5f\n',table');  %'

ng1=0;  ng2=0;
for i=1:kk*(kk-1)/2;
  if table(i,6)<0.01 ng1=ng1+1;
  elseif table(i,6)<0.05  ng2=ng2+1;  end;  end;
if ng1>0
  fprintf ('Note that %d have P<0.01 and %d more have P<0.05.\n',ng1,ng2);
else
  fprintf ('Note that none have P<0.01 and %d have P<0.05.\n',ng2);
end;





