/*prodclin_power.sas
This program automatically finds the sample size necessary for 0.8 power
for the joint significance test of mediation.  Needs prodclin.exe file
available for download at http://www.public.asu.edu/~davidpm/ripl/Prodclin*/
libname newsim 'C:\';
options nonotes noxwait;
title 'PRODCLIN';
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=.59;
%let effb=.26;
%let effcpr=0;
%let sampsize=123;
%let nreps=100000;
%let rho=0;
%let alpha=.05;

%macro newsim;
proc datasets; delete orig2;
/*Generates the data*/
%do j=1 %to &nreps;
	data data1;
 	do i=1 to &sampsize;
   		x=rannor(0);
   		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;

	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) keep=_TYPE_ m )
	OUT2(where=(_NAME_='m') rename=(m=varb) keep=_NAME_ m );
	drop _TYPE_ _NAME_;

	data data2; set orig;
	sea=sqrt(vara);
	seb=sqrt(varb);
	call symput('a',a);
	call symput('sea',sea);
	call symput('b',b);
	call symput('seb',seb);

data data3; 
*Change file address to match the location of the file prodclin.exe;
file "C:/Fritz/raw.txt";
a=&a; sea=&sea; b=&b; seb=&seb; rho=ρ alpha=α
put a @; put sea @; put b @; put seb @; put rho @; put alpha @;data data1;

*Change file address to match the location of the file prodclin.exe; 
X cd C:/Fritz;
*Change file address to match the location of the file prodclin.exe;
X call "C:/Fritz/prodclin.exe";

data data4; 
*Change file address to match the location of the file prodclin.exe;
infile "C:/Fritz/critval.txt";
input lcrit ucrit;
a=&a; sea=&sea; b=&b; seb=&seb; rho=ρ alpha=α
da=a/sea;
db=b/seb;
sedadb=sqrt(da*da+db*db+1);
dadb=da*db;
low=(lcrit-dadb)/(sedadb);
high=(ucrit-dadb)/(sedadb);
ab=a*b;
sobelse=sqrt(a*a*seb*seb+b*b*sea*sea);
prodlow=ab+low*sobelSE;
produp=ab+high*sobelSE;
if lcrit=99999 then prodlow='.';
if ucrit=99999 then produp='.';
if (lcrit<0 and ucrit>0) then acl95=0; else acl95=1;
keep acl95;
	proc append base=orig2 force data=data4;
%end;
run;

proc means data=orig2 noprint;
output out=b mean(acl95)=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;
run;
quit;

%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;