*--- program name:			bootstrap indirect effect;
*--- programmed by:			R Scherzer;
*--- creation date:			15 FEB 2017;
*--- comments:

Using MacKinnon's example from Introduction to Mediation
http://www.slideshare.net/smackinnon/introduction-to-mediation

for a mediation model with a continuous outcome (Y), predictor (X), and mediator (M).
This program estimates the indirect effect (a*b) and percent mediated (Pm = ab/c),
where:
a = effect of X on M
b = effect of M on Y
c' = direct effect of X on Y
c = total effect.

Bias-corrected confidence intervals are estimated using bootstrapping.	

;


*-- import data to reproduce estimates in mediation example;
filename testurl url "http://savvystatistics.com/wp-content/uploads/2015/03/crossroads.2015.data_.csv";
proc import file=testurl out=work.example dbms=csv replace;
run;


%MACRO bootindirect(data=example, alpha=0.05, replicates=1000, xvar=x, yvar=y, mvar=m);

proc surveyselect data= &data out=bootsample
     seed = 1234 method = urs
	 samprate = 1 outhits rep = &replicates;
run;

data bootsample2; set &data(in=a) bootsample;
if a then replicate=0; ** original sample;
run;

ods listing close;
ods output   ParameterEstimates=pe;
proc reg data=bootsample2; by replicate;
model &mvar = &xvar; 
model &yvar = &mvar &xvar; 
model &yvar = &xvar; 
run; quit;
ods listing ;

data a; set pe;
if model='MODEL1' and variable="&xvar";
rename estimate=a;
keep replicate estimate;
run;
data b; set pe;
if model='MODEL2' and variable="&mvar";
rename estimate=b;
keep replicate estimate;
run;
data cprime; set pe;
if model='MODEL2' and variable="&xvar";
rename estimate=cprime;
keep replicate estimate;
run;
data orig ab; merge a b cprime;
by replicate;
ab = a*b; 
c = ab + cprime;
Pm = ab / (ab + cprime);
if replicate=0 then output orig;
if replicate>0 then output ab;
run;

*store the estimated ab, Pm from original sample;
data _null_;
 set orig; 
 if replicate=0;
 call symput('abbar', put(ab,best.));
 call symput('Pmbar', put(Pm,best.));
run;

%let alpha1 = %sysevalf(1 - &alpha/2);
proc sql;
  select sum(ab<=&abbar)/count(ab) into :z0bar
  from ab;
quit;

data _null_;
  z0 = probit(&z0bar);
  zalpha = probit(&alpha1);
  p1 = put(probnorm(2*z0 - zalpha)*100, 3.0);
  p2 = put(probnorm(2*z0 + zalpha)*100, 3.0);
  output;
  call symput('a1', p1); call symput('a2', p2);
run;

* creating confidence interval, bias-corrected method for indirect effect, ab;
proc univariate data = ab alpha = .05 noprint;
  var ab;
  output out=pctl mean = abhat pctlpts=&a1 &a2 pctlpre = p pctlname = _lb _ub ;
run;

proc sql;
  select sum(pm<=&pmbar)/count(pm) into :z0bar_pm
  from ab;
quit;

data _null_;
  z0 = probit(&z0bar_pm);
  zalpha = probit(&alpha1);
  p1 = put(probnorm(2*z0 - zalpha)*100, 3.0);
  p2 = put(probnorm(2*z0 + zalpha)*100, 3.0);
  output;
  call symput('b1', p1); call symput('b2', p2);
run;

* creating confidence interval, bias-corrected method for proportion mediated, Pm;
proc univariate data = ab alpha = .05 noprint;
  var Pm;
  output out=pctl_pm mean = pmhat pctlpts=&b1 &b2 pctlpre = p pctlname = _lb _ub ;
run;

data ab2;
  set pctl;
  bias = abhat - &abbar;
  ab = &abbar;
run;
data pm2;
  set pctl_pm;
  bias = pmhat - &pmbar;
  pm = &pmbar;
run;

ods listing;
title 'Estimated indirect effect (ab) with 95pct bias-corrected confidence interval';
proc print data  = ab2;
  var ab bias p_lb p_ub;
run;

title 'Estimated percent mediated (Pm) with 95pct bias-corrected confidence interval';
proc print data  = pm2;
  var pm bias p_lb p_ub;
run;
title;


%MEND bootindirect;


*
This reproduces estimates for ab and confidence interval here:
http://www.slideshare.net/smackinnon/introduction-to-mediation;

options notes mprint;
%bootindirect(data=example, alpha=0.05, replicates=1000, xvar=bfi_c, yvar=health, mvar=behave);


/*
Estimated indirect effect (ab) with 95pct bias-corrected confidence interval

Obs       ab            bias      p_lb       p_ub

 1     0.20502    .001154822    0.15184    0.26031

Estimated pct mediated (Pm) with 95pct bias-corrected confidence interval

Obs       pm            bias      p_lb       p_ub

 1     0.43640    .004589566    0.31819    0.56095
*/