
% Two-Way Layout with Replications
%
% Data  X(i,q,j)  for 1 le i le nn (blocks)
%           1 le q le cc (replications)
%           1 le j le kk (treatment groups)
%
%   assuming cc replications at each (i,j) combination
%
% Assign Friedman-like ranks R(i,q,j) to all cc*kk observations
%   within each block, so that  1 le R(i,q,j) le cc*kk
%
% The assumption of exactly cc replications for each cell
%   or (i,j) combination is needed only for the large-sample
%   chi-square approximation.
%
% The permutation-test below can be done for any two-way layout
%   with repeated or missing data
%

disp('Two-Way Layout with Multiple Replications per Cell');
disp('Data  X(i,q,j)  where');
disp('   1 le i le nn   Subjects or Blocks');
disp('   1 le q le cc   Replications for (i,j) combination');
disp('   1 le j le kk   Treatments or Treatment Groups');

disp(' ');
disp('Example:');
disp('Test for consistency across 4 laboratories for amounts');
disp('  of niacin detected in enriched bran flakes.');
disp('Each laboratory tests 3 samples of bran flakes for each');
disp('  of 3 levels of known additional enrichment (0mg, 4mg, 8mg per 100g)');
disp('Data source:  Campbell and Pelletier (1962)');
disp('See Table 7.20 in text (p331)');

xxr = [ 1  7.58  7.87  7.71   11.63 11.87 11.40   15.00 15.92 15.58 ;
        2  8.00  8.27  8.00   12.20 11.70 11.80   16.60 16.40 15.90 ;
        3  7.60  7.30  7.82   11.04 11.50 11.49   15.87 15.91 16.28 ;
        4  8.03  7.35  7.66   11.50 10.10 11.70   15.10 14.80 15.70  ];

[kkr bbr] = size(xxr);
xx = xxr(:,2:bbr);       % Data without lab numbers
[kk bb] = size(xx);
nn = 3;     % Number of blocks
cc = bb/3;  % Number of replications per block per treatment group

xxx=zeros(nn,cc,kk);
for i=1:nn;  for j=1:kk;  for q=1:cc;
   xxx(i,q,j) = xx(j,(i-1)*cc+q);  end;  end;  end;

disp(' ');
disp('Data by blocks:');
%% Column headings
fprintf('    ');
for j=1:kk;  fprintf('        Lab%d      ',j);  end;  fprintf('\n');
%% Data
for i=1:nn;
  fprintf('%dmg:',(i-1)*4);
  for j=1:kk;  for q=1:cc;
    if q==1 && j>1
      fprintf('|%5.2f',xxx(i,q,j));
    else
      fprintf(' %5.2f',xxx(i,q,j));
    end; end;  end;
  fprintf('\n');  end;

%% Assign ranks within each block
rrr=zeros(nn,cc,kk);
xvals=zeros(kk*cc,1);  rvals=zeros(kk*cc,1);
for i=1:nn;  r=0;
  %% Copy block #i into a (kk*cc)x1 column vector
  for j=1:kk;  for q=1:cc;
    r=r+1;  xvals(r)=xxx(i,q,j);  end;  end;
  %% Assign midranks for the column vector
  [rvals tiegr] = getmidranks(xvals);  r=0;
  %% Copy midranks for block #i into nnxccxkk matrix
  for j=1:kk;  for q=1:cc;
    r=r+1;  rrr(i,q,j)=rvals(r);  end;  end;
  end;

disp(' ');
disp('Within-block ranks by blocks:');
%% Column headings
fprintf('   ');
for j=1:kk;  fprintf('        Lab%d      ',j);  end;  fprintf('\n');
%% Data
for i=1:nn;
  fprintf('%dmg:',(i-1)*4);
  for j=1:kk;  for q=1:cc;
    if q==1 && j>1
      fprintf('|%4.1f ',rrr(i,q,j));
    else
      fprintf(' %4.1f ',rrr(i,q,j));
  end;  end;  end;
  fprintf('\n');  end;

% Compute values for Summary Table
%% rrr is nnxccxkk
%% This sums out the first two dimensions
ranksum=sum(sum(rrr));
rankavg=ranksum/(nn*cc);
xsum=sum(sum(xxx));
xavg=xsum/(nn*cc);
disp(' ');
disp('Summaries by treatment group:');
disp('         Xavg Ranksum  Rankavg');
for j=1:kk;
  fprintf(' Lab%d: %6.2f  %5.1f   %6.3f\n', ...
    j,xavg(j),ranksum(j),rankavg(j));  end;

Lscore=sum(ranksum.*ranksum);
fprintf('\nLscore = Sum ranksum^2 for permutation test: %g\n',Lscore);

% Compute Mack-Skillings statistic
rankc = ranksum - nn*cc*(kk*cc+1)/2;
mmss = (12/(cc*cc*kk*nn*(kk*cc+1)))*sum(rankc.*rankc);
df=kk-1;
Pval=pvchisq(mmss,df);
disp(' ');
disp('Large-sample approximation:');
fprintf ('MS=%g  Pval=%7.5f  (df=%d)\n', mmss,Pval,df);


%% Find observed rank-average differences for 
%%   multiple-comparison-corrected pairwise differences
obsdiff=zeros(kk*kk,1);
r=0;
for j1=1:kk;  for j2=j1+1:kk;
  r=r+1;  obsdiff(r)=abs(ranksum(j1)-ranksum(j2));
  end;  end;
nobs=r;
obsdiff=obsdiff(1:nobs);


nsims=100000;

fprintf('\nCarrying out Nsims=%d Friedman-like permutations:\n',nsims);

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

disp('Lscore = Sum of Squares over Treatment Groups');
fprintf('  of Treatment-Group Rank Sums:  Observed Lscore=%g\n',Lscore);

% Do nsims permutations for a test of Lscore
rrp=rrr;   % Preserve original rank matrix
rvals=zeros(kk*cc,1);
Ngr=0;  Ngrp=zeros(nobs,1);
for ns=1:nsims;
  %% Permute blocks independently
  %% Consider rrs as an nn x (cc*kk) matrix
  %%   then permute each row (of length cc*kk)
  rrs=reshape(rrp,nn,cc*kk);
  for i=1:nn;  rrs(i,:) = rrs(i,randperm(cc*kk));  end;
  %% Restore shape of rrp
  rrp=reshape(rrs,nn,cc,kk);
  %% Treatment-group rank sums are found by summing out
  %%   the first two dimensions of rrp(i,q,j)
  rsum=sum(sum(rrp));
  %% Sum of squares of the treatment-group rank sums
  Rscore=sum(rsum.*rsum);
  if Rscore >= Lscore  Ngr=Ngr+1;  end;
  %% Find maximum pairwise difference
  rsum=sort(rsum);
  maxdiff=rsum(kk)-rsum(1);  %% Last - First
  %% Apply to counts for all-pair difference
  %% This is applied for each r for 1 <= r <= nobs=kk(kk-1)/2
  Ngrp = Ngrp + (maxdiff>=obsdiff);
  end;

disp('Number of simulations with values >= Lscore and total number:');
disp(' ');
disp([Ngr nsims]);

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

disp(' ');
disp('Rank-sum (and rank-average) differences between pairs of treatments');
disp('  and pairwise P-values multiple-comparison-corrected for all pairs:');
disp('  Based on maximum pairwise rank-sum differences');
fprintf('    in N=%d simulations\n',nsims);
r=0;
for j1=1:kk;  for j2=j1+1:kk;  r=r+1;
  fprintf('  Lab%d v Lab%d:  %6.2f  %6.3f   %5d/%d = %6.4f\n', ...
    j1,j2, obsdiff(r), obsdiff(r)/(nn*cc), Ngrp(r),nsims, Ngrp(r)/nsims);
  end;  end;

fprintf('\nLabs 1,3,4 appear similar.  Lab2 differs from Labs 1,3,4\n');
fprintf('  but the difference is MCC significant only for Lab2 v Labs 3,4\n');


