%let progname=corrbin.sas;
* ref:  Tarone (1979)
        Kupper and Haseman (1978)
        Altham (1978)
* Tarone's Score Test Statistics
   X_c := based on correlated/additive binomial dist.
   X_m := based on multiplicative dist.
   X_v := binomial variance test (df=K-1, K:= # samples).
* Several example datasets.
***********************************************************************;
title1 "&progname: Tests for Extra-binomial Variation (Overdispersion)";
***********************************************************************;
proc iml;
reset noname;
***********************************************************************;
start tarone_c(y, m);
/*
returns Tarone's test based on correlated binomial
signed version ~ N(0,1)
+ve means overdispersion, -ve underdispersion
*/
  m1 = sum(m);
  p  = sum(y)/m1;                  * estimated \pi under H_0;
  A  = sum(m#(m-1));
  S  = ssq(y - m*p)/ (p*(1-p));
  test_n = S - m1;                       * numerator;
  test_d = sqrt(2*A);                    * denominator;

  ts  = test_n/test_d;           * signed version ~ N(0,1);
  return (ts);

finish;
******************************************************************************;
start tarone_m (y, m);
/*
returns Tarone's test based on multiplicative binomial
signed version ~ N(0,1)
+ve means overdispersion, -ve underdispersion
*/

  m1 = sum(m);
  p  = sum(y)/m1;                  * estimated \pi under H_0;
  q  = 1-p;
  A = sum(m#(m-1));
  R = sum(y#(m-y));
  B = sum(m#(m-1)##2);

  I_12=(q-p)*A;
  I_11=p*q* (B + 2*p*q*(A - 2*B));
  I_22=m1/(p*q);

  test_n= A*p*q - R;                     /* numerator, OD -> small R */
  test_d= sqrt(I_11 - I_12*I_12/I_22);   /*denominator */

  ts = test_n/test_d;           * signed version ~ N(0,1);
  return (ts);

finish;

******************************************************************************;
start bvtest (y, m);
/*
returns binomial variance test == Pearson X2
*/

  p  = sum(y)/sum(m);              * estimated \pi under H_0;
  r  = y - m*p;
  ts = sum(r#r/m) / (p*(1-p));
  return (ts);

finish;

***********************************************************************;
start result (y, m);

 n = nrow(y)*ncol(y);
 print "Number of observations =" n;
 x_c = tarone_c (y, m);   /* get the score test based on corr. bin*/
 x_m = tarone_m (y, m);   /* get the score test based on mult. bin*/
 x_v = bvtest   (y, m);   /* Pearson's X2 */

 print "Taron's test - correlated binominal:";
 run pval(x_c);

 print "Taron's test - multiplicative binominal:";
 run pval(x_m);

 print "Pearson's X2 (binomial variance) test:";
 run pvalx2(x_v, n-1);

finish;
***********************************************************************;
start pval (z);

 p1 = 1-probnorm(z);    * 1-sided;
 p2 = 2*min(p1, 1-p1);  * 2-sided;
 x2 = z*z;
 print "z =" z
       ",  p =" p1 "(one-sided)"
       ",  p =" p2 "(two-sided)." ;

finish;
***********************************************************************;
start pvalx2 (x2, df);

 p2 = 1-probchi(x2, df);
 print "x2 =" x2
       ",  df =" df
       ",  p =" p2 "(two-sided)." ;

finish;
***********************************************************************;
/* IML execution starts here */

print "---------------------------------------------------------------------------";

print "CTRL group";
yc={32 34 37 39 39 39}; * CTRL group (1);
mc={40 40 40 40 40 40};
run result (yc, mc);     * CTRL group (1);

print "---------------------------------------------------------------------------";

print "FLU group";
yf={32 34 35 36 40 40}; * FLU  group (2);
mf={40 40 40 40 40 40};
run result (yf, mf);     * FLU  group (2);

print "---------------------------------------------------------------------------";

print "Artificial data. All m_i equal -> x^2_c == x^2_m";
y={1 2 7 8 9}`;
m={9 9 9 9 9}`;
run result (yf, mf);

print "---------------------------------------------------------------------------";

print "Crowder(1978), group 0 0" ;
y={10 23 23 26 17}`;
m={39 62 81 51 39}`;
run result (y, m);

print "---------------------------------------------------------------------------";

print "KH(1978), ctl";
y={0 2 0 0 0 0 0 1 2 1 }`;
m={5 6 7 7 8 8 8 9 9 10 }`;
run result (y, m);

print "---------------------------------------------------------------------------";

print "* KH(1978), rx";
y={0 2 1 0 2 3 0 4 1 6 }`;
m={5 5 7 8 8 8 9 9 10 10 }`;
run result (y, m);

print "---------------------------------------------------------------------------";