/* Friedman and Page tests for twoway layouts */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Prototypes for functions in  statlib.c */
#include "statlib.h"

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

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

/* See OneWayLayout.c for comments about DOUBLE ARRAYS in C
     and how they are initialized.
   For twoway layouts with blocks and/or subjects, we store treatment
     groups as columns so that the data for each subject or block is a row.
     Rows in higher-order arrays are stored in memory as one-dimensional
     arrays, which makes creating Friedman-style midranks and carrying
     out within-block permutations easier.

   The `const' in the declarations below says that the program should not
     change these values without a compiler crash or at least a very
     strong warning. This is intended as a safety feature for
     inattentive programmers. */


/* First data set (for the Friedman test): */

/* Strategies of running to second base for 22 baseball players */

/* The number of treatment groups and the number of subjects (blocks) */
const int kkb=3, nnb=22;

/* Names of treatment groups for displays */
const char *baseball_treatnames[MAXK] =
 { "Round Out",
   "Narrow Angle",
   "Wide Angle" };

/* Two dimensional data with
      i-th row are times for i-th subject
      j-th column is j-th treatment group */

const double baseball_dat[MAXN][MAXK] =
 { { 5.40,  5.50, 5.55 }, 
   { 5.85,  5.70, 5.75 }, 
   { 5.20,  5.60, 5.50 }, 
   { 5.55,  5.50, 5.40 }, 
   { 5.90,  5.85, 5.70 }, 
   { 5.45,  5.55, 5.60 }, 
   { 5.40,  5.40, 5.35 }, 
   { 5.45,  5.50, 5.35 }, 
   { 5.25,  5.15, 5.00 }, 
   { 5.85,  5.80, 5.70 }, 
   { 5.25,  5.20, 5.10 }, 
   { 5.65,  5.55, 5.45 }, 
   { 5.60,  5.35, 5.45 }, 
   { 5.05,  5.00, 4.95 }, 
   { 5.50,  5.50, 5.40 }, 
   { 5.45,  5.55, 5.50 }, 
   { 5.55,  5.55, 5.35 }, 
   { 5.45,  5.50, 5.55 }, 
   { 5.50,  5.45, 5.25 }, 
   { 5.65,  5.60, 5.40 }, 
   { 5.70,  5.65, 5.55 }, 
   { 6.30,  6.30, 6.25 } };


/* The second data set (for the Page test): */

/* Strength of cotton fibers as a function of the amount of potash */
/*   in fertilizer */
/* The number of treatment groups and and number of subjects (blocks) */
const int kkp=5, nnp=3;

/* Names of treatment groups for displays */
/* # is a common abbreviation for pounds */
const char *potash_treatnames[MAXK] = 
 { "144#/acre", "108#/acre", "72#/acre", "54#/acre", "36#/acre" };

const double potash_dat[MAXN][MAXK] =
 { { 7.46, 7.17, 7.76, 8.14, 7.63 },
   { 7.68, 7.57, 7.73, 8.15, 8.00 },
   { 7.31, 7.80, 7.74, 7.87, 7.93 } };


/* Prototypes for the two tests */

void do_friedman(const double xx[MAXN][MAXK],
  int kk, int nn, const char *treatnames[MAXK]);
void do_page(const double xx[MAXN][MAXK],
  int kk, int nn, const char *treatnames[MAXK]);

int main(int argc, char **argv)
 { printf ("Illustration of nonparametric tests for two-way layouts\n");
   printf ("Friedman and Page tests.\n");
   /* Helpful messages on segmentation faults and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();
   /* Enter `TwoWayTests <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");

   /* The Friedman Two-Way layout test */
   printf ("\nThe Friedman test for a two-way layout\n");
   /* Text for this example */
   printf ("Time for Rounding First Base under Three Different Strategies\n");
   printf ("Text (Hollander & Wolfe), Table 7.1, p274");
   do_friedman(baseball_dat,kkb,nnb,baseball_treatnames);

   /* The Page Two-Way layout test for increasing alternatives */
   printf ("\nThe Page (increasing alternative) test for a two-way layout\n");
   /* Text for this example */
   printf ("Cotton strength index as a function of potash in fertilizer\n");
   printf ("Text (Hollander & Wolfe), Table 7.5, p286");
   do_page(potash_dat,kkp,nnp,potash_treatnames);

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


/* Display two-way layout data with Friedman (within-block) ranks */

void showranks (const char *treatname[MAXK],
  const double xx[MAXN][MAXK], double midr[MAXN][MAXK], int kk, int nn)
 { char nbuff[88];
   int i,j; 
   /* 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"); }
   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) */
       { /* `sprintf' writes output to a char buffer, not to screen */
         /* This is necessary to make all columns the same width */
         sprintf(nbuff, "%4.2f (%g)", xx[i][j], midr[i][j]);
         printf ("%-14s", nbuff); }
      printf("\n"); }
   printf ("\nAvData: ");
   for (j=0; j<kk; j++)
    { double rowsum=0.0;
      /* Compute the average of column #j */
      for (i=0; i<nn; i++) rowsum+=xx[i][j];
      sprintf (nbuff, "   %.2f", rowsum/nn);
      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"); }



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

   /* Describe the data */
   printf ("%d treatments, %d subjects\n", kk,nn);
   /* First, find Friedman ranks, which are 1<=R<=kk within each block */
   /* NOTE:  If xx[MAXN][MAXK]  is a double array,  */
   /*   then  xx[i]  is a pointer to the beginning of the i-th row, */
   /*   which can be treated as a one-dimensional array of length kk */
   /* Keep track of the total number of within-block tiegroups */
   /*   and the total tiesum: */
   tiesum=0.0;  totntg=0;
   for (i=0; i<nn; i++)  /* For each block or subject */
    { int ntg, tgsize[MAXK];
      /* Recall that for double array xx[][], xx[i] is the address */
      /*   of the beginning of the i-th row, which is stored in memory */
      /*   as a one-dimensional array. */
      /* Thus the `markranks()' call takes the values in the i-th row */
      /*   of xx[][] and put midranks in the i-th row of midr[][], */
      /*   which is what is required. */
      /* The function  makeranks()  is in statlib */
      makeranks(midr[i], xx[i],kk, tgsize,&ntg);
      /* The tie correction for the Friedman S statistic uses */
      /*   the total cumulative tiesum. The total number of tiegroups */
      /*   is for display. */
      for (k=0; k<ntg; k++)
       { int ts=tgsize[k];  tiesum+=(ts-1)*ts*(ts+1); }
      totntg+=ntg; }

   printf ("\nThe data with Friedman (within-block) ranks:\n");
   showranks(treatnames, xx,midr, kk,nn);

   /* The Friedman rank test: */
   /* The sum of squares of the treatment rank sums */
   fscore=0.0;
   for (j=0; j<kk; j++)
    { double rsum=0.0;
      for (i=0; i<nn; i++) rsum+= midr[i][j];
      fscore += rsum*rsum; }
   /* The Friedman statistic S (scaled to be asymptotically chi-square) */
   ssfrd = ss0 = (12*fscore)/(nn*kk*(kk+1)) - 3*nn*(kk+1);
   /* Is a tie correction needed? */
   if (totntg>0)
     ssfrd /= (1.0 - tiesum/((double)nn*(kk-1)*kk*(kk+1)));
   printf ("\nFriedman statistic S'=%g   P=%6.4f (%d df, %d tiegroup%s)\n",
     ssfrd, pvchisq(kk-1,ssfrd), kk-1, totntg, (totntg==1)?"":"s");
   /* Do simulation for Friedman statistic */
   printf ("\nSimulating the Friedman P-value using %d permutations:\n",
        nsims);
   printf ("  Friedman score (for simulations):  %g\n", fscore);
   /* Carry out the simulations using Friedman permutations */
   /* Warning: This randomizes the midr[][] matrix */
   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>=fscore)  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 = %7.5f  95%% CI (%6.4f, %6.4f, %6.4f)\n",
     nbig, nsims, pp, pp-pwid, pp, pp+pwid);
   if (totntg>0)
    { printf ("\nLarge-sample approximation WITHOUT tie correction:\n");
      printf ("  S=%g   P=%6.4f   (%d d.f.)\n",
        ss0, pvchisq(kk-1,ss0), kk-1); }
         }





/* Page's test for increasing alternatives */

void do_page(const double xx[MAXN][MAXK],
  int kk, int nn, const char *treatnames[MAXK])
 { int i,j,k, ns,nbig, totntg;
   /* Midranks of data in xx[MAXN][MAXK] */
   double midr[MAXN][MAXK];
   /* Miscellaneous doubles */
   double tiesum, pp,pwid, lpage,lpmean,obspage,zpage;

   /* Describe the data */
   printf ("%d treatments, %d subjects\n", kk,nn);
   /* Find Friedman ranks */
   /* See comments in do_freidman() */
   tiesum=0.0;  totntg=0;
   for (i=0; i<nn; i++)  /* For each block or subject */
    { int ntg, tgsize[MAXK];
      makeranks(midr[i], xx[i],kk, tgsize,&ntg);
      for (k=0; k<ntg; k++)
       { int ts=tgsize[k];  tiesum+=(ts-1)*ts*(ts+1); }
      totntg+=ntg; }

   printf ("\nThe data with Friedman (within-block) ranks:\n");
   showranks(treatnames, xx,midr, kk,nn);

   /* Compute the Page score from the midranks */
   lpage=0.0;
   for (j=0; j<kk; j++)
    { double rsum=0.0;
      for (i=0; i<nn; i++) rsum+= midr[i][j];
      lpage += (j+1)*rsum; }
   /* Find the large-sample normal approximation */
   lpmean=((double)nn*kk*(kk+1)*(kk+1))/4;
   zpage = (lpage - lpmean) / (kk*(kk+1)*sqrt(nn*(kk-1))/12.0);
   printf ("\nPage statistic L=%g   E(L)=%g   "
     "Z=%g   P=%6.4f (two-sided)\n",
     lpage, lpmean, zpage, 2*(1-normint(fabs(zpage))));
   /* Simulate the exact Page-test P-value: */
   printf ("  Simulating the Page P-value using %d permutations:\n",
     nsims);
   /* Convert to two-sided statistic */
   obspage=fabs(lpage-lpmean);
   /* Warning: This randomizes the midr[][] matrix */
   for (ns=0; ns<nsims; ns++)
    { double newscore;
      /* Shuffle the midranks for each subject */
      for (i=0; i<nn; i++)  shuffle_dub(midr[i],kk);
      /* Recompute the Page 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 += (j+1)*rsum; }
      /* Convert to two-sided statistic */
      newscore=fabs(newscore-lpmean);
      /* Newscore= int/4, so no floating-point roundoff error */
      if (newscore>=obspage)  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 = %7.5f  95%% CI (%6.4f, %6.4f, %6.4f)\n",
     nbig, nsims, pp, pp-pwid, pp, pp+pwid);
              }