/*Sobel_power.sas
This program automatically finds the sample size necessary for 0.8 power
for the Sobel (1982) First-order test of mediation.*/
libname newsim 'C:\';
options nonotes symbolgen;

Filename nullog dummy 'C:\';
Proc Printto Log=nullog;

title 'Sobel ';

/*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=.59;
%let sampsize=667;
%let nreps=100000;

%macro newsim;
proc datasets; delete power;
run;
%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;

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

data a1; set out1;
if _MODEL_='MODEL1' and _TYPE_='PARMS';
a=x; keep a;
data a2; set out1;
if _MODEL_='MODEL1' and _NAME_='x';
vara=x; keep vara;

data b1; set out1;
if _MODEL_='MODEL2' and _TYPE_='PARMS';
b=m; keep b;
data b2; set out1;
if _MODEL_='MODEL2' and _NAME_='m';
varb=m; keep varb;

data data2; merge a1 a2 b1 b2;
ab=a*b;

/*Calculates the Sobel first order test*/
sobel=sqrt(a*a*varb+b*b*vara);


*confidence limits;
sobl95=ab-probit(.975)*sobel;
sobu95=ab-probit(.025)*sobel;


*If significance test is 2 tailed, then you need the absolute value of the calculated z;
if sobl95 le 0 and sobu95 ge 0 then hsob95=0; else hsob95=1;

proc append base=power force data=data2;
%END;
run;

proc means data=power noprint;
output out=b mean(hsob95)=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;