/*
This macro tests the significance of the three path mediated effect.
It uses methods described in Taylor, A. B., MacKinnon, D. P., and Tein, J. -Y.
(in press). Tests of the three-path mediated effect. Organizational Research Methods.

How to use the macro:
1. Open your data set.
2. Change the input parameters for the macro at the bottom of this file
   to indicate the name of your data set, variable names, desired
   confidence interval percentage, etc. The input parameters are
   described in more detail at the bottom.
3. Run the macro, including both the code that makes up the macro
   (from %macro to %mend) and the macro call (%med3path).
*/

%macro med3path(dataset,x,m1,m2,y,nboot,seed,clpct,showreg);

* Convert confidence limit percent into p value at which to reject H0 and ;
* percentiles for confidence interval. ;
data _NULL_;
	* Save p value as a decimal (e.g., .05) ;
	pcrit=(100-&clpct)/100;
	call symput('pcrit',pcrit);
	* Save percentiles as percentages (e.g., 2.5 and 97.5) ;
	lowerpct=(100-&clpct)/2;
	upperpct=100-lowerpct;
	call symput('lowerpct',lowerpct);
	call symput('upperpct',upperpct);
run;	

* Get data and listwise delete to make sure the set of cases used ;
* is the same across the regression models. ;
data working; set &dataset;
	if (&x ne .) and (&m1 ne .) and (&m2 ne .) and (&y ne .);
	x=&x; m1=&m1; m2=&m2; y=&y;
	keep x m1 m2 y;
run;

* Find sample size ;
proc means data=working noprint;
	output out=meanout n(x)=nobs;
run;
data _NULL_; set meanout;
	call symput('nobs',nobs);
run;

%if &showreg=1 %then %let regprint=;
%else %let regprint=noprint;

* Estimate regression models. ;
* Regress M1 on X. ;
proc reg data=working outest=out1 tableout &regprint;
	model m1= x;
run;
data parm1; set out1;
	if _TYPE_='PARMS';
	b1=x; keep b1;
run;
data se1; set out1;
	if _TYPE_='STDERR';
	seb1=x; keep seb1;
run;
data t1; set out1;
	if _TYPE_='T';
	tb1=x; keep tb1;
run;
data pval1; set out1;
	if _TYPE_='PVALUE';
	pb1=x; keep pb1;
run;

* Regress M2 on X and M1. ;
proc reg data=working outest=out2 tableout &regprint;
	model m2= x m1;
run;
data parm2; set out2;
	if _TYPE_='PARMS';
	b2=m1; keep b2;
run;
data se2; set out2;
	if _TYPE_='STDERR';
	seb2=m1; keep seb2;
run;
data t2; set out2;
	if _TYPE_='T';
	tb2=m1; keep tb2;
run;
data pval2; set out2;
	if _TYPE_='PVALUE';
	pb2=m1; keep pb2;
run;

* Regress Y on X, M1, and M2. ;
proc reg data=working outest=out3 tableout &regprint;
	model y= x m1 m2;
run;
data parm3; set out3;
	if _TYPE_='PARMS';
	b3=m2; keep b3;
run;
data se3; set out3;
	if _TYPE_='STDERR';
	seb3=m2; keep seb3;
run;
data t3; set out3;
	if _TYPE_='T';
	tb3=m2; keep tb3;
run;
data pval3; set out3;
	if _TYPE_='PVALUE';
	pb3=m2; keep pb3;
run;

* Combine regression results and calculate standard errors of the ;
* three path mediated effect. ;
data regresult; merge parm1 parm2 parm3 se1 se2 se3 t1 t2 t3 pval1 pval2 pval3;
	med=b1*b2*b3;
	call symput('med',med);
	* Square coefficients and convert standard errors to variances for convenience. ;
	b1sq=b1*b1;
	b2sq=b2*b2;
	b3sq=b3*b3;
	varb1=seb1*seb1;
	varb2=seb2*seb2;
	varb3=seb3*seb3;
	* Calculate variance estimates. ;
	part1=b1sq*b2sq*varb3+b1sq*b3sq*varb2+b2sq*b3sq*varb1;
	part2=b1sq*varb2*varb3+b2sq*varb1*varb3+b3sq*varb1*varb2;
	part3=varb1*varb2*varb3;
	mvdvar=part1;
	unbivar=part1-part2+part3;
	exacvar=part1+part2+part3;
	* Find limits of confidence intervals. ;
	zcrit=probit(&upperpct/100);
	array variance mvdvar unbivar exacvar;
	array se mvdse unbise exacse;
	array me mvdme unbime exacme;
	array lcl mvdlcl unbilcl exaclcl;
	array ucl mvducl unbiucl exacucl;
	* Perform hypothesis test to get usual p value. ;
	array z mvdz unbiz exacz;
	array pval mvdp unbip exacp;
	do i=1 to 3;
		se[i]=sqrt(variance[i]);
		me[i]=zcrit*se[i];
		lcl[i]=med-me[i];
		ucl[i]=med+me[i];
		z[i]=med/se[i];
		pval[i]=2*(1-probnorm(z[i]));
	end;
run;

* Get random number seed from current time if seed = 0. ;
%if &seed < 1 %then %do;
	data _NULL_;
		seed=round(time()*10,1);
		call symput('seed',seed);
	run;
%end;

* Bootstrap. ;
proc surveyselect data=working seed=&seed out=boot sampsize=&nobs
	rep=&nboot method=urs outhits noprint;
run;

* Estimate regressions for bootstrap samples. ;
* Regress M1 on X. ;
proc reg data=boot outest=bootout1 noprint;
	by Replicate;
	model m1= x;
run;
data bootparm1; set bootout1;
	if _TYPE_='PARMS';
	bootb1=x; keep bootb1;
run;
* Regress M2 on X and M1. ;
proc reg data=boot outest=bootout2 noprint;
	by Replicate;
	model m2= x m1;
run;
data bootparm2; set bootout2;
	if _TYPE_='PARMS';
	bootb2=m1; keep bootb2;
run;
* Regress Y on X, M1, and M2. ;
proc reg data=boot outest=bootout3 noprint;
	by Replicate;
	model y= x m1 m2;
run;
data bootparm3; set bootout3;
	if _TYPE_='PARMS';
	bootb3=m2; keep bootb3;
run;

* Combine bootstrap results and get bootstrap estimates of the mediated effect. ;
data bootdist; merge bootparm1 bootparm2 bootparm3;
	bootmed=bootb1*bootb2*bootb3;
run;

* Get percentile bootstrap confidence limits. ;
proc univariate data=bootdist noprint;
	var bootmed;
	output out=pboot std=bootse pctlpts=&lowerpct &upperpct pctlpre=pb pctlname=lcl ucl;
run;

* Get p and z0 for bias-corrected bootstrap. ;
data findp; set bootdist;
	if bootmed > &med then p=1; else p=0;
run;
proc means data=findp noprint;
	output out=findp2 mean(p)=p;
run;
data _NULL_; set findp2;
	z0=probit(1-p);
	call symput('z0',z0);
run;

* Find percentiles to get for bias-corrected bootstrap. ;
data _NULL_;
	zp=probit(&upperpct/100);
	roundpoint=100/&nboot;
	bcbzlo=(2*&z0)-zp;
	bcbzup=(2*&z0)+zp;
	bcbplo=probnorm(bcbzlo);
	bcbpup=probnorm(bcbzup);
	bcbpctlo=round(bcbplo*100,roundpoint);
	bcbpctup=round(bcbpup*100,roundpoint);
	call symput('bcbpctlo',bcbpctlo);
	call symput('bcbpctup',bcbpctup);
run;

* Get values for bias-corrected bootstrap confidence limits. ;
proc univariate data=bootdist noprint;
	var bootmed;
	output out=bcboot pctlpts=&bcbpctlo &bcbpctup pctlpre=bcb pctlname=lcl ucl;
run;

* Get results for each method. ;
data result0;
	method='Individual coefficients:';
run;
%do i=1 %to 3;
	data result&i; set regresult;
		method="   test of b&i";
		estimate=b&i;
		stderr=seb&i;
		teststat=tb&i;
		pval=pb&i;
		lcl=.;
		ucl=.;
		if pb&i < &pcrit then rejh0='yes'; else rejh0='no';
	run;
%end;
data result4; set regresult;
	method='Joint significance test';
	estimate=.;
	stderr=.;
	teststat=.;
	pval=.;
	if (pb1 < &pcrit) and (pb2 < &pcrit) and (pb3 < &pcrit) then rejh0='yes'; else rejh0='no';
run;
%do i=5 %to 7;
	%if &i=5 %then %do; %let meth=mvd; %let methname=Multivariate delta SE; %end;
	%else %if &i=6 %then %do; %let meth=unbi; %let methname=Unbiased SE; %end;
	%else %if &i=7 %then %do; %let meth=exac; %let methname=Exact SE; %end;
	data result&i; set regresult;
		method="&methname";
		estimate=med;
		stderr=&meth.se;
		teststat=&meth.z;
		pval=&meth.p;
		lcl=&meth.lcl;
		ucl=&meth.ucl;
		if (&meth.lcl > 0) or (&meth.ucl < 0) then rejh0='yes'; else rejh0='no';
	run;
%end;
data result8; merge regresult pboot;
	method='Percentile bootstrap';
	estimate=med;
	stderr=.;
	teststat=.;
	pval=.;
	lcl=pblcl;
	ucl=pbucl;
	if (pblcl > 0) or (pbucl < 0) then rejh0='yes'; else rejh0='no';
run;
data result9; merge regresult bcboot;
	method='Bias-corrected bootstrap';
	estimate=med;
	stderr=.;
	teststat=.;
	pval=.;
	lcl=bcblcl;
	ucl=bcbucl;
	if (bcblcl > 0) or (bcbucl < 0) then rejh0='yes'; else rejh0='no';
run;

* Combine results. ;
data allresult;
	length method $ 24;
	set result0 result1 result2 result3 result4 result5 result6 result7 result8 result9;
	keep method estimate stderr teststat pval lcl ucl rejh0;
	label method='Method' estimate='Estimate' stderr='Standard Error' teststat='Test Statistic'
		pval='p value' lcl="Lower &clpct% limit" ucl="Upper &clpct% limit"
		rejh0="Reject H0 (p < &pcrit)";
run;

* Print results. ;
options label;
run;
proc print data=allresult label;
	title 'Tests of 3 path mediated effect. All are two-tailed.';
	title2 "&nboot bootstrap samples used. Random seed for bootstrapping = &seed.";
	id method;
run;
* Clear titles. ;
title;
title2;

%mend med3path;

/*
These are the macro input parameters:
 	Name of data set in which you want to test three path mediated effect
	Name of X variable
	Name of M1 variable
	Name of M2 variable
	Name of Y variable
	Number of bootstrap samples to take
	Random number seed for taking bootstrap samples (set to 0 to use current time)
	Confidence interval percentage (e.g., 95)
	Show individual regression model output (1=show, 0=hide)

For example, the input parameters might look like this for the variables in
the Keifer (2005) data used in the article:
%med3path(keiferdata,orgchange,workcond,negem,withdraw,1000,123456,95,1);

Replace the values in the line below with the appropriate names from your data. You
can also change the number of bootstrap samples, random number seed, confidence interval
percentage, and whether regression model output should be displayed.
*/
%med3path(datasetname,xvar,m1var,m2var,yvar,1000,123456,95,1);