
disp('Multiple Comparison with a Control for One-Way layouts');
disp('Control-other 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('The first treatment (the control) are normal subjects.');
disp('Hollander & Wolfe (Table 6.4, p201)');
disp('See OneWayMultComp.txt for more detail and more analyses');

%% 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
nn = sum(num);     % Total number of observations
nmax = max(num);   % Maximum number of observations per group

disp(' ');
fprintf('Grp #1 (n=%2d) Normal subjects (Control)\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,nn);

% Put data in nn rows with group numbers
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:nn]'  zz  midr  tiegr];  %'

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');  %'

disp(' ');
disp('See OneWayMultComp for a table of observations and midranks');
disp('  and for the Kruskal-Wallis test.');

disp(' ');
disp('MULTIPLE COMPARISONS WITH A CONTROL:');
disp('Which of the k-1  non-control treatment groups are significantly');
disp('  LARGER than the control, allowing for multiple comparisons?');
disp('(We could also ask for significantly different in absolute value, or');
disp('  even significantly smaller, but here we consider significantly greater.)');
disp('See Section 6.7 in the text.');

disp(' ');
disp('As in OneWayMultComp.m, we use normalized Wilcoxon rank-sum differences,');
disp('  and base a multiple-comparison-corrected P-value for treatment #i on');
disp('  the number of times that the maximum of treatments #j vs. control group');
disp('  normalized comparisons is LARGER THAN the observed treatment #i vs');
disp('  control comparison. Since this is a one-sided procedure, absolute');
disp('  values are not used.');
disp('NOTE THAT this is the method of Steel (1959) (text p259) and NOT');
disp('  the Nemenyi-Damico-Wolfe method discussed in most of Section 6.7.');

disp(' ');
disp('For comparison, the following are MULTIPLE-COMPARISON-UNCORRECTED P-values.');
disp('One-way Wilcoxon W and Z values:');

% We proceed as before, except that we compare all treatment groups
%   with the control, and consider only the probability of being larger
%   than the control.
% We assume that the control group comes first,

ncon=num(1);    % Control-group sample size
wbuff=zeros(2*nmax+3,1);  % Buffer for values from two treatment groups
% Create a table with 5 columns
table=zeros(kk-1,5);  ord=0;
zsobs=zeros(1,kk);
% 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 one-sided pairwise Wilcoxon rank-sum P-values
% Note that maxzstar below can be negative
for i=2:kk;    % All but the control group
  ni=num(i);  nb=ni+ncon;  r=0;
  ord=ord+1;  table(ord,1)=i;
  for u=1:ni;  r=r+1;  wbuff(r)=mdat(i,u);  end;
  % Recall that xx1() is the control
  for u=1:ncon;  r=r+1;  wbuff(r)=xx1(u);  end;
  [midri tiegri] = getmidranks(wbuff(1:nb));
  % Wilcoxon score of the non-control sample
  wx=0;  for u=1:ni;  wx=wx+midri(u);  end;
  table(ord,2)=wx;
  wmean=ni*(nb+1)/2;  table(ord,3)=wmean;
  wvar=ni*ncon*(nb+1)/12;
  % Include the tiesum 
  tsum=0;
  for g=1:nb;
    if tiegri(g)>0  tt=tiegri(g);  tsum=tsum+(tt-1)*tt*(tt+1);
    else break;  end;  end;
  tcorr=tsum/((nb-1)*nb*(nb+1));
  zstar=(wx-wmean)/sqrt(wvar*(1-tcorr));
  zsobs(i)=zstar;
  table(ord,4)=zstar;  table(ord,5)=1-normint(zstar);
  % This initializes maxzstar at zstar if i==2
  if (i==2 || zstar>maxzstar) maxzstar=zstar;  end;
  end;

fprintf ('  #%d vs. #1:  W=%5.1f  E(W)=%5.1f  Z=%6.3f   P=%6.4f (1-sided, UNCORR)\n', table');  %'
fprintf ('Maximum observed Z-value:  %.4f\n', maxzstar);
ng1=0;  ng2=0;
for i=1:kk-1;
  if table(i,5)<0.01 ng1=ng1+1;
  elseif table(i,5)<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;


disp(' ');
disp('Simulating P-values for MULTIPLE-COMPARISON-CORRECTED comparisons');
fprintf('  for groups #2-#%d vs control:\n',kk);

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

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

disp('One-sided MULTIPLE-COMPARISON-CORRECTED P-values:');

% vals are the values
vals=zz(:,2);  Nge=zeros(1,kk);
for ns=1:nsims;
    % Permute the values
    valperm=vals(randperm(nn));
    % 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_i zstar for this permutation
    for i=2:kk;
      ni=num(i);  nb=ni+ncon;  r=0;
      % Copy group #i followed by control group to wbuff
      for u=1:ni;  r=r+1;  wbuff(r)=mdat(i,u);  end;
      for u=1:ncon;  r=r+1;  wbuff(r)=mdat(1,u);  end;
      % Carry out Wilcoxon rank-sum procedure
      [midri tiegri] = getmidranks(wbuff(1:nb));
      % Wilcoxon score of the non-control sample
      wx=0;  for u=1:ni;  wx=wx+midri(u);  end;
      wmean=ni*(nb+1)/2;  wvar=ni*ncon*(nb+1)/12;
      % Include the tiesum 
      tsum=0;
      for g=1:nb;
        if tiegri(g)>0  tt=tiegri(g);  tsum=tsum+(tt-1)*tt*(tt+1);
        else break;  end;  end;
      tcorr=tsum/((nb-1)*nb*(nb+1));
      zstar=(wx-wmean)/sqrt(wvar*(1-tcorr));
      if (i==2 || zstar > mstar) mstar=zstar;  end;
      end;
    %% Which treatment groups are larger than this?
    for i=2:kk;
      if mstar>=zsobs(i)  Nge(i)=Nge(i)+1;  end;  end;
    end;

% Build table for final MC-corrected table

ord=0;
for i=2:kk;
  ord=ord+1;  table(ord,1)=i;  table(ord,2)=zsobs(i);
  table(ord,3)=Nge(i);  table(ord,4)=nsims;  table(ord,5)=Nge(i)/nsims;
  end;

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

ng1=0;  ng2=0;
for i=1:kk-1;
  if table(i,5)<0.01 ng1=ng1+1;
  elseif table(i,5)<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;

