
disp('UMBRELLA TESTS for one-way layouts');
disp('Test whether or not medians of treatment groups increase to that of');
disp('  a particular treatment group and then decrease afterwards');
disp('Two versions: Given peak and unknown peak');

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

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

disp('SECOND DATA SET:');
disp('Fasting metabolic rates in deer by month pairs (Jan-Feb, Mar-Apr, etc)');
disp('Table 6.8 p215 in Hollander & Wolfe 2nd edn');

xx1 = [ 36.0  33.6  26.9  35.8  30.1  31.2  35.3 ];
xx2 = [ 39.9  29.1  43.4 ];
xx3 = [ 44.6  54.4  48.2  55.7  50.0 ];
xx4 = [ 53.8  53.9  62.5  46.6 ];
xx5 = [ 44.3  34.1  35.7  35.6 ];
xx6 = [ 31.7  22.1  30.7 ];

kk=6;  % Number of treatment groups
nn = [length(xx1)  length(xx2)  length(xx3)  length(xx4)  length(xx5)  length(xx6)];
ntot = sum(nn);

fprintf('%d treatment groups with %d total observations:\n',kk,ntot);
% Copy data to rows of a matrix of fixed width
mdat=zeros(kk,ntot);
for j=1:nn(1);  mdat(1,j)=xx1(j);  end;
for j=1:nn(2);  mdat(2,j)=xx2(j);  end;
for j=1:nn(3);  mdat(3,j)=xx3(j);  end;
for j=1:nn(4);  mdat(4,j)=xx4(j);  end;
for j=1:nn(5);  mdat(5,j)=xx5(j);  end;
for j=1:nn(6);  mdat(6,j)=xx6(j);  end;
% Display data by rows of mdat
for i=1:kk;
  fprintf('  #%d (n=%d):', i,nn(i));
  for j=1:nn(i);  fprintf('  %4.1f',mdat(i,j));  end;  fprintf('\n');
  end;

% Put in ntot rows with group numbers
zz = [ xx1'    ones(nn(1),1);
       xx2'  2*ones(nn(2),1);
       xx3'  3*ones(nn(3),1);
       xx4'  4*ones(nn(4),1);
       xx5'  5*ones(nn(5),1);
       xx6'  6*ones(nn(6),1) ];

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

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

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

nn2=zeros(1,kk);  dsum=zeros(1,kk);  ranksum=zeros(1,kk);
for i=1:ntot;  val=zz(i,2);  grp=zz(i,3);  rnk=zz(i,4);
  nn2(grp)=nn2(grp)+1;  dsum(grp)=dsum(grp)+val;
  ranksum(grp)=ranksum(grp)+rnk;  end;
% Data average and rank average for each treatment group
% (Componentwise division, so OK if treatment group sizes are different)
DataAv=dsum./nn2;
RankAv=ranksum./nn2;

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

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

disp(' ');
disp('THE UMBRELLA TEST for KNOWN p=4 (Mack-Wolfe):');
disp('Single-umbrella test:');
disp('Is the data significant for a `peak` in treatment-group medians at p=4?');

disp(' ');
disp('First Step:');
disp('Matrix of Mann-Whitney differences between treatment groups:');
% Make sure main data-set is sorted on group
zz=sortrows(zz,[3 1]);
uu=zeros(kk);  val=zz(:,2);  grp=zz(:,3);
for i=1:ntot;  for j=i+1:ntot;  % Only if grp(i)<grp(j)
  if grp(i)<grp(j)
    uu(grp(i),grp(j)) = uu(grp(i),grp(j)) + (sign(val(j)-val(i))+1)/2;
    end;  end;  end;
% Fill in below the diagonal
for i=1:kk;  for j=i+1:kk;
  uu(j,i)=nn(i)*nn(j)-uu(i,j);  end;  end;
% Put -1 on the diagonal
for i=1:kk;  uu(i,i)=-1;  end;

vv=[ [1:kk]'  uu];  %'
fprintf (' %2d.   %4.1f   %4.1f   %4.1f   %4.1f   %4.1f   %4.1f\n', vv');  %'

disp(' ');
disp('Example of an Ap statistic:');
disp('  Ap(4) = U(1,2)+U(1,3)+U(1,4)+U(2,3)+U(2,4)+U(3,4) + U(5,4)');

disp(' ');
disp('The umbrella statistics Ap with normal approximations Astar:');
Ap=zeros(kk,1);
for k=1:kk;  au=0;  % Compute Ap(k)
  for i=1:k-1;  for j=i+1:k;  au=au+uu(i,j);  end;  end;
  for i=k+1:kk;  for j=k:i-1;  au=au+uu(i,j);  end;  end;
  Ap(k)=au;  end;

% Compute means and variances for Ap(k) and normal approximation Astar(k)
% Constants used for all k:
nn2=0;  nn3=0;
for i=1:kk;  ns=nn(i);
  nn2=nn2+ns*ns;  nn3=nn3+ns*ns*(2*ns+3);  end;

MAp=zeros(kk,1);  Varp=zeros(kk,1);  Astar=zeros(kk,1);
for pp=1:kk;  np=nn(pp);
  tn1=0;  tn2=0;
  for i=1:pp;  tn1=tn1+nn(i);  end;
  for i=pp:kk;  tn2=tn2+nn(i);  end;
  tsq = tn1*tn1 + tn2*tn2;
  tcb = tn1*tn1*tn1 + tn2*tn2*tn2;
  MAp(pp) = (tsq - nn2 - np*np)/4;
  var = 2*tcb + 3*tsq - nn3 - np*np*(2*np+3) + 12*np*(tn1*tn2 - np*ntot);
  Varp(pp) = var/72;
  Astar(pp) = (Ap(pp) - MAp(pp))/sqrt(Varp(pp));
  end;

vv = [ [1:kk]'  Ap  MAp  Varp  Astar  1-normint(Astar)];   %'

disp('(P-values are one-sided normal.)');
fprintf(' p=%d  Ap=%5.1f   Mean=%5.1f  Var=%7.2f   Astar=%6.3f   P=%7.5f\n', vv');  %'

disp(' ');  pp=4;
fprintf('Simulating the P-value of Ap=%g given p=%d using %d permutations:\n',Ap(pp),pp,nsims);

Nge=0;  ApObs=Ap(pp);  val=zz(:,2);  grp=zz(:,3);
for ns=1:nsims;
  % valp are permuted data values, uu is local Mann-Whitney-difference matrix
  valp=val(randperm(ntot));  uu=zeros(kk);
  for i=1:ntot;  for j=i+1:ntot;  % Only if grp(i)<grp(j)
    if grp(i)<grp(j)
      uu(grp(i),grp(j)) = uu(grp(i),grp(j)) + (sign(valp(j)-valp(i))+1)/2;
    end;  end;  end;
  % Fill in below the diagonal
  for i=1:kk;  for j=i+1:kk;
    uu(j,i)=nn(i)*nn(j)-uu(i,j);  end;  end;
  % Compute App=Ap(pp) for permuted data
  App=0;
  for i=1:pp-1;  for j=i+1:pp;  App=App+uu(i,j);  end;  end;
  for i=pp+1:kk;  for j=pp:i-1;  App=App+uu(i,j);  end;  end;
  if App>=ApObs  Nge=Nge+1;  end;
  end;

disp('Number greater-or-equal to old Ap(p) and number of permuations:');
fprintf('    Nge=%d    nsims=%d\n', Nge,nsims);
disp('95% Conf.Int. for one-sided P-value bracketing estimate:');
phat=Nge/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp([phat-pwid  phat  phat+pwid]);

disp(' ');
disp('THE UMBRELLA TEST with UNKNOWN p (Chen-Wolfe):');
disp(' ');
disp('Multiple-umbrella test:');
disp('Is there a `peak` in treatment-group medians at SOME p, and if so,');
disp('  at which treatment group?');
disp(' ');
disp('Note that the one-sided P-values in the table above are NOT');
disp('  multiple-comparison corrected.');

disp(' ');
[Max_Astar pp0] = max(Astar);
fprintf('Max_Astar=%g  is attained at  p=%d\n', Max_Astar,pp0);
fprintf('Simulating the P-value of Max_Astar=%g using %d permutations:\n',Max_Astar,nsims);

Nge=0;  ApObs=Ap(pp);  val=zz(:,2);  grp=zz(:,3);
for ns=1:nsims;
  % valp are permuted data values, uu is local Mann-Whitney-difference matrix
  valp=val(randperm(ntot));  uu=zeros(kk);
  for i=1:ntot;  for j=i+1:ntot;  % Only if grp(i)<grp(j)
    if grp(i)<grp(j)
      uu(grp(i),grp(j)) = uu(grp(i),grp(j)) + (sign(valp(j)-valp(i))+1)/2;
    end;  end;  end;
  % Fill in below the diagonal
  for i=1:kk;  for j=i+1:kk;
    uu(j,i)=nn(i)*nn(j)-uu(i,j);  end;  end;
  MaxRand=0;
  % Compute Ap(p) (for all p) for permuted data
  for p=1:kk;
    App=0;
    for i=1:p-1;  for j=i+1:p;  App=App+uu(i,j);  end;  end;
    for i=p+1:kk;  for j=p:i-1;  App=App+uu(i,j);  end;  end;
    zz=(App-MAp(p))/sqrt(Varp(p));
    if (p==1 || zz>MaxRand) MaxRand=zz;  end;  end;
  if MaxRand>=Max_Astar  Nge=Nge+1;  end;
  end;

disp(' ');
disp('Number greater-or-equal to Max_Astar and number of permuations:');
disp([Nge nsims]);
disp('95% Conf.Int. for one-sided P-value bracketing estimate:');
phat=Nge/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
fprintf('    %7.5f   %7.5f   %7.5f\n',phat-pwid, phat, phat+pwid);



