1 %let progname=corrbin.sas; 2 * ref: Tarone (1979) 3 Kupper and Haseman (1978) 4 Altham (1978) 5 * Tarone's Score Test Statistics 6 X_c := based on correlated/additive binomial dist. 7 X_m := based on multiplicative dist. 8 X_v := binomial variance test (df=K-1, K:= # samples). 9 * Several example datasets. 10 ***********************************************************************; 11 title1 "&progname: Tests for Extra-binomial Variation (Overdispersion)"; 12 ***********************************************************************; 13 proc iml; NOTE: IML Ready 14 reset noname; 15 ***********************************************************************; 16 start tarone_c(y, m); 17 /* 18 returns Tarone's test based on correlated binomial 19 signed version ~ N(0,1) 20 +ve means overdispersion, -ve underdispersion 21 */ 22 m1 = sum(m); 23 p = sum(y)/m1; 23 ! * estimated \pi under H_0; 24 A = sum(m#(m-1)); 25 S = ssq(y - m*p)/ (p*(1-p)); 26 test_n = S - m1; 26 ! * numerator; 27 test_d = sqrt(2*A); 27 ! * denominator; 28 29 ts = test_n/test_d; 29 ! * signed version ~ N(0,1); 30 return (ts); 31 2 The SAS System 17:45 Thursday, August 7, 2003 32 finish; NOTE: Module TARONE_C defined. 33 ******************************************************************************; 34 start tarone_m (y, m); 35 /* 36 returns Tarone's test based on multiplicative binomial 37 signed version ~ N(0,1) 38 +ve means overdispersion, -ve underdispersion 39 */ 40 41 m1 = sum(m); 42 p = sum(y)/m1; 42 ! * estimated \pi under H_0; 43 q = 1-p; 44 A = sum(m#(m-1)); 45 R = sum(y#(m-y)); 46 B = sum(m#(m-1)##2); 47 48 I_12=(q-p)*A; 49 I_11=p*q* (B + 2*p*q*(A - 2*B)); 50 I_22=m1/(p*q); 51 52 test_n= A*p*q - R; 52 ! /* numerator, OD -> small R */ 53 test_d= sqrt(I_11 - I_12*I_12/I_22); 53 ! /*denominator */ 54 55 ts = test_n/test_d; 55 ! * signed version ~ N(0,1); 56 return (ts); 57 58 finish; NOTE: Module TARONE_M defined. 59 60 ******************************************************************************; 61 start bvtest (y, m); 62 /* 63 returns binomial variance test == Pearson X2 64 */ 65 66 p = sum(y)/sum(m); 66 ! * estimated \pi under H_0; 67 r = y - m*p; 68 ts = sum(r#r/m) / (p*(1-p)); 69 return (ts); 70 71 finish; NOTE: Module BVTEST defined. 72 73 ***********************************************************************; 74 start result (y, m); 75 76 n = nrow(y)*ncol(y); 3 The SAS System 17:45 Thursday, August 7, 2003 77 print "Number of observations =" n; 78 x_c = tarone_c (y, m); 78 ! /* get the score test based on corr. bin*/ 79 x_m = tarone_m (y, m); 79 ! /* get the score test based on mult. bin*/ 80 x_v = bvtest (y, m); 80 ! /* Pearson's X2 */ 81 82 print "Taron's test - correlated binominal:"; 83 run pval(x_c); 84 85 print "Taron's test - multiplicative binominal:"; 86 run pval(x_m); 87 88 print "Pearson's X2 (binomial variance) test:"; 89 run pvalx2(x_v, n-1); 90 91 finish; NOTE: Module RESULT defined. 92 ***********************************************************************; 93 start pval (z); 94 95 p1 = 1-probnorm(z); 95 ! * 1-sided; 96 p2 = 2*min(p1, 1-p1); 96 ! * 2-sided; 97 x2 = z*z; 98 print "z =" z 99 ", p =" p1 "(one-sided)" 100 ", p =" p2 "(two-sided)." ; 101 102 finish; NOTE: Module PVAL defined. 103 ***********************************************************************; 104 start pvalx2 (x2, df); 105 106 p2 = 1-probchi(x2, df); 107 print "x2 =" x2 108 ", df =" df 109 ", p =" p2 "(two-sided)." ; 110 111 finish; NOTE: Module PVALX2 defined. 112 ***********************************************************************; 113 /* IML execution starts here */ 114 115 print "---------------------------------------------------------------------------"; 116 117 print "CTRL group"; 118 yc={32 34 37 39 39 39}; 118 ! * CTRL group (1); 119 mc={40 40 40 40 40 40}; 120 run result (yc, mc); 4 The SAS System 17:45 Thursday, August 7, 2003 120 ! * CTRL group (1); 121 122 print "---------------------------------------------------------------------------"; 123 124 print "FLU group"; 125 yf={32 34 35 36 40 40}; 125 ! * FLU group (2); 126 mf={40 40 40 40 40 40}; 127 run result (yf, mf); 127 ! * FLU group (2); 128 129 print "---------------------------------------------------------------------------"; 130 131 print "Artificial data. All m_i equal -> x^2_c == x^2_m"; 132 y={1 2 7 8 9}`; 133 m={9 9 9 9 9}`; 134 run result (yf, mf); 135 136 print "---------------------------------------------------------------------------"; 137 138 print "Crowder(1978), group 0 0" ; 139 y={10 23 23 26 17}`; 140 m={39 62 81 51 39}`; 141 run result (y, m); 142 143 print "---------------------------------------------------------------------------"; 144 145 print "KH(1978), ctl"; 146 y={0 2 0 0 0 0 0 1 2 1 }`; 147 m={5 6 7 7 8 8 8 9 9 10 }`; 148 run result (y, m); 149 150 print "---------------------------------------------------------------------------"; 151 152 print "* KH(1978), rx"; 153 y={0 2 1 0 2 3 0 4 1 6 }`; 154 m={5 5 7 8 8 8 9 9 10 10 }`; 155 run result (y, m); 156 157 print "---------------------------------------------------------------------------"; NOTE: Exiting IML. NOTE: The PROCEDURE IML printed pages 1-4. NOTE: PROCEDURE IML used: real time 0.07 seconds cpu time 0.07 seconds NOTE: SAS Institute Inc., SAS Campus Drive, Cary, NC USA 27513-2414 NOTE: The SAS System used: real time 0.32 seconds cpu time 0.28 seconds corrbin.sas: Tests for Extra-binomial Variation (Overdispersion) 1 17:45 Thursday, August 7, 2003 --------------------------------------------------------------------------- CTRL group Number of observations = 6 Taron's test - correlated binominal: z = 2.5833345 , p = 0.0048925 (one-sided) , p = 0.009785 (two-sided). Taron's test - multiplicative binominal: z = 2.5833345 , p = 0.0048925 (one-sided) , p = 0.009785 (two-sided). Pearson's X2 (binomial variance) test: x2 = 14.836364 , df = 5 , p = 0.0110849 (two-sided). --------------------------------------------------------------------------- FLU group Number of observations = 6 Taron's test - correlated binominal: z = 2.7023439 , p = 0.0034426 (one-sided) , p = 0.0068853 (two-sided). Taron's test - multiplicative binominal: z = 2.7023439 , p = 0.0034426 (one-sided) , p = 0.0068853 (two-sided). Pearson's X2 (binomial variance) test: corrbin.sas: Tests for Extra-binomial Variation (Overdispersion) 2 17:45 Thursday, August 7, 2003 x2 = 15.243438 , df = 5 , p = 0.0093713 (two-sided). --------------------------------------------------------------------------- Artificial data. All m_i equal -> x^2_c == x^2_m Number of observations = 6 Taron's test - correlated binominal: z = 2.7023439 , p = 0.0034426 (one-sided) , p = 0.0068853 (two-sided). Taron's test - multiplicative binominal: z = 2.7023439 , p = 0.0034426 (one-sided) , p = 0.0068853 (two-sided). Pearson's X2 (binomial variance) test: x2 = 15.243438 , df = 5 , p = 0.0093713 (two-sided). --------------------------------------------------------------------------- Crowder(1978), group 0 0 Number of observations = 5 Taron's test - correlated binominal: z = 1.4594887 , p = 0.0722153 (one-sided) , p = 0.1444306 (two-sided). Taron's test - multiplicative binominal: z = 1.9672878 , p = 0.024575 (one-sided) , p = 0.04915 (two-sided). corrbin.sas: Tests for Extra-binomial Variation (Overdispersion) 3 17:45 Thursday, August 7, 2003 Pearson's X2 (binomial variance) test: x2 = 9.7595422 , df = 4 , p = 0.0446789 (two-sided). --------------------------------------------------------------------------- KH(1978), ctl Number of observations = 10 Taron's test - correlated binominal: z = 0.2351618 , p = 0.4070416 (one-sided) , p = 0.8140831 (two-sided). Taron's test - multiplicative binominal: z = -0.136759 , p = 0.5543894 (one-sided) , p = 0.8912213 (two-sided). Pearson's X2 (binomial variance) test: x2 = 11.895436 , df = 9 , p = 0.2192686 (two-sided). --------------------------------------------------------------------------- * KH(1978), rx Number of observations = 10 Taron's test - correlated binominal: z = 2.5749211 , p = 0.0050131 (one-sided) , p = 0.0100263 (two-sided). Taron's test - multiplicative binominal: corrbin.sas: Tests for Extra-binomial Variation (Overdispersion) 4 17:45 Thursday, August 7, 2003 z = 1.8618542 , p = 0.0313118 (one-sided) , p = 0.0626236 (two-sided). Pearson's X2 (binomial variance) test: x2 = 19.029656 , df = 9 , p = 0.0249419 (two-sided). ---------------------------------------------------------------------------