
%   Cochran`s test for two-way layouts with 0,1 (dichotomous) data
%
%   Suppose that we have data  X_ij  where each X_ij is 1 or 0
%     (`success` or `failure`). There are n blocks or subjects (1 le i le n)
%     and k treatments (1 le j le k). 
%     
%   Assume that we want to test if there is a treatment effect.
%     That is, if some treatment groups have more successes than
%     others to a statistically significant extent.
%
%   The row (block) and column (treatment-group) sums are
%
%         L_j = Sum(i=1,n) X_ij
%         B_i = Sum(j=1,k) X_ij
%
%   Define
%
%         L1 = Sum(j=1,k) L_j
%         L2 = Sum(j=1,k) L_j^2
%         B1 = Sum(i=1,n) B_i
%         B2 = Sum(i=1,n) B_i^2
%
%   Then L1=B1 is the total number of successes, and
%   
%         Sum(j=1,k) (L_j-Lbar)^2 = Sum(j=1,k) L_j^2 - k Lbar^2
%
%     where Lbar = (1/k)Sum(j=1,k) L_j = L1/k  is a nature score to measure
%     whether or not there is a treatment effect. This is equivalent to
%
%         L2 = Sum(j=1,k) L_j^2
%
%   A natural permutation structure is to take independent permutations
%     of each block as in the Friedman test. Cochran showed that the
%     statistic
%
%         Q = k(k-1) Sum(j=1,k) (L_j-Lbar)^2 / (kB1 - Sum(i=1,n) B_i^2)
%
%           = k(k-1) (L2 - L1*L1/k) / (k*B1 - B2)
%
%     has a chi-square distribution with k-1 degrees of freedom under
%     Friedman-like permutations for large n, so that Q can be used for
%     a large-sample approximation of the P-value.
%
%   Curiously, this test is exactly equivalent to Friedman`s test, in spite of
%     the fact its statistic is much easier to calculate for dichotomous data.
%     The key is that the within-block midranks of the 0s and 1s in the
%     i-th block are equal to
%
%         R_ij = (1+Z_i)/2  (if X_ij=0)  = (1+Z_i+k)/2  (if X_ij=1)
%
%     where  Z_i = k-B_i  is the number 0s (zeroes) in the i-th block. This is
%     because the 0 scores have tentative ranks 1 to Z_i while the 1 scores have
%     tentative ranks Z_1+1 to k. The Friedman rank sums are then
%
%         R_j = Sum(i=1,n) R_ij = (Sum (i=1,n) (1+Z_i)/2) + (k/2) L_j
%
%             = n(k+1)/2 - (k/2)(L_j - Lbar)   since  Lbar=L1/k
%
%     Thus
%
%         (R_j - Rbar) = (k/2) (L_j - Lbar)
%
%     and
%
%         Sum(j=1,k) (R_j - Rbar)^2 = (k^2/4) Sum(j=1,k) (L_j - L1/k)^2
%
%     A moderate amount of algebra to simplify the tie-correction term in
%     Friedman`s test (see the text) then shows that  S`=Q , where S` is the
%     tie-corrected version of Friedman`s statistic and Q is as above.
%


disp('Cochran`s Test:');
disp('Two-Way Layout with 0,1 data:');
disp('Data  X(i,j)=0,1  where');
disp('   1 le i le nn   Subjects or Blocks');
disp('   1 le j le kk   Treatments or Treatment Groups');
disp('Same model and H_0 as in Friedman`s test');

disp(' ');
disp('Example: Success for four methods of soothing newborn babies');
disp('From: Lehmann ``Nonparametrics: Statistical methods based on ranks``');
disp('  pages 267-270');
disp('See Lehmann p283 for the full data set.');

disp(' ');
disp('k=4 treatments and n=12 newborns');
disp('Treatments:');
disp('  1. Warm Water');
disp('  2. Rocking');
disp('  3. Pacifier');
disp('  4. Sound');

xxr = [ 1     0      0      0      0  ;
        2     0      1      0      0  ;
        3     0      0      1      0  ;
        4     0      0      0      0  ;
        5     1      1      1      1  ;
        6     0      1      1      1  ;
        7     1      1      0      0  ;
        8     0      1      0      1  ;
        9     0      1      1      1  ;
       10     0      0      1      0  ;
       11     1      1      0      1  ;
       12     1      1      1      1   ];

[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('B_i are row sums and L_j are column sums');
disp('Row number, warm water, rocking, pacifier, sound:');

% Add block sums (row sums)
rowsums=zeros(nn,1);  llsums=zeros(1,kk);
for i=1:nn;  for j=1:kk;
   rowsums(i)=rowsums(i)+xx(i,j);
   llsums(j)=llsums(j)+xx(i,j);  end; end;

% Note continuation line
fprintf('  %2d.   %2d   %2d   %2d   %2d    Sum: %d  (B(%2d))\n', ...
  [xxr rowsums xxr(:,1)]'); %'

disp(' ');
disp('Treatment-group (column) sums (L(j)):');
fprintf('        %2d   %2d   %2d   %2d    Sum: %2d\n', ...
  [llsums  sum(llsums)]);

% Cochran statistics
ll1 = sum(llsums);
ll2 = sum(llsums.*llsums);
df = kk-1;
bb1 = sum(rowsums);
bb2 = sum(rowsums.*rowsums);

disp(' ');
disp('Row and column-sum statistics:');
disp('(L1,B1 are sums, L2,B2 are sums of squares)');
disp(' ');
fprintf('   L1=%g   L2=%g   B1=%g   B2=%g\n', ll1,ll2,bb1,bb2);

% Cochran statistic in form for large-sample chi-square approximation
qq = kk*(kk-1)*(ll2 - ll1*ll1/kk)/(kk*bb1 - bb2);

fprintf('\nCochran-test statistics:\n\n');
fprintf('   L2=%g   Q=%5.3f   Pval=%6.4f  (df=%d)\n',...
  ll2, qq, pvchisq(qq,df), df);

disp(' ');
disp('Simulating P-values for the Cochran test:');
disp('Permuting data values (=0,1) 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);

% Do nsims permutations for a test of Fscore
xxp=xx;   % Preserve original data
Ngr=0;
for ns=1:nsims;
  %% Permute rows of xxp independently
  for i=1:nn;  xxp(i,1:kk) = xxp(i,randperm(kk));  end;
  %% Sum of squares of the column sums of xxp
  Rscore=sum(sum(xxp,1).*sum(xxp,1));
  if Rscore >= ll2  Ngr=Ngr+1;  end;
  end;

disp('Number of simulations with values >= L2 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('Applying the Friedman test procedure:');

% Find within-subject ranks and tie groups

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('Finding within-block midranks for each block:');
fprintf(' %2d.  %d(%3.1f)  %d(%3.1f)  %d(%3.1f)  %d(%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);

df=kk-1;
disp(' ');
disp('With no tie correction:');
Pval = pvchisq(ss,df);
fprintf('  S=%5.3f   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=%5.3f   P= %6.4f  (df=%d)\n', ssp,Pval,df);

disp(' ');
disp('Note that the second set of values are exactly');
disp('  the same as in Cochran`s test');


