
disp('Two-Way Layout:');
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) = mu(i,j) plus error  or mu + theta_i + tau_j plus error');
disp('Test H_0:  mu(i,j) = theta_i  and there is no effect of Treatment');

disp(' ');
disp('Example 2:');
disp('Strength of cotton fibers as a function of levels of potash (K_2O) in field');
disp('k=5 potash fields with n=3 randomized blocks');
disp('Potash levels in pounds per acre:');
disp('  144(k=1)   108(k=2)   72(k=3)   54(k=4)  36(k=5)');
disp('Table 7.5 (p286) in text');

xxr = [ 1   7.46   7.17   7.76   8.14   7.63;
        2   7.68   7.57   7.73   8.15   8.00;
        3   7.31   7.80   7.74   7.87   7.93 ];

[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('Cotton fiber strengths at 144, 108, 72, 54, 36 pounds per acre:');
fprintf(' %2d.  %5.2f   %5.2f   %5.2f   %5.2f   %5.2f\n',xxr');  %'

fprintf('\nFRIEDMAN RANKS are WITHIN-BLOCK ranks, here 1 le R(i,j) le %d\n',kk);
fprintf('Cotton strengths for n=%d replications with within-replication ranks\n',nn);

%% Ranks and tie-groups within each of 22 blocks

rr=zeros(nn,kk);  tg=zeros(nn,kk);
sbuff=zeros(kk,1);  % Column-vector buffer
for i=1:nn;  % For all replications
  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;

%%disp(xxd);
fprintf(' %2d.  %5.2f(%3.1f)  %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 by treatment group:');

%% sum(rr,1) is a row vector of column sums of matrix rr
rsum=sum(rr,1);
ravg=rsum/nn;    % ravg is also a row vector

fprintf('   %2d. %7.2f   %6.2f\n', [ [1:kk]'   rsum'   ravg' ]');


%% Friedman statistic
%% Centered rank averages

ravgc=ravg - (kk+1)/2;

ss = (12*nn/(kk*(kk+1)))*sum(ravgc.*ravgc);

% Alternatively:
%ss=0;
%for j=1:kk;  ss=ss+(ravg(j)-(kk+1)/2)*(ravg(j)-(kk+1)/2);  end;
%ss = (12*nn/(kk*(kk+1)))*ss;

df=kk-1;
disp(' ');
disp('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);

disp(' ');
disp('Simulating P-values for the Friedman test:');
disp('Permuting rank values independently for each block');

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

nsims=10000;

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

% Note that S = A_n * FScore + C_n   for  A_N=12n/(k*(k+1)), C_N=-3n(k+1)
%   for  FScore = Sum(j=1,k) ranksum^2
% Thus simulating  S  is the same as simulating  FScore

Fscore=sum(rsum.*rsum);
% Alternatively:
% Fscore=0;
% for j=1:kk;  Fscore=Fscore+rsum(j)*rsum(j);  end;
disp('Fscore = Sum of Squares over Treatment Groups');
fprintf('  of Treatment-Group Rank Sums: %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;
  %% 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(' ');
disp('PAGE TEST for monotonic alternatives');
disp('L is essentially the sum over blocks');
disp('  of the Spearman correlation coefficient of');
disp('  within-block midranks and the treatment-group number');
disp('Specifically, L = sum(i=1,n) sum(j=1,k) j*R(i,j)');
disp('  = sum(j=1,k) j*sum(i=1,n) R(i,j)');
disp(' ');
disp('Observed score and large-sample approximation:');
ll=0;
for j=1:kk;  ll = ll + j*rsum(j);  end;
% Large-sample normal approximation, ignoring ties
lmean = nn*kk*(kk+1)*(kk+1)/4;
lvar = (nn*(kk-1)*kk*kk*(kk+1)*(kk+1))/144;
lstdev = sqrt(lvar);
lz = (ll-lmean)/lstdev;
Pval=1-normint(lz);
fprintf('  L=%g   Lmean=%g   Lstdev=%g\n', ll,lmean,lstdev);
fprintf('  Z=%g   P=%6.4f (one-sided)\n', lz,Pval);


%% Note continuation line
fprintf(...
 '\nCarrying out Nsims=%d Friedman permutations for Page`s test:\n',nsims);

Lobserv=ll;
fprintf('Observed L = %g\n',Lobserv);
 
% Do nsims permutations for a test of Lobserv
rrp=rr;   % Preserve original ranks
Ngr=0;  jj=1:kk;   % Row vector of integers 1-kk
for ns=1:nsims;
  %% Permute rows of rrp independently
  for i=1:nn;  rrp(i,1:kk) = rrp(i,randperm(kk));  end;
  %% Calculate Rll from the column (treatment-group) sums of rrp
  Rll = sum(jj.*sum(rrp,1));
%  Rll=0;
%  for j=1:kk;  sum=0;
%    for i=1:nn;  sum=sum+rrp(i,j);  end;
%    Rll=Rll+j*sum;  end;
  if Rll >= Lobserv  Ngr=Ngr+1;  end;
  end;

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

phat=Ngr/nsims;  pwid=1.960*sqrt(phat*(1-phat)/nsims);
disp('95% CI for true one-sided P-value bracketing estimate of true pvalue:');
disp(' ');
disp([phat-pwid  phat  phat+pwid]);



