/* How to handle nonparametric One-Way layouts: */
/* Kruskal-Wallis and Jonckheere-Terpstra tests. */

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

/* Compile with e.g.   */
/*    gcc -Wall OneWayLayout.c o OneWayLayout statlib.c   */
/* with both statlib.h and statlib.c in the same directory as OneWayLayout.c */

/* Upper bounds for
      MAXK: The number of treatment groups (k)
      MAXN: The number of observations in a treatment group (nn[i])
      MAXTOT:  The total number of observations: */

#define MAXK 5
#define MAXN 20
#define MAXTOT 100

/* The starting random-numbr seed. */

unsigned long wseed=408408;

/* The basic data is
     kk    : the number of treatment groups
     nn[i] : an integer array with the treatment-group sizes
     xx[i][j]:  the data itself, arranged as kk rows (0 <= i < kk)
       where the i-th rows has  nn[i] entries.
       Here  xx[i][j]  is the j-th entry in the i-th row.
     Here the rows are all of the same length, but this is not true
       in general. */

int kk=3;
int nn[MAXK] = { 5, 5, 5 };

double xx[MAXK][MAXN] = {
  { 53,  31, 139,    5,  16 },
  { 75,  64, 205,   97, 139 },
  { 54, 139, 194, 1253, 275 } };

/* NOTE:
   double xx[K][N];    (for example)
     sets aside room for K rows, where each row has room
     for N doubles.

 Syntactically, xx[][] is initializing by setting

   double xx[5][20] = {
     Row1,
     Row2,
     Row3 }; 

   with (for example)                           

     Row1 = { 53, 31, 139, 5, 16 }

   Note that
     (i) Within a row, all numbers except the last are followed
       by commas. The last number is NOT followed by a comma.
       (This is the same as initialized a one-dimensional array.)
     (ii) All rows except the last are followed by commas
       (OUTSIDE OF the braces that define the row). The last rom
       is NOT followed by a comma.
     (iii) The last row is followed by a closing brace and a
       semicolon.

   In global arrays (that is, defined outside of any function),
     uninitialed values are set to 0.0. The sets of braces around
     rows tell C to fill out the rows with zeroes.
   If you give values for all matrix entries, then you don't have
     to enclose the rows in braces. For example,

   double xx[3][3] = { 1, 2, 4, 2, 4, 8, 4, 8, 17 };

   has the same effect as
                                                
   double xx[3][3] = { 
     { 1, 2, 4 },
     { 2, 4, 8 },
     { 4, 8, 17 } };
                                 */


/* Prototypes: */

void do_kruskal_wallis (double *mrank, int ntot,
   int *nfirst, int *nmax, int kk, int *tgsize, int ntg);
double get_hobs (double *mrank,
   int *nfirst, int *nmax, int kk);
void do_sim_kruskal_wallis (double hobs, double *mrank, int ntot,
   int *nfirst, int *nmax, int kk, int nsims);

void do_jonckheere_terpstra (double *xdat, int ntot,
   int *nfirst, int *nmax, int kk, int *tgsize, int ntg);
void do_sim_jonckheere_terpstra (double *xdat, int ntot,
   int *nfirst, int *nmax, int kk, int nsims);

int main(int argc, char *argv[])
 { int i,j,r, ntot, ntg;
   int tgsize[MAXTOT];     /* Non-singlet tie-group sizes */
   double hobs;
   double xdat[MAXTOT];    /* The data in a linear order */
   double mrank[MAXTOT];   /* KW midranks of xdat[][] in a linear order */
   /* The range of the i-th treatment group in xdat[] and mrank[] */
   int nfirst[MAXK], nmax[MAXK];

   printf ("Nonparametric one-way layouts:\n");
   printf ("Kruskal-Wallis, Jonckheere-Terpstra, and "
     "multiple comparison tests.\n");

   /* Argc is the number of words on the command line.
      The arguments themselves are  argv[0], argv[1], ..., argv[argc-1]
      By convention, argc >=1 and argv[0] is the program name, so
      that a command-line argument was entered only if argc<1

    If a number is entered on the command line, use it as the
      starting seed. Recall that 0 means to initializes the seed
      from the system clock.
    If a number is not entered on the command line, use the
      previous value of wseed, which is as above. */

   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 the data */
   printf ("\nThe data:\n");
   /* Count the total number of observations */
   for (ntot=i=0; i<kk; i++)  ntot+=nn[i];
   printf ("Number of treatments: %d    "
     "Total number of observations: %d\n", kk,ntot);
   for (i=0; i<kk; i++)
    { printf ("Trtment #%d (n=%d):", i+1, nn[i]);
      for (j=0; j<nn[i]; j++) printf (" %6g", xx[i][j]);
      printf("\n"); }

   /* Carry out the Kruskal-Wallis test and compare with
        the parametric one-way ANOVA test
     Use a function in statlib.c to fill in midranks
     To do this, first
       (i) Copy all observations to a linear array xdat[] and
       (ii) fill in offsets noff[i] of the i-th row in xdat
       (iii) Find offsets [nfirst[i],nmax[i])  for i-th treatment group */
   r=0;
   for (i=0; i<kk; i++)
    { nfirst[i]=r;   /* Offset of row #i in xdat[] */
      for (j=0; j<nn[i]; j++)  xdat[r++]=xx[i][j];
      nmax[i]=r; }
   ntot=r;  /* The total number of observations */

   /* Call makeranks() in  statlib.c  to put midranks of the array  xdat[]
        in the array  mrank[]
      The midranks of the i-th treatment group will be at offsets
        nfirst[i] <= r < nmax[i]  in the array  mrank[]
      The function makeranks() reads an array of `ntot' doubles at xdat[]
        and writes the corresponding midranks at the same offsets in
        mrank[].  The array tgsize[] is filled in with the tie-group sizes
        for tie groups with more than one element, which are needed for
        tie corrections in large-sample approximations.
      The 5-th argument of makeranks() is a POINTER to a integer,
        which is filled in as the number of tie groups. */
   makeranks(mrank, xdat,ntot, tgsize,&ntg);

   /* Display the data with midranks: */
   printf ("\nData with midranks (RA=Rank-Average of treatment group):\n");
   for (i=0; i<kk; i++)       /* Sum over treatment groups */
    { /* Find the rank average */
      double rsum=0.0;
      if (nn[i]>0)
       { for (j=nfirst[i]; j<nmax[i]; j++) rsum += mrank[j];
         rsum=rsum/nn[i]; }
      printf (" #%d (RA=%4.1f): ", i+1, rsum);
      for (j=nfirst[i]; j<nmax[i]; j++)  /* Display the row */
        printf (" %4g(%4.1f)", xdat[j], mrank[j]);
      printf("\n"); }

   /* Carry out the large-sample Kruskal-Wallis test */
   /* This assume that there are `kk' treatment groups, that the midranks
        of the i-th treatment group are at offsets nfirst[i] <= r < nmax[i]
        in the array  mrank[], and that there are `ntg' tie-group sizes
        for tie groups with more than than one element at tgsize[] */    
   do_kruskal_wallis(mrank,ntot, nfirst,nmax,kk, tgsize,ntg);

   /* Simulate the Kruskal-Wallis P-values for 3 values of Nsims */
   printf (
   "\nNow SIMULATING the exact P-value of the Kruskal-Wallis test:\n");
   /* A simplified, equivalent form of the Kruskal-Wallis statistic
        is  Hobs = Sum RankSum_i^2/nn_i
      That is, Hobs is the same as H and HPrime except for constants
        that do not change under permutations, so that Hobs can
        be assumed to be the observed score */
   hobs=get_hobs(mrank,nfirst,nmax,kk);
   printf ("  Hobs = Sum RankSum_i^2/n_i = %g\n", hobs);
   /* The simulations randomize the order of mrank[], so that `hobs'
        cannot be calculated again */
   do_sim_kruskal_wallis(hobs, mrank,ntot, nfirst,nmax,kk, 1000);
   do_sim_kruskal_wallis(hobs, mrank,ntot, nfirst,nmax,kk, 10000);
   do_sim_kruskal_wallis(hobs, mrank,ntot, nfirst,nmax,kk, 100000);

   /* Are the treatment groups increasing in treatment-group order? */
   /* Is the data significant under a different criterion? */
   /* Carry out the Jonckheere-Terpstra test to find out */

   do_jonckheere_terpstra(xdat,ntot, nfirst,nmax,kk, tgsize,ntg);

   /* Simulate the Jonckheere-Terpstra P-value */
   printf (
   "\nNow SIMULATING the exact P-value of the Jonckheere-Terpstra test:\n");
   do_sim_jonckheere_terpstra(xdat,ntot, nfirst,nmax,kk, 100000);

   /* How many random numbers did we use? */
   /* That is, how many times did we update `rseed'? */
   /* Recall that we should worry about cycling for a LCRNG */
   /*   after around 100,000,000 random numbers. */

   show_randnums_used();

   return 0; }  /* Return with no errors */


/* Data arrays used in the Kruskal-Wallis test: */

void do_kruskal_wallis (double *mrank, int ntot,
   int *nfirst, int *nmax, int kk, int *tgsize, int ntg)
 { int i,j;
   double hsum,hh,hhnocorr, tiesum,tiecorr;
   printf ("\nCarrying out the Kruskal-Wallis test:\n");

   /* Compute Kruskal-Wallis H and (tie-corrected) HP */
   printf ("\nRecall that for the Kruskal-Wallis test:\n");
   printf ("  H = (12/(N*(N+1)))*Sum n_i*(RankAverage - (N+1)/2)^2\n");
   printf ("    (with no tie correction)\n");
   printf ("  HPrime = Hstat WITH tie correction.\n");

   /* Compute H  before tie correction */
   /* Uses the arrays  mrank[], nfirst[], and nmax[] */
   hsum=0.0;
   /* Sum over treatment groups, ignoring any empty groups */
   for (i=0; i<kk; i++) if (nn[i]>0)
    { /* Find the average rank */
      double rsum=0.0, rdif;
      for (j=nfirst[i]; j<nmax[i]; j++) rsum += mrank[j];
      rsum=rsum/nn[i];
      rdif=rsum-(double)(ntot+1)/2;
      hsum+=nn[i]*rdif*rdif; }
   hhnocorr=(12.0/(ntot*(ntot+1)))*hsum;

   /* Use  (tgsize[],ntg)  to find the Kruskal-Wallis tie correction */
   tiesum=0.0;
   for (i=0; i<ntg; i++)
    { int ts=tgsize[i];  tiesum += (ts-1)*ts*(ts+1); }
   /* Tie correction for variance */
   tiecorr = 1 - tiesum/((double)(ntot-1)*ntot*(ntot+1));
   hh=hhnocorr/tiecorr;

   printf ("Here:\n");
   printf ("  Rank-Average=%.1f  (Ntot=%d, ntreatments=%d)\n",
     (double)(ntot+1)/2, ntot, kk);
   printf ("  H(Prime)=%g   Tiesum=%g   Tiecorr=%g   H(NoCorr)=%g\n",
     hh, tiesum,tiecorr, hhnocorr);
   /* Note the function call to return a Chi-square P-value */
   printf ("  Large-sample Chisquare approx.:  P=%6.4f (HP=%g, df=%d)\n",
     pvchisq(kk-1,hh), hh, kk-1);
   printf ("  Without the tie correction:      P=%6.4f (H=%g,  df=%d)\n",
     pvchisq(kk-1,hhnocorr), hhnocorr, kk-1);
             }


/* Now, SIMULATED the P-value of the Kruskal-Wallis test: */

/* Compute Hobs = Sum RankSum_i^2/nn_i  given offsets and nsizes */

double get_hobs (double *mrank, int *nfirst, int *nmax, int kk)
 { int i,j;  double hobs=0.0;
   for (i=0; i<kk; i++) if (nmax[i]>nfirst[i])
    { double rsum=0.0;
      for (j=nfirst[i]; j<nmax[i]; j++)  rsum+=mrank[j];
      hobs+=rsum*rsum/(nmax[i]-nfirst[i]); }
   return hobs; }


void do_sim_kruskal_wallis (double hobs, double *mrank, int ntot,
   int *nfirst, int *nmax, int kk, int nsims)
 { int ns,nbig;  double pp,pwid;
   printf ("\nSimulating the true KW P-value with "
     "n=%d permutations:\n", nsims);
   /* Carry out permutations for the Kruskal-Wallis test */
   nbig=0;   /* `nbig' is the number of `successes' */
   for (ns=0; ns<nsims; ns++)
    { double hran;
      shuffle_dub(mrank,ntot);  /* Shuffle `mrank' array */
      hran=get_hobs(mrank,nfirst,nmax,kk);
      /* Allow for tiny roundoff error */
      if (hran>=hobs-1e-8) nbig++; }
   pp=(double)nbig/nsims;
   pwid = 1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P = %d/%s = %6.4f   95%% CI: (%6.4f, %6.4f)\n",
     nbig,commnum(nsims), pp, pp-pwid, pp+pwid); }


/* Return J = value of Jonckheere-Terpstra statistic */
/* Note that J is essentially the same as the Kendall (correlation)
     statistic between data values and the treatment group number */

double get_jtstat(double *xdat, int *nfirst, int *nmax, int kk)
 { int i,j;  double jj=0.0;
   /* Sum over pairs of treatment groups */
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
    { /* Add the number of positive Mann-Whitney differences */
      int r, ni1=nfirst[i], ni2=nmax[i];
      int s, nj1=nfirst[j], nj2=nmax[j];
      for (r=ni1; r<ni2; r++)  for (s=nj1; s<nj2; s++)
       { if (xdat[r]<xdat[s]) jj=jj+1;
         else if (xdat[r]==xdat[s]) jj=jj+0.5; }}
   return jj; }


double get_jmax(int *nfirst, int *nmax, int kk)
 { int i,j;  double jsum=0.0;
   for (i=0; i<kk; i++) for (j=i+1; j<kk; j++)
    { /* Add the maximum number of possible Mann-Whitney differences */
      int ni=nmax[i]-nfirst[i], nj=nmax[j]-nfirst[j];
      jsum += ni*nj; }
   return jsum; }


void do_jonckheere_terpstra (double *xdat, int ntot,
   int *nfirst, int *nmax, int kk, int *tgsize, int ntg)
 { int i;  double jtobs, jmax,jmean, jvar, zz;
   printf (
   "\nAre the treatment groups increasing in treatment-group order?\n"
   "Carry out the Jonckheere-Terpstra test to find out:\n");
   jtobs=get_jtstat(xdat,nfirst,nmax,kk);
   /* Find the P-value from the large-sample approximation: */
   /* Calculate the maximum possible J value (depends on nfirst,nmax) */
   jmax=get_jmax(nfirst,nmax,kk);
   jmean=jmax/2;
   printf ("\nStatistic J=%.1f   Min: %.2f  Mean: %.2f  Max: %.2f\n",
     jtobs, 0.0, jmean, jmax);
   /* Now the variance of J */
   jvar=ntot*ntot*(2*ntot+3);
   /* Tie correction for within treatment groups */
   for (i=0; i<kk; i++)
    { int ns=nmax[i]-nfirst[i];
      jvar-=ns*ns*(2*ns+3); }
   jvar/=72;  /* Means jvar=jvar/72 */
   printf ("  J = %g   E(J)=%g   Var(J)=%g (NOT tie-corrected)\n", 
     jtobs,jmean,jvar);
   /* There are three additional Jonckheere-Terpstra tiesums */
   if (ntg>0)
    { int g;
      double tiesum1, tiesum2, tiesum3, nsum, nsum2;
      tiesum1=tiesum2=tiesum3=0.0;
      for (g=0; g<ntg; g++)
       { int ts=tgsize[g];
         tiesum1 += ts*(ts-1)*(2*ts+5);
         tiesum2 += ts*(ts-1)*(ts-2);  tiesum3 += ts*(ts-1); }
      jvar -= tiesum1/72;
      nsum=nsum2=0.0;
      for (i=0; i<kk; i++)
       { int ns=nmax[i]-nfirst[i];
         nsum+=ns*(ns-1)*(ns-2);  nsum2+=ns*(ns-1); }
      nsum /= 36*ntot*(ntot-1)*(ntot-2);
      nsum2 /= 8*ntot*(ntot-1);
      jvar += (nsum*tiesum2 + nsum2*tiesum3); }
   printf ("  J = %g   E(J)=%g   Var(J)=%g (tie-corrected)\n", 
     jtobs,jmean,jvar);
   zz=(jtobs-jmean)/sqrt(jvar);
   printf ("  Large-sample approximation:  Z=%6.4f  P=%7.5f (2-sided, T.C.)\n",
     zz, normpval(zz));
        }


void do_sim_jonckheere_terpstra (double *xdat, int ntot,
    int *nfirst, int *nmax, int kk, int nsims)
 { int ns,nbig;
   double jtobs, jmax,jmean, jscore, pp,pwid;
   /* The observed Jonckheere-Terpstra J value: */
   jtobs=get_jtstat(xdat,nfirst,nmax,kk);
   /* Recalculate the maximum and mean of J under permutations */
   jmax=get_jmax(nfirst,nmax,kk);
   jmean=jmax/2;
   jscore=fabs(jtobs-jmean);
   printf ("\nSimulating the true JT P-value "
     "(n=%d simulations, Jscore=%.1f):\n",nsims,jscore);
   /* Since JT is symmetrically distrubuted about E(JT), */
   /*   |JT-E(JT)| is the appropriate score for a two-sided test. */
   nbig=0;    /* The number of successes */
   for (ns=0; ns<nsims; ns++)
    { double jj;
      shuffle_dub(xdat,ntot);  /* Shuffle `xdat' array */
      jj=get_jtstat(xdat,nfirst,nmax,kk);
      /* jj is an integer or half-integer, so no roundoff error */
      if (fabs(jj-jmean)>=jscore) nbig++; }
   pp=(double)nbig/nsims;
   pwid = 1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P = %d/%s = %7.5f   95%% CI: (%7.5f, %7.5f) (2-sided)\n",
     nbig,commnum(nsims), pp, pp-pwid, pp+pwid);
                 }