/* Test for Balanced Incomplete Block design */
/* See Section 7.6 in the text (Hollander & Wolfe, 2nd edn)

   A two-way layout with k treatment groups and n blocks
     is an INCOMPLETE BLOCK DESIGN (IBD) if it has cells with no data
   It is a BALANCED INCOMPLETE BLOCK DESIGN (BIBD) if there exist
     constants s, p, and lambda such that
      (i) each block has values for s treatments
     (ii) each treatment group has values for p blocks
    (iii) each PAIR of treatments shares exactly lamba values
     (iv) p*(s-1) = lambda*(k-1)

   These conditions are necessary for the large-sample chi-square
     approximation. They are not necessary for the permutation test,
     but it may be difficult to find a reasonable score for a
     permutation test for an extremely unbalanced design.
         */

#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 20       /* Allow for up to this many treatments */
#define MAXN 20       /* Allow for up to this many blocks */

/* Text for this example */
/* Note \n characters at the end of each line except for the last line */
/* The `const' says that the program should not change these values */
/*   without a compiler crash or at least a very strong warning. */
const char *rat_title_data =
 "Reactions of rats to chemical substances on their skin\n"
 "Text (Hollander & Wolfe, 2nd ed) Table 7.13, p315\n"
 "Data Source: Mendenhall (1968)";

/* Number of treatment groups and subjects */
const int rat_kk=7, rat_nn=7;

const char *rat_treatname[MAXK] =
 { "A", "B", "C", "D", "E", "F", "G" };

const char *rat_blockprefix = "Rat";

/* See OneWayLayout.c for comments about DOUBLE ARRAYS in C
     and how they are initialized.

   The program checks BIBD assumptions as a check for data entry.
   Nonobserved values (empty cells) are entered as 0
   If 0 (or negative values) could be legitimate observations, you
     could use a second array with integer values:

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

   with 1s for observed values and 0s for omitted observations.
   
   In the following, 0s in a data array indicate missing data. */

const double rat_xx[MAXN][MAXK] =
 { { 10.2,   6.9,     0,  14.2,    0,     0,     0  },
   {    0,     0,   9.9,  12.9,    0,  14.1,     0  },
   {    0,  12.1,  11.7,     0,  8.6,     0,     0  },
   {    0,     0,     0,  14.3,  9.1,     0,   7.7  },
   {    0,   8.8,     0,     0,    0,  16.3,   8.6  },
   { 13.1,     0,     0,     0,  9.2,  15.2,     0  },
   { 11.3,     0,   9.7,     0,     0,    0,   6.2  } };

/* printf() format for displaying data within 9 horizontal spaces */
const char *rat_format = "  %4.1f   ";

/* A second data set : */

char *aphid_title_data=
 "Logarithms of toxic dosages of seven chemicals to kill 95% of aphids\n"
 "Text (Hollander & Wolfe, 2nd ed) Table 7.12, p311\n"
 "Moore and Bliss (1942)";

/* Number of treatment groups and subjects */
int aphid_kk=7, aphid_nn=7;

const char *aphid_treatname[MAXK] =
 { "A", "B", "C", "D", "E", "F", "G" };

const char *aphid_blockprefix = "Run";

const double aphid_xx[MAXN][MAXK] =
 { { 0.465, 0.343,    0,  0.396,   0,     0,      0  },
   { 0.602,    0,  0.873,    0,  0.634,   0,      0  },
   {    0,     0,  0.875, 0.325,   0,     0,   0.330 },
   { 0.423,    0,     0,     0,    0,   0.987, 0.426 },
   {    0,  0.652, 1.142,    0,    0,   0.989,    0  },
   {    0,  0.536,    0,     0,  0.409,   0,   0.309 },
   {    0,     0,     0,  0.609, 0.417, 0.931,    0  } };

/* printf() format for displaying data within 9 horizontal spaces */
const char *aphid_format = "  %5.3f  ";




/* Number of simulations and default starting random-number seed */
unsigned long wseed=123456;
int nsims=100000;

/* Prototype of function to display and analyze a BIBD */
void do_twoway_bibd(const double xx[MAXN][MAXK], int kk, int nn,
   const char *title_data, const char *treatname[MAXK], 
   const char *blockprefix, const char *dataformat);

int main(int argc, char **argv)
 { /* Helpful message on segmentation fault and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();

   printf ("A two-way layout with missing values:\n");
   printf ("  Balanced Incomplete Block Designs (BIBD)\n");

   /* Enter `ProgramName <Number>' to customize 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");

   /* Display and analyze the first example BIBD block design */
   printf ("\nFIRST DATASET:\n");
   do_twoway_bibd(rat_xx, rat_kk, rat_nn, rat_title_data, 
     rat_treatname, rat_blockprefix, rat_format);

   /* Display and analyze the first example BIBD block design */
   printf ("\nSECOND DATASET:\n");
   do_twoway_bibd(aphid_xx, aphid_kk, aphid_nn, aphid_title_data, 
     aphid_treatname, aphid_blockprefix, aphid_format);

   show_randnums_used();

   return 0; }



/* Routines used by do_twoway_bibd() */

void show_twoway_bibd_data(const double xx[MAXN][MAXK], int kk, int nn,
   const char *treatname[MAXK], const char *blockprefix,
   const char *dataformat);
void show_twoway_bibd_ranks(double mr[MAXN][MAXK], int kk, int nn,
   const char *treatname[MAXK], const char *blockprefix);

/* Display and analyze a BIBD block design */

void do_twoway_bibd(const double xx[MAXN][MAXK], int kk, int nn,
   const char *title_data, const char *treatname[MAXK],
   const char *blockprefix, const char *dataformat)
 { int i,j,r, nbig,ns;
   int ss,rr,la, n1,n2;
   double dd, avdd, dscore, pp,pwid;
   /* Buffers for finding within-block ranks */
   double xvals[MAXK], mranks[MAXK], ranksum[MAXK];
   /* Array of within-block midranks */
   double midr[MAXN][MAXK];

   /* Description of data */
   printf ("\n%s:\n", title_data);
   printf ("  (%d treatments, %d subjects)\n", kk,nn);

   /* Display the data */
   printf ("\nData:\n");
   show_twoway_bibd_data(xx,kk,nn, treatname, blockprefix, dataformat);

   /* Find within-block ranks */
   for (i=0; i<nn; i++)  /* Sum over blocks */
    { /* Copy observed values in i-th block to xvals[] */
      for (r=j=0; j<kk; j++) if (xx[i][j]>0.0)
        xvals[r++]=xx[i][j];
      /* Put midranks in mvals[] */
      /* This functions is in statlib.c */
      /* Don't collect tie-group information */
      makeranks(mranks, xvals,r, NULL,NULL);
      /* Copy the midranks for i-th block to midr[][] */
      /* Make sure midr[i][j]=0.0 for missing values */
      for (r=j=0; j<kk; j++)
       { if (xx[i][j]>0.0)  midr[i][j]=mranks[r++];
         else midr[i][j]=0.0; }}

   /* Display midranks */
   printf ("\n\nMidranks:\n");
   show_twoway_bibd_ranks(midr,kk,nn, treatname, blockprefix);

   /* Check to see if all blocks have the same number (ss) of observations */
   /*    all treatment groups have the same number (rr) of observations */
   /* We use rr for what is p in the text to preserve pp for P-value */
   ss=rr=-1;  /* Starting values */
   for (i=0; i<nn; i++)
    { int nvalues=0;   /* Number of data values in this block */ 
      /* Count the number of entries in the i-th row */ 
      for (j=0; j<kk; j++) if (xx[i][j]>0.0) nvalues++;
      if (ss<0) ss=nvalues;  /* First block */
      else if (ss!=nvalues)
       { printf ("\nERROR!  Number of values in block#1: %d  Block#%d: %d\n"
           "  Must be the same!\n", ss, i+1,nvalues);  exit(1); }}
   /* Check to make sure that all treatments have the same number of values */
   /* Calculate ranks sums (for each treatment group) at the same time */
   for (j=0; j<kk; j++)
    { int nvalues=0;  double sum=0.0;
      /* Count the number of entries in the j-th column */
      /* Also such the ranks in that column */ 
      for (i=0; i<nn; i++)
        if (midr[i][j]>0.0) { nvalues++;  sum+=midr[i][j]; }
      if (rr<0) rr=nvalues;
      else if (rr!=nvalues)
       { printf ("\nERROR!  Values in Trt#1: %d  Trt#%d: %d\n"
           "  Must be the same!\n", rr, j+1,nvalues);  exit(1); }
      ranksum[j]=sum; }
   /* Lets check the other two conditions */
   la=-1;  /* For every pair of treatments */
   for (i=0; i<kk; i++) for (j=i+1; j<kk; j++)
    { int nvalues=0; /* For this pair of treatment groups */
      for (r=0; r<nn; r++)
        if (xx[r][i]>0.0 && xx[r][j]>0.0)  nvalues++;
      if (la<0) la=nvalues;
      else if (la!=nvalues)
       { printf ("\nERROR!  Values in Trts#1,#2: %d  Trts#%d,#%d: %d\n"
           "  Must be the same!\n", la, i+1,j+1, nvalues);  exit(1); }}

   printf ("\nData is a BIBD design with\n"
     "  k=%d, s=%d,  n=%d, p=%d,  la=%d, ", kk,ss,nn,rr,la);
   n1=rr*(ss-1);  n2=la*(kk-1);
   if (n1!=n2)
    { printf ("\nERROR!  But with  p*(s-1)=%d != la*(kk-1)=%d\n", 
         n1,n2);  exit(1); }
   else
     printf (" p*(s-1) = la*(kk-1) = %d\n", n1);

   /* Compute rank-sum score and permutation score: */
   dd=0.0;  avdd=(double)(rr*(ss+1))/2;
   for (j=0; j<kk; j++)  dd+=(ranksum[j]-avdd)*(ranksum[j]-avdd);
   dd*=12;  dd/=(la*kk*(ss+1));
   printf ("\nDurbin-Skillings-Mack statistic  D=%g\n", dd);
   printf ("  Large-sample chi-square approx:  P=%7.5f  df=%d\n",
     pvchisq(kk-1,dd), kk-1);

   /* NOTE: The function `commnum(nn)' in statlib.c allows you to display
       a number with commas, as in
           nn:        100   1000   10000    100000    1000000     10000000
        commnum(nn):  100   1000   10,000   100,000   1,000,000   10,000,000
      WARNING:  commnum()  returns a pointer to a character string,
        so should be `%s' in a printf() function, not `%d'  */

   printf ("\nCarrying out %s permutations to estimate\n"
     "  the exact P-value for D=%g:\n", commnum(nsims), dd);
   dscore=0.0;
   for (j=0; j<kk; j++)  dscore+=ranksum[j]*ranksum[j];
   printf ("Observed permutation score = Sum of squares of\n"
     "  treatment-group rank sums = %g\n", dscore); 
   for (nbig=ns=0; ns<nsims; ns++)
    { double rscore;
      /* Permute the ranks of observed values within each block: */
      for (i=0; i<nn; i++)
       { double mvals[MAXK];
         /* Copy midranks for block #i to mvals[], shuffle, and return them */
         for (r=j=0; j<kk; j++) if (midr[i][j]>0.0)
           mvals[r++]=midr[i][j];
         shuffle_dub(mvals,r);
         for (r=j=0; j<kk; j++) if (midr[i][j]>0.0)
           midr[i][j] = mvals[r++]; }
      /* Recompute the Dscore permutation statistic */
      rscore=0.0;
      for (j=0; j<kk; j++)
       { double rsum=0;
         for (i=0; i<nn; i++) if (midr[i][j]>0.0) rsum+=midr[i][j];
         rscore+=rsum*rsum; }
      /* Count how many times rscore is bigger than the observed? */
      if (rscore>=dscore) 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/%s = %7.5f  95%% CI (%7.5f, %7.5f, %7.5f)\n",
     nbig,commnum(nsims), pp, pp-pwid,pp,pp+pwid);
       }


void show_twoway_bibd_data(const double xx[MAXN][MAXK], int kk, int nn,
   const char *treatname[MAXK], const char *blockprefix,
   const char *dataformat)
 { char nbuff[388];  int i,j;
   for (i=0; i<kk*9; i++) nbuff[i]='-';  nbuff[i]=0;
   /* Column headings are treatment names */
   printf ("     ");
   for (j=0; j<kk; j++) printf ("     %-4s", treatname[j]);  printf("\n");
   printf ("-----%s\n", nbuff);
   /* Display data with blanks for missing values */
   for (i=0; i<nn; i++)
    { printf ("%s#%d ", blockprefix,i+1);
      for (j=0; j<kk; j++)
       { if (xx[i][j]<=0.0) printf ("         ");
         else  printf (dataformat, xx[i][j]); }
      printf("\n"); }
   /* Treatment averages */
   printf ("-----%s\n", nbuff);
   printf ("Avg:  ");
   for (j=0; j<kk; j++)
    { double sum=0.0;  int n=0;
      for (i=0; i<nn; i++) if (xx[i][j]>0.0) { sum+=xx[i][j];  n++; }
      printf (dataformat, sum/n); }
   printf ("\n");
      }


void show_twoway_bibd_ranks(double midr[MAXN][MAXK], int kk, int nn,
   const char *treatname[MAXK], const char *blockprefix)
 { char nbuff[388];  int i,j;
   for (i=0; i<kk*9; i++) nbuff[i]='-';  nbuff[i]=0;
   /* Column headings are treatment names */
   printf ("     ");
   for (j=0; j<kk; j++) printf ("     %-4s", treatname[j]);  printf("\n");
   printf ("-----%s\n", nbuff);
   /* Display data with blanks for missing values */
   for (i=0; i<nn; i++)
    { printf ("%s#%d", blockprefix,i+1);
      for (j=0; j<kk; j++)
       { double mr=midr[i][j];
         if (mr<=0.0) printf ("         ");
         else if (fabs(mr-floor(mr))<0.10)
             printf ("   %3g   ", mr);
         else  printf ("  %3.1f  ", mr); }
      printf("\n"); }
   printf ("-----%s\n", nbuff);
   printf ("Sum: ");
   for (j=0; j<kk; j++)
    { double sum=0.0;
      for (i=0; i<nn; i++) if (midr[i][j]>0.0) sum+=midr[i][j];
      if (fabs(sum-floor(sum))<0.10) printf ("   %3g   ", sum);
      else printf ("   %3.1f   ", sum); }
   printf ("\n");
      }