
%  Balanced Incomplete Block Design
%
%  Data  X(i,j)  for 1 le i le nn (blocks)  1 le j le kk (treatments)
%    with missing data such that
%
%  (i) Each block has exactly ss observations (1 le ss le kk)
%  (ii) Each treatment group has exactly pp observations (1 le pp le nn)
%  Thus there are  nn*ss = kk*pp  total observations, not nn*kk
%  (iii) Also need:
%     Each pair of treatment groups has lam observations in common
%
%  Assign Friedman-like ranks to the observed data
%    by assigning ranks R(i,j) (1 le R(i,j) le ss) to non-missing
%    values within each block.
%
%  Assumptions (i,ii,iii) are 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('Balanced Incomplete Block Designs');
disp('Data  X(i,j)  J(i,j)=0,1  where');
disp('   Observation X(i,j) is present if J(i,j)=1, otherwise missing');
disp('   1 le i le nn   Subjects or Blocks');
disp('   1 le j le kk   Treatments or Treatment Groups');
disp('   ss observations in each block (Sum(j=1,k) J(i,j)=ss, all i)');
disp('   pp observations for each treatment (Sum(i=1,n) J(i,j)=pp, all j)');
disp('   lam observations common between every pair of treatments');

disp(' ');
disp('Example:');
disp('Percent conversion of Methyl Glucoside to Monovinyl Isomers');
disp('  in presence of acetylene and a strong base under high pressure');
disp('See Table 7.14 (p316) in text');

%% -1 means missing data
%% First column is block number
xxr = [ 1   16   18   -1   32   -1 ; 
        2   19   -1   -1   46   45 ; 
        3   -1   26   39   -1   61 ; 
        4   -1   -1   21   35   55 ; 
        5   -1   19   -1   47   48 ; 
        6   20   -1   33   31   -1 ; 
        7   13   13   34   -1   -1 ; 
        8   21   -1   30   -1   52 ; 
        9   24   10   -1   -1   50 ; 
       10   -1   24   31   37   -1  ];

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

% Display data
disp(' ');
disp('Data:');
disp('Run    250psi   325psi   400psi   475psi   550psi');
for i=1:nn;
  fprintf(' %2d. ',i)
  for j=1:kk;
    if xx(i,j)>0  fprintf('%6d   ', xx(i,j));
    else fprintf('         ');  end;  end;  fprintf('\n');
  end;

fprintf('\nk=%d treatments and n=%d blocks (runs)\n', kk,nn);

disp(' ');
%% Find how many observations in each block or treatment group
numb=zeros(1,nn);  numtr=zeros(1,kk);
for i=1:nn;  for j=1:kk;
  if xx(i,j)>0  numb(i)=numb(i)+1;  numtr(j)=numtr(j)+1;  end;
  end;  end;
%% Find ss=number of observations per block and check
ss=numb(1);
for i=2:nn;
  if numb(i)~=ss 
    fprintf ('ERROR!  First block: ss=%d observations  Block #%d: %d\n',...
       ss,i,numb(i));
    fprintf ('Number of observations per block is not constant !!!\n');
    return;  end;  end;
fprintf ('ss=%d observations in each block\n',ss);
%% Find pp=number of observations per treatment and check
pp=numtr(1);
for j=2:kk;
  if numtr(j)~=pp
    fprintf ('ERROR!  First group: pp=%d observations  Group #%d: %d\n',...
       pp,j,numtr(j));
    fprintf ('Number of observations per treatment group is not constant !!!\n');
    return;  end;  end;
fprintf ('pp=%d observations in each treatment group\n',pp);
%% Find lam=number of observations in common for each pair of treatments
larr=zeros(kk,kk);
lam=-1;
for j1=1:kk;  for j2=j1+1:kk;
  ll=0;  for i=1:nn;  if xx(i,j1)>0 && xx(i,j2)>0  ll=ll+1;  end;  end;
  if lam<0 lam=ll;
  elseif lam~=ll
    fprintf ('ERROR!  First two groups: lam=%d common observations\n',lam);
    fprintf ('  Groups #%d and #%d:  ll=%d common observations\n', j1,j2,ll);
    fprintf ('Number of observations per group pair is not constant !!!\n');
    return;
  end;  end;  end;
fprintf ('lam=%d observations in each pair of treatment groups\n',lam);

mm=nn*ss;
fprintf ('Total number of observations: %dx%d = %d*%d = %d\n',...
  nn,ss, kk,pp, mm);
m1=pp*(ss-1);  m2=lam*(kk-1);
disp ('These two numbers should be the same:');
fprintf ('  pp*(ss-1)=%d  lam*(kk-1)=%d\n', m1,m2);
if m1~=m2
  fprintf('  THESE NUMBERS MUST BE THE SAME !!!\n');  return;  end;


% Assign ranks 1:ss for observations within each block
rr=zeros(nn,kk);  xvals=zeros(ss,1);  rvals=zeros(ss,1);
for i=1:nn;  r=0;
  for j=1:kk;  if xx(i,j)>0 r=r+1;  xvals(r)=xx(i,j);  end;  end;
  [rvals tiegr] = getmidranks(xvals);  r=0;
  for j=1:kk;  if xx(i,j)>0 r=r+1;  rr(i,j)=rvals(r);  end;  end;
  end;

% Display ranks
disp(' ');
disp('Within-block ranks at observed values:');
disp('Run    250psi   325psi   400psi   475psi   550psi');
for i=1:nn;
  fprintf(' %2d. ',i)
  for j=1:kk;
    if xx(i,j)>0  fprintf('%6.1f   ', rr(i,j));
    else fprintf('         ');  end;  end;  fprintf('\n');
  end;
disp('Observed treatment-group rank sums:');
ranksum=sum(rr,1);   % ranksum is a 1xkk column vector
fprintf('     ');
for j=1:kk;  fprintf('%6.1f   ',ranksum(j));  end;  fprintf('\n');

% Lscore = sum(j=1,kk) ranksum(j)^2  for permutations
Lscore = sum(ranksum.*ranksum);
fprintf('\nObserved rank-sum score for permutations:  Lscore=%g\n',Lscore);

% Durbin etal D statistic
rankc = ranksum - pp*(ss+1)/2;  % Mean of each ranksum(j) under permutations
dd = (12/(lam*kk*(ss+1)))*sum(rankc.*rankc);
df=kk-1;  Pval=pvchisq(dd,df);
disp('Large-sample approximation:');
fprintf('Durbin etal test:  D=%g   P=%7.5f  (df=%d)\n', dd,Pval,df);

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=rr;   % Preserve original rank matrix
Ngr=0;
for ns=1:nsims;
  %% Permute rows of rrp independently
  for i=1:nn;  r=0;
    for j=1:kk;  if xx(i,j)>0  r=r+1;  rvals(r)=rrp(i,j);  end;  end;
    r2vals=rvals(randperm(ss));  r=0;
    for j=1:kk;  if xx(i,j)>0  r=r+1;  rrp(i,j)=r2vals(r);  end;  end;
  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 >= Lscore  Ngr=Ngr+1;  end;
  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);




