/* TwoProportionRandomizationTest.sas */
/* Leading Question re gun permits example */
/* size1 success1 sample size and number of successes for group 1 */
/* size2 success2 sample size and number of successes for group 2 */
/* G1 data for group 1, G2 data for group 2 */
proc iml;
size1 = 615;                       /* enter group one sample size   */
success1 = 463;                    /* enter group one success count */
size2 = 585;                       /* enter group two sample size   */
success2 = 403;                    /* enter group two success count */
fail1 = size1 - success1;
fail2 = size2 - success2;
G1=j(success1,1,1) // j(fail1,1,0);
G2=j(success2,1,1) // j(fail2,1,0);
phat1 = success1/size1;
phat2 = success2/size2;
obsratio = mean(G1) / mean(G2);
* print obsratio;            /* the observed ratio of the sample means */ 
* call randseed(12345);                      /* set random number seed */
alldata = G1 // G2;                   /* stack data in a single vector */
N1 = nrow(G1);  N = N1 + nrow(G2);            /* find the sample sizes */ 
NRepl = 9999;                       /* number of randomizations less 1 */
nulldist = j(NRepl,1);                  /* initialize a vector of ones */
x = sample(alldata, N//NRepl, "WOR");          /* create all resamples */
nulldist = x[,1:N1][,:] / x[,(N1+1):N][,:]; /* compute all mean ratios */
                    /* x[,1:N1] submatrix containing rows 1 through N1 */
           /* for a matrix A, a[:] is the average of the elements of A */
title "Histogram of Null Distribution of the ratio phat1/phat2";
title2 "Observed ratio phat1/phat2 at line, leading question re gun permit example";
refline = "refline " + char(obsratio) + " / axis=x lineattrs=(color=red);";
call Histogram(nulldist) other=refline;
pvalue = (1 + sum(abs(nulldist) >= abs(obsratio))) / (NRepl+1);
print "Observed ratio phat1/phat2, Pvalue for H_1:p1 notequal p2";
print phat1 size1 phat2 size2 obsratio pvalue;
/* ------------------------------------------------------------- */
/* HIV Vaccine example */
/* ------------------------------------------------------------- */
size1 = 3598;                       /* enter group one sample size   */
success1 = 241;                    /* enter group one success count */
size2 = 1805;                       /* enter group two sample size   */
success2 = 126;                    /* enter group two success count */
fail1 = size1 - success1;
fail2 = size2 - success2;
G1=j(success1,1,1) // j(fail1,1,0);
G2=j(success2,1,1) // j(fail2,1,0);
phat1 = success1/size1;
phat2 = success2/size2;
obsratio = mean(G1) / mean(G2);
* print obsratio;            /* the observed ratio of the sample means */ 
* call randseed(12345);                      /* set random number seed */
alldata = G1 // G2;                   /* stack data in a single vector */
N1 = nrow(G1);  N = N1 + nrow(G2);            /* find the sample sizes */ 
NRepl = 9999;                       /* number of randomizations less 1 */
nulldist = j(NRepl,1);                  /* initialize a vector of ones */
x = sample(alldata, N//NRepl, "WOR");          /* create all resamples */
nulldist = x[,1:N1][,:] / x[,(N1+1):N][,:]; /* compute all mean ratios */
                    /* x[,1:N1] submatrix containing rows 1 through N1 */
           /* for a matrix A, a[:] is the average of the elements of A */
title "Histogram of Null Distribution of the ratio phat1/phat2";
title2 "Observed ratio phat1/phat2 at line, HIV vaccine example";
refline = "refline " + char(obsratio) + " / axis=x lineattrs=(color=red);";
call Histogram(nulldist) other=refline;
pvalue = (1 + sum(abs(nulldist) >= abs(obsratio))) / (NRepl+1);
print "Observed ratio phat1/phat2, Pvalue for H_1:p1 notequal p2";
print phat1 size1 phat2 size2 obsratio pvalue;