%let f1_prevalence          = 0.20; ** given that prevelance for RF1 is 20%;
%let f2_prevalence          = 0.30; ** given that prevelance for RF2 is 30%;
%let f1pos_f2pos_prevalance = 0.05; ** given that prevelance for RF1 and RF2 is 5%;


%macro simlogistic(numits,n);

data summary1; run;  ** initializes dataset where summary of hypothesis tests will be stored;

ods listing close;  ** stop results from displaying in the output window temporarily;

%do i = 1 %to &numits;  ** repeat the process a total of &numits times;
data new1;

** solve for other prevalence estimates;
f1pos_f2neg_prevalance = &f1_prevalence - &f1pos_f2pos_prevalance;
f1neg_f2pos_prevalance = &f2_prevalence - &f1pos_f2pos_prevalance;
f1neg_f2neg_prevalance = 1.00 - &f1pos_f2pos_prevalance - f1pos_f2neg_prevalance - f1neg_f2pos_prevalance;

do i = 1 to &n;  ** do for each of the n subjects;
z = ranuni(-9);  ** generates a random number from a U(0,1) distribution;

** randomly puts a subject into a RF group based on prevalence estimates;
if 0 <= z < &f1pos_f2pos_prevalance  then group='1: F1+ F2+';
if &f1pos_f2pos_prevalance <= z < &f1pos_f2pos_prevalance + f1pos_f2neg_prevalance then group='2: F1+ F2-';
if &f1pos_f2pos_prevalance + f1pos_f2neg_prevalance <= z < &f1pos_f2pos_prevalance + f1pos_f2neg_prevalance + f1neg_f2pos_prevalance then group='3: F1- F2+';
if 1-f1neg_f2neg_prevalance <= z then group='4: F1- F2-';

if group in('1: F1+ F2+' '2: F1+ F2-') then f1=1; else f1=0;
if group in('1: F1+ F2+' '3: F1- F2+') then f2=1; else f2=0;

if group = '1: F1+ F2+' then risk=5.0*0.10;   **RR for this group was given as 5.0;
if group = '2: F1+ F2-' then risk=0.333*0.10; **I had to solve for this RR using the other given parameters;
if group = '3: F1- F2+' then risk=0.800*0.10; **I had to solve for this RR using the other given parameters;
if group = '4: F1- F2-' then risk=1.0*0.10;   **Baseline risk is 10 percent;

** randomly assigns disease based on risk;
z2 = ranuni(-9);
if z2 < risk then disease=1; else disease=0;

output;
end;

run;

/*
proc freq data=new1;
tables group f1 f2;
run;
proc means data=new1;
class group;
var risk;
run;
proc means data=new1;
class f1;
var risk;
run;
proc means data=new1;
class f2;
var risk;
run;
proc freq data=new1;
tables group*disease / nopercent nocol;
run;
*/

proc freq data=new1;
tables f1*f2*disease;
run;

proc logistic data=new1 descending;
class f1 f2;
model disease = f1 f2 f1*f2;
ods output type3=result1;
run;

** save the result of the interaction p-value from the logistic regression model;
data result1; set result1;
if effect='f1*f2';
if ProbChiSq < 0.05 then reject=1; else reject=0;
run;

data summary1; set summary1 result1;
if reject^=.;
run;
%end;

ods listing; ** turns back on the output window;

proc freq data=summary1;
tables reject;
run;

%mend;


%simlogistic(1000,1000); ** 99.4% power;
%simlogistic(1000,800);  ** 98.1% power;
%simlogistic(1000,400);  ** 83.0% power;
%simlogistic(1000,300);  ** 74.5% power;
%simlogistic(1000,200);  ** 47.5% power;