/* Cochran test for twoway layouts with 0,1 (dichotomous) data

   Suppose that we are given data  X_ij  where each X_ij is 1 or 0
     (`success' or `failure') for n blocks or subjects (1 le i le n)
     and k treatments (1 le j le k). Suppose 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.

   Let L_j = Sum(i=n) X_ij be the number of successes for the j-th
     treatment.  Since

         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 = (1/k)Sum(j=1,k) Sum(i=1,n) X_ij
     = M/k   where M is the total number of successes, a natural score
     for different effects of treatments is

         T = 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 / (kM - Sum(i=1,n) B_i^2)

     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.  In the formula for Q,
     Q, B_i = Sum(j=1,k) X_ij  is the number of successes in the i-th
     block, M = Sum (i=1,n) B_i = Sum (j=1,k) L_j is the total number
     of successes, and Lbar = Sum(j=1,k) L_j/k = M/k  as above.

     Curiously, this test is exactly equivalent to Friedman's test except
     that T and Q are easier to compute 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  (X_ij=0)  = (1+Z_i+k)/2  (X_ij=1)

     where  Z_i = k-B_i  is the number of failures (or 0s) 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. Then the Friedman
     rank sums

         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=M/k

      so that

         (R_+j - Rbar) = (k/2) (L_j - Lbar)

     Similarly, the usual large-sample approximation with tie correction
     for the Friedman statistic S' is exactly equal to Q. The denominator
     in Q results from the Friedman-test tie correction and a little bit
     of algebra. */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Prototypes for functions in  statlib.c
   Compile by e.g.   gcc -o ProgName Progname.c statlib.c   with
     ProgName.c, statlib.c, and statlib.h in the default directory. */
#include "statlib.h"

#define MAXK 10       /* Maximum possible number of treatments */
#define MAXN 50       /* Maximum possible number of subjects/blocks */

char *title_data =
 "Success for four methods of soothing newborn babies\n"
 "Source: Lehmann ``Nonparametrics: Statistical methods based on ranks''\n"
 "  pages 267-270\n"
 "See Lehmann p283 for the full data set.";

/* The number of treatment groups and the number of blocks (subjects) */
const int kk=4, nn=12;

/* The default starting random-number seed */
unsigned long wseed=408408;
/* Number of simulations to estimate P-values */
int nsims=100000;

/* See OneWayLayout.c for comments about ARRAYS in C
     and how they are initialized.
   As in TwoWayFriedman.c, we store treatment groups as columns
     and blocks (subjects) as rows
   Note that we use an integer array here, not a double array */

/* Treatment group name */
const char *treatname[MAXK] = 
  { "WarmWater", "Rocking", "Pacifier", "Sound" };
/* Shorter column heading names */
const char *treatname_short[MAXK] =
  { "WrmWa", "Rockg", "Pacif", "Sound" };

const int xx[MAXN][MAXK] = {
  { 0,  0,  0,  0 },
  { 0,  1,  0,  0 },
  { 0,  0,  1,  0 },
  { 0,  0,  0,  0 },
  { 1,  1,  1,  1 },
  { 0,  1,  1,  1 },
  { 1,  1,  0,  0 },
  { 0,  1,  0,  1 },
  { 0,  1,  1,  1 },
  { 0,  0,  1,  0 },
  { 1,  1,  0,  1 },
  { 1,  1,  1,  1 } };


/* WARNING: Data here are in an integer array, NOT a double array */
/*   as in other handouts */

void show_cochran_data(const int xx[MAXN][MAXK],
  int kk, int nn, const char *treatname[MAXK]);
void do_cochran(const int xx[MAXN][MAXK], int kk, int nn);

void show_friedman_data(const int xx[MAXN][MAXK], 
  int kk, int nn, const char *treatname[MAXK]);
void do_friedman(const int xx[MAXN][MAXK], int kk, int nn);

int main(int argc, char **argv)
 { int i;
   printf ("\nIllustration of Cochran test for two-way layouts\n"
       "  with 0,1 (dichotomous) data\n");
   /* Helpful messages on segmentation faults and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();

   /* Text for this example */
   printf ("\n%s\n", title_data);
   printf ("k=%d treatment groups and n=%d subjects:\n", kk,nn);

   /* Enter `ProgramName <Number>' to set starting RN seed=<Number> */
   if (argc>=2) wseed=atol(argv[1]);
   wseed=setrand(wseed);  /* If wseed=0, replace by true starting seed */
   printf ("\nStarting random-number seed: %lu\n",wseed);
   printf ("(Enter a number on the command line "
     "to specify a different seed.)\n");

   printf ("\nTable of treatment group headings:\n");
   for (i=0; i<kk; i++)  printf("  %5s - %s\n", treatname_short[i], treatname[i]);
   printf("\n");
   show_cochran_data(xx,kk,nn,treatname_short);
   /* Carry out the Cochran test */
   do_cochran(xx, kk,nn);

   printf ("\nApply the Friedman procedure as a comparison:\n");
   show_friedman_data(xx, kk,nn, treatname);
   /* Carry out the Friedman test */
   do_friedman(xx, kk,nn);

   /* Show how many random numbers were used */
   show_randnums_used();
   /* No errors to report */
   return 0; }


void show_cochran_data(const int xx[MAXN][MAXK],
  int kk, int nn, const char *treatname[MAXK])
 { int i,j,msum;
   printf ("Dichotomous data in a two-way layout:\n");
   printf ("     ");
   for (j=0; j<kk; j++) printf ("%5s  ", treatname[j]);
   /* The last column heading */
   printf ("Block sum\n");
   msum=0;
   for (i=0; i<nn; i++)
    { int bb;
      printf (" %2d. ", i+1);
      for (j=0; j<kk; j++) printf ("  %d    ", xx[i][j]);
      bb=0; for (j=0; j<kk; j++) bb+=xx[i][j];
      printf ("  %2d\n",bb);
      msum+=bb; } /* Keep track of total sum */
   printf ("Treatment group sums:\n");
   printf ("     ");
   for (j=0; j<kk; j++)
    { int ll=0;
      for (i=0; i<nn; i++) ll+=xx[i][j];
      printf (" %2d    ", ll); }
   printf ("  %2d\n",msum); }


void do_cochran(const int xx[MAXN][MAXK], int kk, int nn)
 { int i,j, ns,nbig;
   int bb[MAXN], ll[MAXK], xtest[MAXN][MAXK], msum, l2sum,b2sum;
   double lbar,cc, pp,pwid;
   printf ("\nApplying the Cochran test:\n");
   /* Fill bb[] (block sums), ll[] (treatment sums), msum (total sum) 
     l2sum = Sum(j,k) ll[j]*ll[j] is the observed Cochran statistic
        for permutations
     b2sum = Sum(i,n) bb[i]*bb[i] is needed for tie correction */
   msum=l2sum=b2sum=0;
   for (j=0; j<kk; j++) 
    { int lsum=0;
      /* Treatment-group sum of successes */
      for (i=0; i<nn; i++) lsum+=xx[i][j];
      ll[j]=lsum;  msum+=lsum;  l2sum+=lsum*lsum; };
   /* Repeat displays as a check: */
   printf ("Treatment sums:");
   for (j=0; j<kk; j++) printf ("  %d",ll[j]);
   printf("  Total: %d\n",msum);
   /* Now the block sums */
   for (i=0; i<nn; i++)
    { int bsum=0;  for (j=0; j<kk; j++) bsum+=xx[i][j];
      bb[i]=bsum;  b2sum += bb[i]*bb[i]; }
   printf ("Block sums:");
   for (i=0; i<nn; i++) printf ("  %d",bb[i]);  printf("\n");
   /* Calculate the Cochran statistic cc: */
   lbar=(double)msum/kk;
   cc=0;  for (j=0; j<kk; j++) cc+=(ll[j]-lbar)*(ll[j]-lbar);
   cc = (kk*(kk-1)*cc)/(kk*msum-b2sum);
   printf ("Cochran statistics:   L2sum=%d   CC=%g\n", l2sum,cc);
   printf ("Large sample approximation:   C=%7.5f   P=%6.4f   df=%d\n",
     cc, pvchisq(kk-1,cc), kk-1);
   printf ("Permutation test for L2sum=%d\n"
      "  using %d sets of within-block permutations:\n", l2sum, nsims);
   /* Make a copy of the data for permutations */
   /*   (since xx[][] is protected by `const int') */
   for (i=0; i<nn; i++) for (j=0; j<kk; j++)
      xtest[i][j]=xx[i][j];
   /* Initialize permutation success count */
   nbig=0;
   for (ns=0; ns<nsims; ns++)
    { int testsum;
      /* Permute all rows */
      for (i=0; i<nn; i++) shuffle(xtest[i],kk);
      /* Calculate l2sum statistic for shuffled data */
      testsum=0;
      for (j=0; j<kk; j++)
       { int ll=0;
         for (i=0; i<nn; i++) ll+=xtest[i][j];
         testsum += ll*ll; }
      /* Success or failure? */
      if (testsum >= l2sum) nbig++; }

   /* Display the simulated P-value with a 95% confidence interval */
   pp=(double)nbig/nsims;
   pwid=1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P = %d/%d = %6.4f   95%% CI (%5.3f, %5.3f, %5.3f)\n",
     nbig, nsims, pp, pp-pwid, pp, pp+pwid);
           }



void show_friedman_data(const int xx[MAXN][MAXK],
  int kk, int nn, const char *treatname[MAXK])
 { char nbuff[88];
   int i,j;
   double midr[MAXN][MAXK];
   /* Find Friedman (within block) ranks (so 1 <= R(i,j) <= k) */
   printf ("\nData with Friedman (within-block) ranks:\n");
   for (i=0; i<nn; i++)  /* For each block */
    { double xdub[MAXK];
      /* WARNING: Midrank routine needs data as doubles, not integers */
      for (j=0; j<kk; j++)  xdub[j]=xx[i][j];
      /* We don't need tiegroup information, so tell makeranks()
           not to return it */
      makeranks(midr[i], xdub,kk, NULL,NULL); }
   /* Use treatnames[] to give column (treatment) headings */
   /* Skip if user says no treatnames by entering treatname=NULL */
   if (treatname!=NULL)
    { printf ("        ");
      for (j=0; j<kk; j++)
        printf ("%-14s", treatname[j]);  printf ("\n"); }
   /* Display data and within-block midranks */
   for (i=0; i<nn; i++)  /* For each block (subject) */
    { printf (" %2d.    ", i+1);
      for (j=0; j<kk; j++)  /* For each treatment group (column) */
       { /* This is necessary to make all columns the same width */
         sprintf(nbuff, "   %d  (%g)", xx[i][j], midr[i][j]);
         printf ("%-14s", nbuff); }
      printf("\n"); }
   printf ("AvRanks:");
   for (j=0; j<kk; j++)
    { double rowsum=0.0;
      /* Compute the average of the midranks in column #j */
      for (i=0; i<nn; i++) rowsum+=midr[i][j];
      sprintf (nbuff, "  (%.2f)", rowsum/nn);
      printf("%-14s", nbuff); }
   printf ("\n"); }



/* Adapted from TwoWayFriedman.c
   See TwoWayFriedman.c for comments */

void do_friedman(const int xx[MAXN][MAXK], int kk, int nn)
 { int i,j, ns,nbig, totntg;
   /* Midranks of data in xx[MAXN][MAXK] */
   double midr[MAXN][MAXK];
   /* Miscellaneous doubles */
   double tiesum, obscore,ss0,sstc, pp,pwid;

   /* Describe the data */
   printf ("\nCarrying out classical Friedman test:\n");
   printf ("%d treatments, %d subjects\n", kk,nn);
   /* Find Friedman (within block) ranks (so 1 <= R(i,j) <= k)
       and sum of within-block tiesums */
   tiesum=0.0;  totntg=0;
   for (i=0; i<nn; i++)  /* For each block */
    { double xdub[MAXK];
      int ntg, tiegroup[MAXK];  double tsum;
      /* WARNING: Midrank routine needs data as doubles, not integers */
      for (j=0; j<kk; j++)  xdub[j]=xx[i][j];
      makeranks(midr[i], xdub,kk, tiegroup,&ntg);
      /* Add to cumulative tiesum */
      tsum=0.0;
      for (j=0; j<ntg; j++) tsum += (tiegroup[j]-1)*tiegroup[j]*(tiegroup[j]+1);
      tiesum+=tsum;  totntg+=ntg; }

   /* The Friedman rank test: */
   /* The sum of squares of the treatment rank sums */
   obscore=0.0;
   for (j=0; j<kk; j++)
    { double rsum=0.0;
      for (i=0; i<nn; i++) rsum+= midr[i][j];
      obscore += rsum*rsum; }
   /* The Friedman statistic S (scaled to be asymptotically chi-square) */
   sstc = ss0 = ((12*obscore)/(nn*kk*(kk+1))) - 3*nn*(kk+1);
   /* Is a tie correction needed? */
   if (totntg>0)
     sstc /= (1.0 - tiesum/((double)nn*(kk-1)*kk*(kk+1)));
   printf ("\nFriedman statistic S'=%7.5f   P=%6.4f (%d df, %d tiegroup%s)\n",
     sstc, pvchisq(kk-1,sstc), kk-1, totntg, (totntg==1)?"":"s");
   if (totntg>0)
    { printf ("WITH NO TIE CORRECTION:\n");
      printf ("Friedman statistic S=%7.5f   P=%6.4f (%d df)\n",
        ss0, pvchisq(kk-1,ss0), kk-1); }
   /* Do simulation for Friedman statistic */
   printf ("\nSimulating the Friedman P-value using %d permutations:\n",
        nsims);
   printf ("  Friedman score (for simulations):  %g\n", obscore);
   /* Carry out the simulations using Friedman permutations */
   nbig=0;
   for (ns=0; ns<nsims; ns++)
    { double newscore;
      /* Shuffle the midrank values for each subject */
      /* Recall that midr[i] points to the start of the i-th row */
      for (i=0; i<nn; i++)  shuffle_dub(midr[i],kk);
      /* Recompute the rank-sum score */
      for (newscore=0.0, j=0; j<kk; j++)
       { double rsum=0.0;
         for (i=0; i<nn; i++)  rsum+=midr[i][j];
         newscore += rsum*rsum; }
      /* Newscore= int/4, so no floating-point roundoff error */
      if (newscore>=obscore)  nbig++; }
   /* Display the simulated P-value with a 95% confidence interval */
   pp=(double)nbig/nsims;
   pwid=1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P = %d/%d = %6.4f   95%% CI (%5.3f, %5.3f, %5.3f)\n",
     nbig, nsims, pp, pp-pwid, pp, pp+pwid);
         }