/* H0:  I population proportions are equal
   Ha: not all equal  */
/* This program produces power estimates for ANOM and Pearson Chi-Squared Tests */
/* Required values to be filled in:
reps = XXXXX
I = XXXX
balanced = XXXX  lower case 'unbalanced' or 'balanced'
ng(i) for i = 1 to I, the sample sizes
pi(i) for i = 1 to I, the true population proportion
*/



data binom;
label reject05 = 'ANOM Power alpha = .05';
label reject10 = 'ANOM Power alpha = .10';
label reject01 = 'ANOM Power alpha = .01';
label chireject05 = 'P ChiSq Power alpha = .05';
label chiPowerAsp = 'Asymtotic P ChiSq Power';

label r051 = 'Reject Rate Pop 1';
label r05U1 = 'Reject Rate Pop 1 > UDL';
label r05L1 = 'Reject Rate Pop 1 < LDL';

label r052 = 'Reject Rate Pop 2';
label r05U2 = 'Reject Rate Pop 2 > UDL';
label r05L2 = 'Reject Rate Pop 2 < LDL';

label r053 = 'Reject Rate Pop 3';
label r05U3 = 'Reject Rate Pop 3 > UDL';
label r05L3 = 'Reject Rate Pop 3 < LDL';

label r054 = 'Reject Rate Pop 4';
label r05U4 = 'Reject Rate Pop 4 > UDL';
label r05L4 = 'Reject Rate Pop 4 < LDL';

label r055 = 'Reject Rate Pop 5';
label r05U5 = 'Reject Rate Pop 5 > UDL';
label r05L5 = 'Reject Rate Pop 5 < LDL';

label r056 = 'Reject Rate Pop 6';
label r05U6 = 'Reject Rate Pop 6 > UDL';
label r05L6 = 'Reject Rate Pop 6 < LDL';

label r057 = 'Reject Rate Pop 7';
label r05U7 = 'Reject Rate Pop 7 > UDL';
label r05L7 = 'Reject Rate Pop 7 < LDL';

label r058 = 'Reject Rate Pop 8';
label r05U8 = 'Reject Rate Pop 8 > UDL';
label r05L8 = 'Reject Rate Pop 8 < LDL';

label r059 = 'Reject Rate Pop 9';
label r05U9 = 'Reject Rate Pop 9 > UDL';
label r05L9 = 'Reject Rate Pop 9 < LDL';

label r0510 = 'Reject Rate Pop 10';
label r05U10 = 'Reject Rate Pop 10 > UDL';
label r05L10 = 'Reject Rate Pop 10 < LDL';

label r0511 = 'Reject Rate Pop 11';
label r05U11 = 'Reject Rate Pop 11 > UDL';
label r05L11 = 'Reject Rate Pop 11 < LDL';

label r0512 = 'Reject Rate Pop 12';
label r05U12 = 'Reject Rate Pop 12 > UDL';
label r05L12 = 'Reject Rate Pop 12 < LDL';

label r0513 = 'Reject Rate Pop 13';
label r05U13 = 'Reject Rate Pop 13 > UDL';
label r05L13 = 'Reject Rate Pop 13 < LDL';

label r0514 = 'Reject Rate Pop 14';
label r05U14 = 'Reject Rate Pop 14 > UDL';
label r05L14 = 'Reject Rate Pop 14 < LDL';

label r0515 = 'Reject Rate Pop 15';
label r05U15 = 'Reject Rate Pop 15 > UDL';
label r05L15 = 'Reject Rate Pop 15 < LDL';
reps = 100000;
I = 4;
balanced = 'unbalanced';
array ng(15);
ng(1) = 92;
ng(2) = 95;
ng(3) = 93;
ng(4) = 175;
ng(5) = 100;
ng(6) = 200;
ng(7) = 200;
ng(8) = 200;
ng(9) = 200;
ng(10) = 200;
ng(11) = 150;
ng(12) = 150;
ng(13) = 150;
ng(14) = 150;
ng(15) = 150;
array pi(15);
pi(1) = .70;
pi(2) = .73;
pi(3) = .69;
pi(4) = .85;
pi(5) = .75;
pi(6) = .668;
pi(7) = .75;
pi(8) = .69;
pi(9) = .69;
pi(10) = .70;
pi(11) = .1;
pi(12) = .1;
pi(13) = .1;
pi(14) = .1;
pi(15) = .1;
array nsize(15);
array p(15);
do iset = 1 to 15;
	if iset < = I then nsize(iset) = ng(iset);else nsize(iset) = 0;
	if iset < = I then p(iset) = pi(iset);else p(iset) = .;
end;
/*********  This next section calculates the asymptotic power for Pearson Chi Square */

xsum = 0;
bigNpow = 0;
array np(15);
do ix = 1 to I;
	np(ix) = ng(ix)*pi(ix); 
	xsum = xsum + np(ix);
	bigNpow = bigNpow + ng(ix);
	
end;
array Ep(15);
chisqPow = 0;

do wx = 1 to I ;
	Ep(wx) = xsum*ng(wx)/BigNpow;

	chisqPow = chisqPow + ((Ep(wx)-np(wx))**2)/Ep(wx) + ((Ep(wx)-np(wx))**2)/(ng(wx)-Ep(wx));
end;



CritChi05 = cinv(.95,I-1);
lambda = chisqPow;
chiPowerAsp = 1 - probchi(critchi05,I-1,lambda);

/***********************************************
This begins the series of replications for each scenario that has been established */

do rep = 1 to REPS;

reject10 = 0;reject05 = 0;reject01 = 0; ChiReject05 = 0;
/*below are critical values for the ANOM */
if balanced = 'balanced' then goto BAL100;
/*UnBalanced Critical Values*/
if I = 15 then do;
h01 = 3.40;
m05 = 2.93;
h10 = 2.70;
end;
if I = 14 then do;
h01 = 3.38;
m05 = 2.91;
h10 = 2.67;
end;
if I = 13 then do;
h01 = 3.36;
m05 = 2.88;
h10 = 2.65;
end;
if I = 12 then do;
h01 = 3.34;
m05 = 2.86;
h10 = 2.62;
end;
if I = 11 then do;
h01 = 3.32;
m05 = 2.83;
h10 = 2.59;
end;
if I = 10 then do;
h01 = 3.29;
m05 = 2.80;
h10 = 2.56;
end;
if I = 9 then do;
h01 = 3.26;
m05 = 2.77;
h10 = 2.52;
end;
if I =8 then do;
h01 = 3.23;
m05 = 2.73;
h10 = 2.48;
end;
if I =7 then do;
h01 = 3.19;
m05 = 2.68;
h10 = 2.43;
end;
if I =6 then do;
h01 = 3.14;
m05 = 2.63;
h10 = 2.38;
end;
if I = 5 then do;
h01 = 3.09;
m05 = 2.57;
h10 = 2.31;
end;
if I = 4 then do;
h01 = 3.02;
m05 = 2.49;
h10 = 2.23;
end;
if I = 3 then do;
h01 = 2.93;
m05 = 2.39;
h10 = 2.11;
end;
if I = 2 then do;
h01 = 2.81;
m05 = 2.24;
h10 = 1.95;
end;
if balanced = 'unbalanced' then goto SKIP100;
Bal100:
if I = 15 then do;
h01 = 3.40;
m05 = 2.93;
h10 = 2.69;
end;
if I = 14 then do;
h01 = 3.38;
m05 = 2.90;
h10 = 2.67;
end;
if I = 13 then do;
h01 = 3.36;
m05 = 2.88;
h10 = 2.65;
end;
if I = 12 then do;
h01 = 3.34;
m05 = 2.86;
h10 = 2.62;
end;
if I = 11 then do;
h01 = 3.32;
m05 = 2.83;
h10 = 2.59;
end;
if I = 10 then do;
h01 = 3.29;
m05 = 2.80;
h10 = 2.55;
end;
if I = 9 then do;
h01 = 3.26;
m05 = 2.76;
h10 = 2.52;
end;
if I =8 then do;
h01 = 3.22;
m05 = 2.72;
h10 = 2.47;
end;
if I =7 then do;
h01 = 3.19;
m05 = 2.68;
h10 = 2.42;
end;
if I =6 then do;
h01 = 3.14;
m05 = 2.62;
h10 = 2.36;
end;
if I = 5 then do;
h01 = 3.08;
m05 = 2.56;
h10 = 2.29;
end;
if I = 4 then do;
h01 = 3.01;
m05 = 2.47;
h10 = 2.19;
end;
if I = 3 then do;
h01 = 2.91;
m05 = 2.34;
h10 = 2.05;
end;
if I = 2 then do;
h01 = 2.58;
m05 = 1.96;
h10 = 1.65;
end;
skip100:
array x(15);
array phat(15);
avg = 0;
sumX = 0;
bigN = 0;
do gen = 1 to I;
	
	x(gen) = ranbin(-1,ng(gen),pi(gen));
	phat(gen) = x(gen)/ng(gen);
	sumX =sumX + x(gen);
	bigN = bigN + ng(gen);
end;
pbar = sumX/bigN;avg = pbar;
array UDL01(15);
array LDL01(15);
array UDL05(15);
array LDL05(15);
array UDL10(15);
array LDL10(15);
array r05(15);
array r05U(15);
array r05L(15);
do irj = 1 to 15;
	r05(irj) = 0;
	r05U(irj) = 0;
	r05L(irj) = 0;
end;
reject01 = 0;reject05 = 0;reject10 = 0;
do idl = 1 to I;
	udl01(idl) = avg + h01*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	ldl01(idl) = avg - h01*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	if phat(idl) > udl01(idl) then reject01 = 1;
	if phat(idl) < ldl01(idl) then reject01 = 1;

	udl05(idl) = avg + m05*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	ldl05(idl) = avg - m05*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	if phat(idl) > udl05(idl) then reject05 = 1;
	if phat(idl) < ldl05(idl) then reject05 = 1;
	if phat(idl) > udl05(idl) then r05(idl) = 1;
	if phat(idl) < ldl05(idl) then r05(idl) = 1;
	if phat(idl) > udl05(idl) then r05U(idl) = 1;
	if phat(idl) < ldl05(idl) then r05L(idl) = 1;

	udl10(idl) = avg + h10*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	ldl10(idl) = avg - h10*sqrt(avg*(1-avg))*sqrt((bigN-ng(idl))/(bigN*ng(idl)));
	if phat(idl) > udl10(idl) then reject10 = 1;
	if phat(idl) < ldl10(idl) then reject10 = 1;


end;
/*perform Chi Squared test */
array E(15);
chisqC = 0;

do w = 1 to I ;
	E(w) = sumX*ng(w)/BigN;

	chisqC = chisqC + ((E(w)-x(w))**2)/E(w) + ((E(w)-x(w))**2)/(ng(w)-E(w));
end;

CritChi05 = cinv(.95,I-1);
if chisqC > CritChi05 then ChiReject05 =1;

output;
end;

/*proc print data = binom; */
data IDdata; set binom; if rep > 1 then delete;
proc print data = IDdata; var reps I balanced h01 m05 h10;
proc print data = IDdata; var 
nsize1-nsize15
p1-p15 ;
title 'sampling scenario';run;
proc means data = binom mean; var reject05 ChiReject05 chiPowerAsp reject10 reject01
r051 r052 r053 r054 r055
r056 r057 r058 r059 r0510
r0511 r0512 r0513 r0514 r0515

r05U1 r05U2 r05U3 r05U4 r05U5
r05U6 r05U7 r05U8 r05U9 r05U10
r05U11 r05U12 r05U13 r05U14 r05U15

r05L1 r05L2 r05L3 r05L4 r05L5
r05L6 r05L7 r05L8 r05L9 r05L10
r05L11 r05L12 r05L13 r05L14 r05L15;output out = outsummary /* mean = mean*/;
title 'Power Results (Mean = Power; N = Reps)';run;

proc freq data = binom;
	tables reject05*ChiReject05 / agree;title 'Agreement between ANOM & Pearson ChiSquared';run;
/*data outsummary;set outsummary; if _STAT_ ^= 'MEAN' then delete;
proc print data = outsummary; title 'Power Summary';run;*/
quit;