/*joint_power.sas
This program automatically finds the sample size necessary for 0.8 power
for the joint significance test of mediation.*/
libname newsim 'C:\';
options nonotes;
title 'Joint';
Filename nullog dummy 'C:\';
Proc Printto Log=nullog;

/*Input values of interest
effa=alpha path, effb=beta path, effcpr=direct effect (tau prime path),
sampsize=sample size, and nreps=number of times the simulation should
run per iteration*/
%let effa=.14;
%let effb=.14;
%let effcpr=0;
%let sampsize=530;
%let nreps=100000;

%macro newsim;
proc datasets; delete orig2;
/*Generates the data*/
%do j=1 %to &nreps;
	data data1;
 	do i=1 to &sampsize;
   		x=rannor(0);
   		xsq=x*x;
   		m=&effa*x+rannor(0);
   		y=&effb*m+&effcpr*x+rannor(0);
   		output;
 	end;
 	run;

	/*Find parameter estimates based upon the generated data*/
	proc reg data=data1 covout outest=out1 noprint;
	model m = x;

	proc reg data=data1 covout outest=out2 noprint;
	model y = x m;

	proc reg data=data1 covout outest=out3 noprint;;
	model y = x;

	data orig;
	merge 
	OUT1(where=(_TYPE_='PARMS') rename=(x=a) keep=_TYPE_ x)
	OUT1(where=(_NAME_='x') rename=(x=vara) keep=_NAME_ x)
	OUT2(where=(_TYPE_='PARMS') rename=(m=b x=g) keep=_TYPE_ m x )
	OUT2(where=(_NAME_='m') rename=(m=varb) keep=_NAME_ m )
	OUT3(where=(_TYPE_='PARMS') rename=(x=xonly) keep=_TYPE_ x)
	OUT3(where=(_NAME_='x') rename=(x=varxonly) keep=_NAME_ x);
	drop _TYPE_ _NAME_;

	proc append base=orig2 force data=orig;
%end;
run;

/*Calculates the Joint test of mediation*/
*Confidence limits;
data orig3; set orig2;
al95=a-tinv(.975,(&sampsize-2))*sqrt(vara);
au95=a-tinv(.025,(&sampsize-2))*sqrt(vara);
bl95=b-tinv(.975,(&sampsize-3))*sqrt(varb);
bu95=b-tinv(.025,(&sampsize-3))*sqrt(varb);

*Sig of a;
if al95 le 0 and au95 ge 0 then ha95=0; else ha95=1;
*Sig of b;
if bl95 le 0 and bu95 ge 0 then hb95=0; else hb95=1;
*Overall sig;
hBK95=0;
if ha95=1 and hb95=1 then hBK95=1;


proc means data=orig3 noprint;
output out=b mean(hBK95)=mpower;

data c; set b;
npower=.8-mpower;
opower=abs(npower);
if npower>0 then csampsize=ceil((npower)*&sampsize);
else csampsize=floor((npower)*&sampsize);
	call symput('dpower',npower);
	call symput('power',mpower);
run;
data q;
	mpower=&power; dpower=&dpower; sampsize=&sampsize;
	proc print data=q;	
proc append base=iterate force data=q;
%mend;

%macro converge;
proc datasets; delete iterate;
%let itnum=0;
%do %until(&repeat=1);
	%newsim;
	run;
	%let itnum=%eval(&itnum+1);
	%let rpt=0;
	%if &itnum>=4 %then %do;
		data iterate2; set iterate;
		if _N_ ^= &itnum then delete;one=sampsize;i=1;keep one i;
		data iterate3; set iterate;
		if _N_^=(&itnum-2) then delete; two=sampsize;i=1;keep two i;
		data iterate4; set iterate;
		if _N_^=(&itnum-1) then delete; three=sampsize;i=1;keep three i;
		data iterate5; set iterate;
		if _N_^=(&itnum-3) then delete; four=sampsize;i=1;keep four i;
		data iterate6; merge iterate2 iterate3 iterate4 iterate5; by i;
		if one=two and three=four then rpt=1; else rpt=0;
		call symput('rpt',rpt);
		run;
	%end;
	data d; set c;
		call symput('opower',opower);
		call symput('csampsize',csampsize);
	run;
	%let apower=%sysevalf(&opower*&nreps);
	/*This is where the rule comes in that changes the sample
	size for the next iteration.*/
	%let sampsize=%sysevalf(&sampsize+&csampsize);
	%if &apower<100 %then %let repeat=1;
	%else %let repeat=%eval(&rpt*1);
	run;
	
%end;
run;
%mend;

%converge;
run;
quit;