/* One sample nonparametric tests: */
/* Wilcoxon signed-rank test for paired differences with associated */
/*   Hodges-Lehmann estimators and nonparametric confidence interval */

#include <stdio.h>
#include <stdlib.h>  /* Needed for exit(1) for errors */
#include <math.h>    /* Needed for sqrt() and pow() */
#include "statlib.h"

#define PROG_DATE "January 23, 2007"
#define PROG_NAME "ONESAMP_RANKS"

/* An array large enough to hold nn*nn values */
double xval[200];

/* Values from Table A4 for a WSR HL-like confidence interval */
#define NN     9
#define TALF   (2*0.027)
#define TTLEV  39

/* Not the easiest way to enter paired data: */
/* 0 acts as a `sentinel' to flag the end of the data */

double bef[20] =
   { 1.83,  0.50,  1.62,  2.48, 1.68, 1.88, 1.55, 3.06, 1.30, 0 };
double aft[20] =
   { 0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29, 0 };

double xdif[20];  /* Array of differences */
int ord[20], rank[20];  /* Two integer arrays */

/* Prototypes for functions used in this program: */
/* Return  - if x<0,  0 if x=0,  or + if x>0 */
char *xsgn (double x);
/* Display an array with one or two distinguished values */
void showarray (const double *xx, int num, int b1, int b2);

int main(void)
 { int i,j,r, b1,b2,k1,k2, nn,npos,nzero,nwalsh, ncases;
   int wsrp,wsrn,rasum, tlarge;
   double sum1,sum2, xmean,xvar,xstdev,xstderr;
   double tstat,tfac, zz, denom;
   double wsrmean, pval,p1val, twid;
   /* Hodges-Lehmann estimator and CI based on Wilcoxon signed ranks */
   double xmed,xwsrlo,xwsrhi,alpha, cover;

   printf (
   "\nWilcoxon signed-rank test and associated Hodges-Lehmann estimator\n"
   "  and confidence interval - written by XXXXX.\n");

   printf ("\nData: Hamilton Depression Scale Factor IV (Salsburg 1970)\n"
      "See the text (Hollander+Wolfe 2nd ed) Table 3.1, page 39\n"
      "Scale values before and after treatment for 9 patients.\n");
   printf ("Program %s - %s\n", PROG_NAME, PROG_DATE);

   /* Collect and display the data in arrays bef[] and aft[] */
   sum1=sum2=0.0;  npos=nzero=0;
   /* Recall i++ says to increase i by 1 */
   /* The loop stops when we hit the sentinel (value 0.0) */
   for (i=0; bef[i]>0.0; i++)
    { double xd = xdif[i] = aft[i]-bef[i];
      ord[i]=i+1;
      if (xd>0.0)  npos++;
      else if (xd<0.0)  ;  /* Do nothing */
      else nzero++;        /* Difference is zero */
      /* Collect sum of differences and sum of squares of differences */
      sum1+=xd;  sum2+=xd*xd; }
   nn=i;   /* The number of values collected */

   /* Display the data as read */
   printf ("\nN=%d observations:\n", nn);
   printf (" No.  Before  After  Diff(After-Before)  Sign:\n");
   for (i=0; i<nn; i++)
     printf (" %2d.  %6.3f %6.3f   %6.3f  %s\n",
       ord[i], bef[i],aft[i], xdif[i], xsgn(xdif[i]));
   printf ("Is there a significant decrease in the Depression Scale\n"
     "  after treatment?  A significant difference?\n");
   if (nn<=2)
    { printf ("N=%d: Not enough data!\n",nn);  exit(1); }
   /* Ties should be ignored with the signed-rank test */
   if (nzero>0)
    { printf ("%d difference value%s are zero:\n"
        "  Rerun the program without those values!\n",
        nzero, (nzero>1) ? "s" : "");  exit(1); }

   /* Sample mean and variance: */
   xmean=sum1/nn;
   xvar=(sum2-sum1*sum1/nn)/(nn-1);
   xstderr=sqrt(xvar/nn);  tstat=xmean/xstderr;

   /* Classical one-sample t-test */
   /* pvttest() is in statlib.c */
   printf ("\nClassical one-sample t-test (n=%d) for X=After-Before:\n", nn);
   printf ("  Mean=%.3f   Stdev(Mean)=%.3f\n", xmean,xstderr);
   printf ("  T=%.3f   P=%.4f   2-sided t-test, df=%d\n", 
     tstat, pvttest(nn-1,tstat), nn-1);

   /* Carry out 2 versions of the Wilcoxon signed rank test */
   printf ("\nWilcoxon signed rank analysis:\n");
   printf ("Assuming that there are no entries with difference=0\n"
     "  and that there are no tied absolute values.\n");
   printf ("Sorting by absolute differences (n=%d)\n",nn);
   printf ("  and assigning ranks based on abs.value of differences:\n");

   /* Bubble sort on absolute values of 4 arrays in parallel */
   for (i=0; i<nn; i++)  for (j=i+1; j<nn; j++)
     /* If out of order */
     if (fabs(xdif[j])<fabs(xdif[i]))
      { /* Switch values in 4 arrays in parallel */
        /* ord[i]=i+1 (originally) will be used to reverse the sort */
        /* rank[i] is still zero */
        int temp;  double xtemp;
        xtemp=bef[i];  bef[i]=bef[j];  bef[j]=xtemp;
        xtemp=aft[i];  aft[i]=aft[j];  aft[j]=xtemp;
        xtemp=xdif[i];   xdif[i]=xdif[j];    xdif[j]=xtemp;
        temp=ord[i];   ord[i]=ord[j];  ord[j]=temp; }
   /* Assign ranks based on absolute values of differences */
   for (i=0; i<nn; i++)  rank[i]=i+1;
   /* Display signed ranks in order of ranks or absolute values */
   for (i=0; i<nn; i++)
     printf (" %2d.  %6.3f %6.3f   %6.3f   %s %2d\n",
       i+1, bef[i], aft[i], xdif[i], xsgn(xdif[i]), rank[i]);

   /* Reverse bubble sort by a bubble sort on ord[] */
   printf ("Restore original order by bubble sort on  ord[]:\n");
   for (i=0; i<nn; i++)  for (j=i+1; j<nn; j++)
     /* If out of order by ord[] */
     if (ord[j]<ord[i])
      { /* Switch values in 5 arrays in parallel, now including ranks */
        int temp;  double xtemp;
        xtemp=bef[i];  bef[i]=bef[j];  bef[j]=xtemp;
        xtemp=aft[i];  aft[i]=aft[j];  aft[j]=xtemp;
        xtemp=xdif[i];   xdif[i]=xdif[j];    xdif[j]=xtemp;
        temp=ord[i];   ord[i]=ord[j];  ord[j]=temp;
        temp=rank[i];   rank[i]=rank[j];  rank[j]=temp; }
   /* Display the signed ranks */
   for (i=0; i<nn; i++)
     printf (" %2d.  %6.3f %6.3f   %6.3f   %s %2d\n",
       i+1, bef[i], aft[i], xdif[i], xsgn(xdif[i]), rank[i]);

   /* Calculate Wilcoxon signed ranks */
   printf ("Wilcoxon signed rank statistics:\n");
   wsrn=wsrp=0;
   for (i=0; i<nn; i++)
    { if (xdif[i]<0) wsrn+=rank[i];
      if (xdif[i]>0) wsrp+=rank[i]; }
   rasum=wsrn+wsrp;
   printf ("  T+=%d   T-=%d   Sum = %d = %d*%d/2\n", 
     wsrp,wsrn, rasum, nn,nn+1);
   /* Find the exact Wilcoxon-signed-rank P-value by counting cases */
   printf ("Exact P-value for Wilcoxon signed rank "
     "T+=%d by counting cases:\n", wsrp);
   ncases=1+5+4;  denom=pow(2,nn);  p1val=ncases/denom;  pval=2*p1val;
   printf ("  P = (%d(0)+%d(1)+%d(2))/%g = %d/%g = %.4f (one-sided)\n",
     1,5,4,denom, ncases,denom, p1val);
   printf ("  P = %.4f (two-sided)\n", pval);

   /* Find the P-value by a normal approximation to T+ or T- */
   printf ("Large-sample Z-test for WSR statistic (assuming no ties):\n");
   wsrmean=nn*(nn+1)/4.0;  /* E(T+)=E(T-) */
   printf ("  E(T-) = E(T+) = %d*%d/4 = %g\n", nn,nn+1,wsrmean);
   /* Now the Z-test variance */
   xvar=nn*(nn+1)*(2*nn+1)/24.0;  xstdev=sqrt(xvar);
   zz=(wsrn-wsrmean)/xstdev;
   printf ("  Z=(%d-%g)/sqrt(%d*%d*%d/24)=%g   P=%.4f (two-sided)\n", 
    wsrn,wsrmean, nn,nn+1,2*nn+1, zz, normpval(zz));

   /* The Hodges-Lehmann estimator for the median of xdif[] */
   /*   is the sample median of the Walsh averages */
   printf("\nHodges-Lehmann Wilcoxon signed-rank values:\n");
   nwalsh=(nn*(nn+1))/2;
   /* Put Walsh averages in xval[] */
   /* Here  xval[r++]=yy  stores xval[r]=yy  then increase r by 1 */
   /* In contrast,  xval[++r]=yy  would increase r FIRST and then */
   /*   store  xval[r+1]=yy;  */
   for (r=i=0; i<nn; i++)  for (j=i; j<nn; j++)
     xval[r++]=(xdif[i]+xdif[j])/2;
   /* Sort the Walsh averages in increasing order */
   for (i=0; i<nwalsh; i++)  for (j=i+1; j<nwalsh; j++)
     if (xval[j]<xval[i])
      { double temp=xval[i];  xval[i]=xval[j];  xval[j]=temp; } 
   /* Find the median of the Walsh averages */
   k1=k2=nwalsh/2;
   if ( (nwalsh%2)==0)  k1--;
   xmed=(xval[k1]+xval[k2])/2;
   printf("Median of Walsh averages(n=%d)=%.4f\n", nwalsh,xmed);
   printf("Sorted Walsh averages (n=%d):\n", nwalsh);
   showarray(xval,nwalsh,k1,k2);

   /* Now find HL-like Wilcoxon-signed-rank confidence interval */
   printf ("\nCritical value tt=T+(alpha/2,%d) for "
     "1-alpha closest to 0.95:\n", NN);
   cover=1-TALF;
   printf ("  P[T+(%d) >= %d] = %.4f   (%.1f%% CI)\n", 
     nn,TTLEV,TALF/2, 100*cover);
   /* Offsets and values for (1-alpha)*100% confidence interval */
   /* See Section 3.3 in text */
   b1=nwalsh+1-TTLEV;  b2=TTLEV;  alpha=TALF;
   printf ("WSR test CI offsets in (1,2,...,%d): (%d,%d)\n",
     nwalsh, b1,b2);
   /* Convert ordinals to offsets */
   b1--;  b2--;
   xwsrlo=xval[b1];  xwsrhi=xval[b2];
   printf ("WSR test CI offsets in (0,1,...,%d): (%d,%d)\n",
     nwalsh-1, b1,b2);
   /* Display the sorted values */
   showarray(xval,nwalsh,b1,b2);
   printf ("Nonparametric %.1f%% CI: (%.3f, %.3f)\n",
     100*cover, xval[b1], xval[b2]);

   /* Find Student-t  confidence intervals */
   printf ("\nTwo-sided Student-t critical values:\n");
   tfac=crit_tdist(nn-1, 0.025);
   tlarge=300;
   printf ("Tcrit(%d,0.05)=%.4f  Tcrit(%d,0.05)=%.4f\n",
     nn-1,tfac, tlarge,crit_tdist(tlarge,0.05));
   printf ("Mean=%.1f  Stderr=%.1f  Twid=%.1f\n",
     xmean, xstderr, tfac*xstderr);

   printf ("\nCIs of `average' value of differences:\n");
   twid=tfac*xstderr;
   printf ("  Sample mean(T):    %.4f   (%.4f, %.4f)  95%% CI\n", 
     xmean, xmean-twid, xmean+twid);
   printf ("  HL for WSR:        %.4f   (%.4f, %.4f)  %.1f%% CI\n",
     xmed, xwsrlo, xwsrhi, 100*(1-alpha));

   /* End of the program */
   /* Return 0 to say that everything went well. */
   return 0; }


char *xsgn (double x)
 { if (x<0.0) return "-";
   else if (x>0.0) return "+";
   else return "0"; }


void showarray (const double *xx, int num, int b1, int b2)
 { int i;
   /* Display the sorted values */
   for (i=0; i<num; i++)
    { char *str="";
      /* At most 7 values on a line */
      if (i>0 && (i%7)==0)  printf("\n");
      if (i==b1 || i==b2) str="(*)";
      printf ("  %.4f%s", xx[i],str); }  printf("\n"); }