
%  Lehmann`s method of aligned ranks for two-way layouts:
%
%  For data  X_ij  for n blocks or subjects (1 le i le n)
%    and k treatments (1 le j le k):
%
%  Friedman`s test makes no assumptions about the distribution of data
%    within blocks, and loses power because it doesn`t use within-block
%    comparisons. Thus it is closer in spirit to the sign test than to
%    the Wilcoxon signed-rank test. (In fact, the Friedman test reduces
%    to the sign test for k=2.)
%
%  Lehmann`s procedure effectively assumes that both blocks and treatment
%    groups vary by additive shifts so that between-block variation
%    can be taken into account, as in the Wilcoxon signed-rank test.
%    It make more assumptions than the Friedman test, but has greater
%    asymptotic efficiency for most distributions of data.
%
%  The method is as follows:
%     (i) CENTER X_ij WITHIN BLOCKS by setting  Xc_ij = X_ij - Xbar_i
%        (That is, by subtracting block means. Lehmann says that you can
%          also use within-block medians)
%    (ii) Use  Xc_ij  to define Kruskal-Wallis ranks (1 <= R_ij <= nk)
%   (iii) The score is essentially  L = Sum(j=1,k) R_j^2  for 
%       treatment-group sums  R_j = Sum(i=1,nn) R_ij , exactly as in
%       the Friedman and Kruskal-Wallis tests. The significance of the
%       observed score is determined by sets of within-block permutations
%       as in the Friedman test. That is, the aligned-rank procedure
%       uses Kruskal-Wallis ranks with Friedman permutations.
%    (iv) There is a large-sample approximation based on the Hajek Central
%       Limit Theorem: If
%
%          Q = Sum(j=1,k) (R_j - n(nk+1)/2)^2 / (Sum(i=1,n) S^2_i)
%
%       where S^2_i = (1/(k-1))Sum(j=1,k) (R_ij - Rbar_i)^2  is the
%       sample variance of the Kruskal-Wallis ranks within the i-th block,
%       then Q is asymptotically chi-square with k-1 degrees of freedom
%       under Friedman permutations. Note that the denominator of Q is
%       invariant under within-block permutations, so that Q and L are
%       permutation-equivalent scores. No tie correction is needed since
%       it is already built into the denominator of Q.
%


disp('Lehmann`s mathod of aligned ranks for two-way layouts');
disp('Data  X(i,j)  where');
disp('   1 le i le nn   Subjects or Blocks');
disp('   1 le j le kk   Treatments or Treatment Groups');
disp('Model:  x(i,j) = theta_i + tau_j + error');
disp('Testing H_0: tau_j=tau_1,  so no treatment effect');

disp(' ');
disp('Example:  Responses to suggestion under hypnosis as measured by');
disp('  electrical skin potential in millivolts');
disp('From:  Lehmann, ``Nonparametrics: Statistical Methods Based on Ranks``');
disp('  page 264 (see also pp264-271)');
disp('k=4 treatment groups and n=8 subjects');
disp('Treatment groups:');
disp('   1. Fear');
disp('   2. Happiness');
disp('   3. Depression');
disp('   4. Calmness');

xxr = [ 1   23.1   22.7   22.5   22.6 ;
        2   57.6   53.2   53.7   53.1 ;
        3   10.5    9.7   10.8    8.3 ;
        4   23.6   19.6   21.1   21.6 ;
        5   11.9   13.8   13.7   13.3 ;
        6   54.6   47.1   39.2   37.0 ;
        7   21.0   13.6   13.7   14.8 ;
        8   20.3   23.6   16.3   14.8   ];

[nnr kkr] = size(xxr);
xx = xxr(:,2:kkr);       % Data without row numbers
[nn kk] = size(xx);

fprintf('\nn=%d blocks (subjects) and k=%d treatment groups\n',nn,kk);

disp(' ');
disp('Subject, Fear, Happiness, Depression, Calmness');
fprintf('  %2d.   %5.2f   %5.2f   %5.2f   %5.2f\n',xxr');  %'

disp(' ');
disp('First, let`s apply Friedman`s test:');
disp('Data with (within-block) Friedman ranks:');

rr=zeros(nn,kk);  tg=zeros(nn,kk);
sbuff=zeros(kk,1);  % Column-vector buffer
for i=1:nn;  % For all replications
  % Put values for that block into sbuff
  for j=1:kk;  sbuff(j)=xx(i,j);  end;
  [midr  tiegr] = getmidranks(sbuff);
  for j=1:kk;  rr(i,j)=midr(j);  tg(i,j)=tiegr(j);  end;
  end;

% Buffer for a table with row numbers then kk sets of data then ranks
xxd = zeros(nn,2*kk+1);

xxd(:,1) = [1:nn]';  %'
for j=1:kk;  xxd(:,2*j)=xx(:,j);  xxd(:,2*j+1)=rr(:,j);  end;

% Note continuation line
fprintf(' %2d.  %5.2f(%3.1f)  %5.2f(%3.1f)  %5.2f(%3.1f)  %5.2f(%3.1f)\n',...
  xxd');  %'

disp(' ');
disp('Rank sums and rank averages for each treatment group:');
rsum=sum(rr,1);   % 1xkk row vector of column sums of nnxkk matrix rr
ravg=rsum/nn;     % Also a 1xkk row vector
fprintf('   %2d. %7.2f   %6.2f\n',[ [1:kk]'  rsum'  ravg']');

%% Carry out Friedman test for comparison
%% Centered rank averages
ravgc = ravg - (kk+1)/2;
ss = (12*nn/(kk*(kk+1)))*sum(ravgc.*ravgc);
df=kk-1;
disp(' ');
disp('Friedman test: With no tie correction:');
Pval = pvchisq(ss,df);
fprintf('  S=%7.5f   P= %6.4f  (df=%d)\n', ss,Pval,df);

% Compute tie correction
tiesum=0;
for i=1:nn;  tgg=0;
  for j=1:kk;  tt=tg(i,j);
    if tt>1 tgg=tgg+(tt-1)*tt*(tt+1);  end;  end;
  tiesum=tiesum+tgg;  end;
tiecorr = tiesum/(nn*(kk-1)*kk*(kk+1));
ssp = ss/(1-tiecorr);
disp(' ');
disp('With tie correction:');
fprintf('  Tiesum=%g   Tiecorr=%g\n', tiesum,tiecorr);
Pval = pvchisq(ssp,df);
fprintf('  SP=%7.5f   P= %6.4f  (df=%d)\n', ssp,Pval,df);

nsims=10000;

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

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

Fscore = sum(rsum.*rsum);
disp('Fscore = Sum of Squares over Treatment Groups');
fprintf('  of Treatment-Group Rank Sums:  Observed Fscore=%g\n',Fscore);
 
% Do nsims permutations for a test of Fscore
rrp=rr;   % Preserve original ranks
Ngr=0;
for ns=1:nsims;
  %% Permute rows of rrp independently
  for i=1:nn;  rrp(i,1:kk) = rrp(i,randperm(kk));  end;
  %% For a matrix rrp, sum(rrp,1) is the row vector of column sums
  %% Sum of squares of the sum of the columns of rrp
  Rscore=sum(sum(rrp,1).*sum(rrp,1));
  if Rscore >= Fscore  Ngr=Ngr+1;  end;
  end;

disp('Number of simulations with values >= Fscore 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 H and Fscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);

disp(' ');
disp('Is the Friedman test the most appropriate?');
disp('Data again with RANGE (Max-Min for each subject):');
disp('Note that the RANGE is highly variable across subjects.');
disp('The Friedman test does not take RANGE into account.');
disp(' ');
disp('Subject, Fear, Happiness, Depression, Calmness, (RANGE, rank)');
range=zeros(nn,1);
for i=1:nn;  range(i)=max(xx(i,:))-min(xx(i,:));  end;
[midr tiegr] = getmidranks(range);
fprintf('  %2d.   %5.2f   %5.2f   %5.2f   %5.2f   (%5.2f) (%g)\n',...
  [xxr range midr]');  %'

disp(' ');
disp('Perhaps we should weight higher those subjects with higher range?');
disp('The Wilcoxon signed-rank test effectively does this.');
disp('The following procedure accomplishes this, but does not reduce');
disp('  to the Wilcoxon signed-rank test if k=2.');

disp(' ');
disp('Data centered at block (subject) means:');

%% xmean is an  nnxkk  matrix whose rows are identically equal
%%   to the row means of xx
xmean=xx*ones(kk,kk)/kk;
xxc = xx-xmean;   %% xxc is xx centered at row means

%% Find the range of the rows of xxc for th next display
%% The values should be the same as for xx
for i=1:nn;  range(i)=max(xxc(i,:))-min(xxc(i,:));  end;

disp('Subject, Fear, Happiness, Depression, Calmness, RANGE, Subject mean');
fprintf('  %2d.   %5.2f   %5.2f   %5.2f   %5.2f   (%5.2f)  (%5.2f)\n', ...
  [ [1:nn]' xxc range xmean(:,1)]');

disp(' ');
disp('Data centered at block means with Kruskal-Wallis ranks and RANK RANGE:');

% Put centered data in a single Nx1 column vector where N=nn*kk
% MATLAB reshapes matrices by columns
xxc_tall = reshape(xxc,nn*kk,1);

% Find midranks for the Nx1 column vector
% The tie-groups will not be needed
[rkw_tall tiegr] = getmidranks(xxc_tall);

% Reshape the Nx1 column vector of ranks to an nnxkk matrix
%   of Kruskal-Wallis ranks
rkw = reshape(rkw_tall,nn,kk);

% Display a table of centered values with Kruskal-Wallis ranks
xxcd=zeros(nn,2*kk+1);
% The first column is row numbers
% Pairs of columns after the first are data, then ranks
xxcd(:,1)=xxr(:,1);
for j=1:kk;  xxcd(:,2*j)=xxc(:,j);  xxcd(:,2*j+1)=rkw(:,j);  end;

% Also, find ranges of the Kruskal-Wallis ranks within each block
for i=1:nn;  range(i)=max(rkw(i,:))-min(rkw(i,:));  end;

% The final table - note continuation line
fprintf(...
  '  %2d.  %5.2f(%4.1f)  %5.2f(%4.1f)  %5.2f(%4.1f)  %5.2f(%4.1f)  (%4.1f)\n', ...
  [xxcd range]'); %'

disp(' ');
disp('Note that RANGE is still highly variable, but not as variable for the');
disp('  Kruskal-Wallis ranks as for the original data.');

disp(' ');
disp('Rank sums and rank averages for each treatment group:');
disp('Rank averages appear to be more variable than with Friedman ranks');
rsum=sum(rkw,1);   % 1xkk row vector of column sums of nnxkk matrix rkw
ravg=rsum/nn;      % Also a 1xkk row vector
fprintf('   %2d. %7.2f   %6.2f\n',[ [1:kk]'  rsum'  ravg']');

% Compute Q and find large-sample P-value
rmean=nn*(nn*kk+1)/2;   % Expected value of each rsum(j)
rsumc=rsum-rmean;
qq = sum(rsumc.*rsumc);

% Denominator of qq is rkw-sum of sample variances of blocks
bmean=rkw*ones(kk,kk)/kk;  % nnxkk matrix of row means
rkwc=rkw-bmean;
% sum(sum(...)) is necessary for a double sum over a matrix
lvar=sum(sum(rkwc.*rkwc))/(kk-1);
qq = qq/lvar;
Pval=pvchisq(qq,df);

disp(' ');
disp('Large-sample Aligned Rank test:');
fprintf('\n  Q=%g   P=%6.4f  (df=%d)\n', qq,Pval,df);
  
% The sums of squares of Aligned-Rank treatment rank sums
ALobs = sum(rsum.*rsum);
fprintf('\nCarrying out Nsims=%d Friedman permutations\n',nsims);
fprintf('  using ALscore = Sum of Squares over Treatment Groups\n');
fprintf('  of Treatment-Group Rank Sums:  Observed ALobs=%g\n', ALobs);
 
% Do nsims permutations for a test of Fscore
rrp=rkw;   % Preserve original ranks
Ngr=0;
for ns=1:nsims;
  %% Permute rows of rrp independently
  for i=1:nn;  rrp(i,1:kk) = rrp(i,randperm(kk));  end;
  %% For a matrix rrp, sum(rrp,1) is the row vector of column sums
  %% Sum of squares of the sums of the columns of rrp
  Rscore=sum(sum(rrp,1).*sum(rrp,1));
  if Rscore >= ALobs  Ngr=Ngr+1;  end;
  end;

disp('Number of simulations with values >= ALobs 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 H and Fscore >= 0, P-values are inherently two-sided.)');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);

disp(' ');
disp('Thus these data are significant for Lehmann`s Aligned Rank test,');
disp('  but not for Friedman`s test.');


