/* One sample nonparametric tests: */
/* Sign 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_SIGN"

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

/* Not the easiest way to enter paired data: */
/* 0 acts as a `sentinel' to flag the end of the data */
/* This program will count the number of pairs */

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);
/* Show 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,k,r, k1,k2, b1,b2;
   int nn,npos,nzero, ncases, tlarge;
   double sum1,sum2, xmean,xvar,xstderr;
   double tstat,tfac, zz, zdif,denom;
   double pval,p1val, twid, xsmean,xsvar,xsstdev;
   /* Hodges-Lehmann estimator and CI based on sign test */
   double xmed,alpha,cover, xsignlo,xsignhi;

   printf ("\nSign 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 sign 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;

   /* Carry out two classical tests */
   /* 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 3 versions of the sign test */
   printf ("Sign test:\n");
   printf ("  npos=%d positive out of nn=%d\n", npos,nn);
   /* Find the exact binomial P-value by counting cases */
   printf ("Exact sign test Pvalue by counting cases:\n");
   /* Assuming npos=2, the exact one-sided P-value is the sum */
   /*   for k=0 to 2 of (nn choose k), all divided by 2^nn */
   ncases=1+nn+(nn*(nn-1))/2;  /* Sum of combinatorial coefficients */
   /* pow(x,y) is x^y */
   denom=pow(2,nn);  p1val=ncases/denom;  pval=2*p1val;
   printf (
   "  P = (%d(0+)+%d(1+)+(%d*%d/2)(2+))/%g = %d/%g = %.4f (one-sided)\n",
      1, nn, nn,nn-1,denom, ncases,denom, p1val);
   printf ("  P = %.4f (two-sided)\n", pval);
   /* The exact binomial P-value using a function in statlib.c */
   r=npos;  if (nn-r>r) r=nn-r;  /* Find the right tail */
   p1val=binomtail(r,nn,0.50);   /* One-sided P-value */
   printf ("Using the `binomtail' function in statlib.c:\n");
   printf ("Exact Pval=%.4f (one-sided)  = %.4f (two-sided)\n",
     p1val, 2*p1val);
   /* Finally, the normal approximation to the sign-test P-value */
   printf ("Sign test (normal approximation):\n");
   xsmean=nn/2.0;  xsvar=nn/4.0;  xsstdev=sqrt(xsvar);
   printf ("  Npos-Xmean = %d-%g = %g  StdDev(Npos)=%.4f\n", 
     npos,xsmean, npos-xsmean, xsstdev);
   zz=(npos-xsmean)/xsstdev;
   printf ("  Z=%.4f (p0=1/2)  P=%.4f (two-sided)\n",  zz, normpval(zz));
   printf ("NOTE HOWEVER that the sign test does not take into account\n"
    "  the fact that the (+) values are smaller in absolute value.\n");

   /* Now Hodges-Lehmann estimators for the median difference: */
   /* The sign-test HL estimator is just the sample median */
   printf("\nHodges-Lehmann sign-test values:\n");
   for (i=0; i<nn; i++) xval[i]=xdif[i];
   /* Bubble sort xval[] in increasing order */
   for (i=0; i<nn; i++)  for (j=i+1; j<nn; j++)
     /* If out of order, switch the two values */
     if (xval[j]<xval[i])
      { double temp=xval[i];  xval[i]=xval[j];  xval[j]=temp; } 
   /* Find the median value or two median values of nn sorted values */
   /* If  0,1,2,3,4 (nn=5), the median is 2 at offset nn/2 */
   /* If  0,1,2,3 (nn=4), the median is (1+2)/2 at offsets (nn/2-1,nn/2) */
   k1=k2=nn/2;
   if ( (nn%2)==0) k1--;  /* Take k1=k2-1 if nn is even */
   xmed=(xval[k1]+xval[k2])/2;
   printf ("Median of %d values=%.3f\n", nn,xmed);
   /* Display with median value(s): */
   printf("Sorted differences (n=%d):\n", nn);
   showarray(xval,nn,k1,k2);

   printf ("\nFor HL-like confidence interval, see text p75-76:\n");
   /* Now find HL-like sign-test confidence interval for median */
   /* xval[] is still sorted in increasing order */
   /* Prob(B(n,p)>=j) = binomtail(j,n,p);   */
   /* Find k=j so that (1-alpha) is closest to 0.95 */
   printf ("Critical value b=b(alpha/2,%d,1/2) "
     "for 1-alpha closest to 0.95:\n", nn);
   zdif=2.0;  alpha=-1.0;  k=-1;
   for (j=nn-nn/2; j<=nn; j++)
    { double talpha2=binomtail(j,nn,0.50);
      printf ("  P[B(%d) >= %d] = %.4f   (%.1f%% CI)\n", 
        nn,j, talpha2, 100*(1-2*talpha2));
      if ( (zz=fabs(talpha2-0.025))<zdif)
       { k=j;  alpha=2*talpha2;  zdif=zz; }}
   printf ("Selecting closest to 95%%:\n");
   cover=1-alpha;  /* Coverage probability */
   printf ("  P[B(%d) >= %d] = %.4f   (%.1f%% CI)\n", 
     nn,k, alpha/2, 100*cover);
   /* Offsets and values for (1-alpha)*100% confidence interval */
   /* See Section 3.6 in text */
   b1=nn+1-k;  b2=k;
   printf ("Sign test CI offsets in (1,2,...,%d): (%d,%d)\n", nn,b1,b2);
   /* Convert ordinals to offsets */
   b1--;  b2--;
   printf ("Sign test CI offsets in (0,1,...,%d): (%d,%d):\n", nn-1,b1,b2);
   xsignlo=xval[b1];  xsignhi=xval[b2];
   /* Display the sorted values */
   showarray(xval,nn,b1,b2);
   printf ("Nonparametric %.1f%% CI: (%.3f, %.3f)\n",
     100*cover, xval[b1], xval[b2]);

   /* Find Student-t and normal 95% confidence intervals */
   printf ("\nTwo-sided Student-t critical values used for Student-t CI:\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 ("\nCIs for `average' value of differences:\n");
   twid=tfac*xstderr;
   printf ("  Sample mean(T):    %.4f   (%.4f, %.4f)  95%% CI\n", 
     xmean, xmean-twid, xmean+twid);
   twid=1.960*xstderr;
   printf ("  Sample mean(Z):    %.4f   (%.4f, %.4f)  95%% CI\n", 
     xmean, xmean-twid, xmean+twid);
   printf ("  HL for sign test:  %.4f   (%.4f, %.4f)  %.1f%% CI\n",
     xmed, xsignlo, xsignhi, 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 ("  %.3f%s", xx[i],str); }  printf("\n"); }