SPSS and SAS Macros for Bootstrapping Indirect Effects
Andrew F. Hayes, The
Kristopher J. Preacher,
Department of Psychology,
Click here for the SPSS macro command set
Click here for the SAS macro command set
The macro described here is used to generate bootstrap confidence intervals for indirect effects in a multiple mediator model of the form in Figure 1, where c is the total effect of X on Y, c′ is the direct effect of X on Y, and the specific indirect effect of X on Y through mediatior Mi is defined as aibi. For a paper currently undergoing review describing this application of bootstrapping, click here. This macro can do everything that our simple mediation macro can do (Preacher & Hayes, 2004, Behavior Research Methods, Instruments, and Computers, 36, 717-731), except this one allows for multiple mediators, statistical control of covariates, and it produces bias-corrected and bias-corrected and accelerated confidence intervals in addition to percentile-based confidence intervals. Thus, this macro is far superior to our earlier macro.
Figure 1. A multiple mediator model
As with all macros, the macro syntax must first be run in order to “activate” the macro. Once the macro is activated, it will remain active until you quit the SPSS or SAS session. The macro procedure is executed with the syntax and arguments described below.
The macro will use listwise deletion to remove cases with system missing data (generally represented in the data file as a period (“.”). Listwise deletion cannot be overridden. The macro does not accept user-defined missing data, so the data file should be purged of cases with user-defined missing data prior to running the macro
Syntax
for SPSS version
The SPSS version of the macro has the following syntax
INDIRECT y = dv/x = iv/m = mlist covlist/c = cov/boot = z/conf = ci/
percent = p/bc = b/bca = d
where dv is the name of the dependent variable (Y in figure 1), iv is the name of the independent variable (X in Figure 1), and mlist is a list of mediator variables (Mi in figure 1). These are the only arguments that are necessary for execution of the macro. If no other options are provided, the macro will estimate the paths in the model (using OLS) assuming no control variables, the number of bootstraps will be set to 1000, the confidence level will default to 95, and only bias corrected and adjusted confidence intervals for the indirect effects will be printed. For more options, the following arguments are used:
covlist is a list of covariate variables, and cov is the number of covariate variables in the covlist list, z is the number of desired bootstrap resamples desired in increments of 1000 (e.g., boot = 2000 yields 2000 bootstrap resamples), ci is the desired confidence level, p is set to 1 to print percentile-based confidence intervals, b is set to 1 to print bias-corrected confidence intervals, and d is set to 0 to disable printing of bias corrected and adjusted confidence intervals. If any of these arguments are not provided, default values will be used (the defaults are c = 0, z = 1000, ci = 95, p = 0, b =0, d = 1). Note that if c is set to 0 it is assumed that there are no variables listed in covlist, and all variables listed after “m =” are treated as potential mediators.
For example, the command below will estimate the total, direct, and indirect effects of TREAT on SATIS through SESSIONS and AFFECT, controlling for SEX and AGE, using 5000 bootstrap resamples to generate 90% confidence intervals, printing percentile and bias-corrected confidence intervals only.
INDIRECT Y = SATIS/X = TREAT/M = SESSIONS AFFECT SEX AGE/C = 2/BOOT = 5000/CONF = 90/PERCENT = 1/BC = 1/BCA = 0.
The z bootstrap estimates will be saved to a data file named INDIRECT.SAV. You can use this to examine the bootstrap sampling distributions.
The macro command set is below. Do not cut and paste the macro command set on this web page into SPSS, as this may not work. Instead, download the syntax file by clicking here.
Syntax
for SAS version
The SAS version of the macro functions identically to the SPSS
version. The macro accepts missing data
as “.” and will use listwise deletion to exclude cases with missing data. After the SAS macro command set is run, the
macro syntax is
%indirect (data=filename,y=dv,x=iv,m=mlist
covlist,c=cov,boot=z,conf=ci,
percent=p,bc=b,bca=d);
where filename is the name of an SPSS file name and the
remaining arguments are defined as in the SPSS instructions above. The SAS macro command corresponding to the
SPSS example above is
%indirect (data=example,y=satis,x=treat,m=sessions
affect sex age,c=2,boot=5000,conf=90,
percent=1,bc=1,bca=0);
where example is the name of the SAS datafile. The
z bootstrap estimates will be saved in a SAS data file named bootstp.
The macro command set is below. Do not cut and paste the macro command set on this web page into SAS, as this may not work. Instead, download the syntax file by clicking here.
Macro Command Sets
SPSS (v1.0,
DEFINE INDIRECT (y = !charend('/')/x = !charend('/')/m = !charend('/')/c=!charend('/')
!default(0)/boot =!charend('/') !default(1000)/conf = !charend('/')!default(95)/
percent = !charend('/') !default(0)/bc = !charend('/') !default(0)/bca = !charend('/')
!default(1)).
PRESERVE.
SET LENGTH = NONE.
SET MXLOOPS = 10000001.
SET SEED = RANDOM.
MATRIX.
get dd/variables = !y !x !m/names = nm/MISSING = OMIT.
compute nm = t(nm).
compute n = nrow(dd).
compute nv = ncol(dd).
compute nc = !c.
compute con = make(n,1,1).
compute dat2 = dd.
compute dat = dd.
compute bzx = make(nv-2-nc,1,0).
compute bzxse = make(nv-2-nc,1,0).
compute b=make((nv-1-nc),(nv-1-nc),0).
compute p0=-.322232431088.
compute p1 = -1.
compute p2 = -.342242088547.
compute p3 = -.0204231210245.
compute p4 = -.0000453642210148.
compute q0 = .0993484626060.
compute q1 = .588581570495.
compute q2 = .531103462366.
compute q3 = .103537752850.
compute q4 = .0038560700634.
compute conf = rnd(!conf).
compute lowalp = 0.5*(1-(conf/100)).
compute upalp = 0.5*(1+(conf/100)).
compute zbca = {lowalp; upalp}.
do if (!boot > 999).
compute btn = trunc(!boot/1000)*1000.
else.
compute btn = 1000.
end if.
compute blowp = trunc(lowalp*btn).
do if (blowp < 1).
compute blowp = 1.
end if.
compute bhighp = trunc((upalp*btn)+1).
do if (bhighp > btn).
compute bhighp = btn.
end if.
compute indeff = make(n+1+btn,nv-1-nc,0).
loop #d = 1 to (n+1+btn).
do if (#d = (n+2)).
compute dat = dat2.
compute con = make(n,1,1).
end if.
do if (#d > 1 and #d < (n+2)).
do if (#d = 2).
compute con = make((n-1),1,1).
compute dat = dat2(2:n,:).
else if (#d = (n+1)).
compute dat = dat2(1:(n-1),:).
else.
compute dat = {dat2(1:(#d-2),:);dat2((#d:n),:)}.
end if.
end if.
do if (#d > (n+1)).
compute v=trunc(uniform(n,1)*n)+1.
compute dat(:,1:nv) = dat2(v,1:nv).
end if.
compute x = dat(:,2).
compute m = dat(:,3:(nv-nc)).
compute y = dat(:,1).
compute xz = dat(:,2:nv).
compute xo = {con,x}.
do if (nc > 0).
compute c = dat(:,(nv-nc+1):nv).
compute xo = {xo, c}.
end if.
loop #k = 3 to (nv-nc).
compute ytmp = dat(:,#k).
compute bzxt = inv(t(xo)*xo)*t(xo)*ytmp.
compute bzx((#k-2),1)=bzxt(2,1).
do if (#d = 1).
compute mse=csum((ytmp-(xo*bzxt))&**2)/(n-2-nc).
compute olscm=(mse*inv((t(xo)*xo))).
compute bzxse((#k-2),1)=sqrt(olscm(2,2)).
end if.
end loop.
do if (#d = 1).
do if (nc > 0).
compute cnt = dd(:,(
compute xo = {con,x,cnt}.
else.
compute xo = {con,x}.
end if.
compute byx = inv(t(xo)*xo)*t(xo)*y.
compute mse=csum((y-(xo*byx))&**2)/(n-2-nc).
compute olscm=(mse*inv((t(xo)*xo))).
compute byxse = sqrt(olscm(2,2)).
compute byx = byx(2,1).
end if.
compute xzo = {con,xz}.
compute byzx = inv(t(xzo)*xzo)*t(xzo)*y.
compute byzx2 = byzx(3:(nv-nc),1).
do if (#d = 1).
compute mse=csum((y-(xzo*byzx))&**2)/(n-nv).
compute covmat=mse*inv(t(xzo)*xzo).
compute olscm=diag(covmat).
compute sse = mse*(n-nv).
compute sst = csum((y-(csum(y)/n))&**2).
compute r2 = 1-(sse/sst).
compute ar2 = 1-(mse/(sst/(n-1))).
compute fr = ((n-nv)*r2)/((1-r2)*ncol(xz)).
compute pfr = 1-fcdf(fr,ncol(xz),(n-nv)).
do if (nc > 0).
compute bcon = byzx((nv-nc+1):nv,1).
compute bconse = sqrt(olscm((nv-nc+1):nv,1)).
end if.
compute byzx2se = sqrt(olscm(3:(nv-nc),1)).
compute cprime = byzx(2,1).
compute cprimese = sqrt(olscm(2,1)).
end if.
compute indeff2 = (bzx&*byzx2).
compute zs = (bzx&/bzxse)&*(byzx2&/byzx2se).
compute temp = t({csum(indeff2); indeff2}).
compute indeff(#d,:) = temp.
do if (#d = 1).
compute vs = nm(1:(nv-nc),1).
print vs/title = "Dependent, Independent, and Proposed Mediator Variables:"/
rlabels = "DV =" "IV = " "MEDS = "/format a8.
do if (nc > 0).
compute vs = nm((nv-nc+1):nv,1).
print vs/title = "Statistical Controls:"/rlabels = "CONTROL="/format a8.
end if.
compute nms = nm(3:(nv-nc),1).
compute te = bzx&/bzxse.
compute df = n-2-nc.
compute p = 2*(1-tcdf(abs(te), df)).
compute bzxmat = {bzx, bzxse,te,p}.
compute b(2:(nv-1-nc),1)=bzx.
compute se2 = bzxse&*bzxse.
print bzxmat/title = "IV to Mediators (a paths)"/rnames = nms/clabels "Coeff"
"se" "t" "p"/format f9.4.
compute te = byzx2&/byzx2se.
compute df = n-nv.
compute p = 2*(1-tcdf(abs(te), df)).
compute byzx2mat={byzx2, byzx2se, te, p}.
compute df = n-nv.
print byzx2mat/title = "Direct Effects of Mediators on DV (b paths)"/
rnames = nms/clabels "Coeff" "se" "t" "p"/format f9.4.
compute te = byx&/byxse.
compute df = n-2-nc.
compute p = 2*(1-tcdf(abs(te), df)).
compute byxmat = {byx, byxse, te, p}.
compute xnm = nm(2,1).
compute df = n-2-nc.
print byxmat/title = "Total Effect of IV on DV (c path)"/rnames = xnm/
clabels "Coeff" "se" "t" "p"/format f9.4.
compute te = cprime&/cprimese.
compute df = n-nv.
compute p = 2*(1-tcdf(abs(te), df)).
compute cprimmat = {cprime, cprimese, te, p}.
compute df = n-nv.
print cprimmat/title = "Direct Effect of IV on DV (c' path)"/rnames = xnm/
clabels "Coeff" "se" "t" "p"/format f9.4.
do if (nc > 0).
compute df = n-nv.
compute nms = nm((
compute te = bcon&/bconse.
compute p = 2*(1-tcdf(abs(te), df)).
compute bconmat = {bcon, bconse,te,p}.
compute df = n-nv.
print bconmat/title = "Partial Effect of Control Variables on DV"/
rnames = nms/clabels "Coeff" "se" "t" "p"/format f9.4.
end if.
compute dvms = {r2, ar2, fr, ncol(xz), (n-nv), pfr}.
print dvms/title = "Fit Statistics for DV Model"/clabels "R-sq" "Adj R-sq"
"F" "df1" "df2" "p"/format F9.4.
end if.
end loop.
compute lvout = indeff(2:(n+1),:).
compute tdotm = csum(lvout)/n.
compute tm = (make(n,ncol(lvout),1))*mdiag(tdotm).
compute topa = csum((((n-1)/n)*(tm-lvout))&**3).
compute bota = 6*sqrt((csum((((n-1)/n)*(tm-lvout))&**2)&**3)).
compute ahat = topa&/bota.
compute nms = {"Total"; nm(3:(nv-nc),1)}.
compute indsam = t(indeff(1,:)).
compute boot = indeff((n+2):nrow(indeff),:).
compute mnboot = t(csum(boot)/btn).
compute se = (sqrt(((btn*cssq(boot))-(csum(boot)&**2))/((btn-1)*btn))).
do if (btn > 1).
save boot/outfile = estimate.sav/names = nms.
compute nnn = make(1,(nv-1-nc),-999).
compute boot = {nnn;boot}.
loop #e = 1 to (nv-1-nc).
loop #i = 2 to (btn+1).
compute ix = boot(#i,#e).
loop #k= #i to 2 by -1.
compute k = #k.
do if (boot(#k-1,#e) > ix).
compute boot(#k,#e)=boot(#k-1,#e).
else if (boot(#k-1,#e) <= ix).
BREAK.
end if.
end loop.
compute boot(k,#e)=ix.
end loop.
end loop.
compute boot = boot(2:(btn+1),:).
end if.
compute xp = make((nrow(mnboot)+2),1,0).
loop i = 1 to (nrow(mnboot)+2).
do if (i <= nrow(mnboot)).
compute pv = (boot(:,i) < indsam(i,1)).
compute pv = csum(pv)/btn.
else.
compute pv = zbca((i-nrow(mnboot)),1).
end if.
compute p = pv.
do if (pv > .5).
compute p = 1-pv.
end if.
compute y5=sqrt(-2*ln(p)).
compute xp(i,1)=y5+((((y5*p4+p3)*y5+p2)*y5+p1)*y5+p0)/
((((y5*q4+q3)*y5+q2)*y5+q1)*y5+q0).
do if (pv <= .5).
compute xp(i,1) = -xp(i,1).
end if.
end loop.
compute bbb = nrow(mnboot).
compute zz = xp(1:bbb,1).
compute zlo = zz + ((zz+xp((bbb+1),1))&/(1-t(ahat)&*(zz+xp((bbb+1),1)))).
compute zup = zz + ((zz+xp((bbb+2),1))&/(1-t(ahat)&*(zz+xp((bbb+2),1)))).
compute ahat = 0.
compute zlobc = zz + ((zz+xp((bbb+1),1))&/(1-t(ahat)&*(zz+xp((bbb+1),1)))).
compute zupbc = zz + ((zz+xp((bbb+2),1))&/(1-t(ahat)&*(zz+xp((bbb+2),1)))).
compute zlo = cdfnorm(zlo).
compute zup = cdfnorm(zup).
compute zlobc = cdfnorm(zlobc).
compute zupbc = cdfnorm(zupbc).
compute blow = trunc(zlo*(btn+1)).
compute bhigh = trunc(zup*(btn+1))+1.
compute blowbc = trunc(zlobc*(btn+1)).
compute bhighbc = trunc(zupbc*(btn+1))+1.
compute lowbca = make(nrow(blow),1,0).
compute upbca = lowbca.
loop i = 1 to nrow(blow).
do if (blow(i,1) < 1).
compute blow(i,1) = 1.
end if.
compute lowbca(i,1)=boot(blow(i,1),i).
do if (bhigh(i,1) > btn).
compute bhigh(i,1) = btn.
end if.
compute upbca(i,1)=boot(bhigh(i,1),i).
end loop.
compute lowbc = make(nrow(blow),1,0).
compute upbc = lowbca.
loop i = 1 to nrow(blowbc).
do if (blowbc(i,1) < 1).
compute blowbc(i,1) = 1.
end if.
compute lowbc(i,1)=boot(blowbc(i,1),i).
do if (bhighbc(i,1) > btn).
compute bhighbc(i,1) = btn.
end if.
compute upbc(i,1)=boot(bhighbc(i,1),i).
end loop.
print/title = "********************************************************************".
print/title = " BOOTSTRAP RESULTS FOR INDIRECT EFFECTS".
compute res = {indsam, mnboot,(mnboot-indsam), t(se)}.
print res/title = "Total and Specific Indirect Effect(s) of IV on DV via MEDS"/rnames = nms/clabels "Data" "Boot" "Bias" "SE"/format f9.4.
compute lowperc = boot(blowp,:).
compute upperc = boot(bhighp,:).
compute ci = {lowbca, upbca}.
do if (!bca <> 0).
print ci/title = "Bias Corrected and Accelerated Confidence Intervals"/rnames = nms/clabels "Lower" "Upper"/format F9.4.
end if.
do if (!bc <> 0).
compute ci = {lowbc, upbc}.
print ci/title = "Bias Corrected Confidence Intervals"/rnames = nms/
clabels "Lower" "Upper"/format F9.4.
end if.
do if (!percent <> 0).
compute ci = {t(lowperc), t(upperc)}.
print ci/title = "Percentile Confidence Intervals"/rnames = nms/
clabels "Lower" "Upper"/format F9.4.
end if.
print conf/title = "Level of Confidence for Confidence Intervals:".
print btn/title = "Number of Bootstrap Resamples:".
END MATRIX.
RESTORE.
!END DEFINE.
SAS (v1.0,
%macro
indirect(data=,y=,x=,m=,c=0,boot=100,conf=95,percent=0,bc=0,bca=1);
proc
iml;
use
&data;
read
all var{&y &x &m} into dd;
nm={&y
&x &m};
xx=(dd
= .);xx=xx[,+];
j=1;do
i = 1 to nrow(dd);if xx[i,1]=0 then;do;dd[j,]=dd[i,];j=j+1;end;end;
dd=dd[1:j-1,];
nm
= nm`;
n
= nrow(dd);
nv
= ncol(dd);
nc
= &c;
con=j(n,1,1);
dt2
= dd;
dt
= dd;
bzx
= j(nv-2-nc,1,0);
bzxse
= j(nv-2-nc,1,0);
b=j((nv-1-nc),(nv-1-nc),0);
p0
= -0.322232431088;
p1
= -1;
p2
= -0.342242088547;
p3
= -0.0204231210245;
p4
= -.0000453642210148;
q0
= 0.0993484626060;
q1
= 0.588581570495;
q2
= 0.531103462366;
q3
= 0.103537752850;
q4
= 0.0038560700634;
conf=round(&conf);
lowalp
= 0.5*(1-(conf/100));
upalp
= 0.5*(1+(conf/100));
zbca
= lowalp//upalp;
btn
= 1000;
if
(&boot > 999) then;
do;
btn = floor(&boot/1000)*1000;
end;
blowp
= floor(lowalp*btn);
if
(blowp < 1) then;
do;
blowp = 1;
end;
bhighp
= floor((upalp*btn)+1);
if
(bhighp > btn) then;
do;
bhighp = btn;
end;
indeff
= j((n+1+btn),(nv-1-nc),0);
do
d = 1 to (n+1+btn);
if (d = (n+2)) then;
do;
dt = dt2;
con = j(n,1,1);
end;
if (d > 1) then if (d < (n+2))
then;
do;
if (d = 2) then;
do;
con = j((n-1),1,1);
dt = dt2[2:n,];
end;
if (d = (n+1)) then;
do;
dt = dt2[1:(n-1),];
end;
if (d > 2) then if (d < (n+1))
then;
do;
dt = dt2[1:(d-2),]//dt2[(d:n),];
end;
end;
if (d > (n+1)) then;
do;
do nn = 1 to n;
v = int(ranuni(0)*n)+1;
dt[nn,1:nv]=dt2[v,1:nv];
end;
end;
x
= dt[,2];
m
= dt[,3:(nv-nc)];
y
= dt[,1];
xz
= dt[,2:nv];
xo
= con||x;
if
(nc > 0) then;
do;
c = dt[,(nv-nc+1):nv];
xo = xo||c;
end;
do
k = 3 to (nv-nc);
ytmp = dt[,k];
bzxt = inv(xo`*xo)*xo`*ytmp;
bzx[(k-2),1]=bzxt[2,1];
if (d = 1) then;
do;
mse = sum((ytmp-(xo*bzxt))##2)/(n-2-nc);
olscm = (mse*inv(xo`*xo));
bzxse[(k-2),1]=sqrt(olscm[2,2]);
end;
end;
if
(d = 1) then;
do;
if (nc > 0) then;
do;
cnt = dd[,(
xo = con||x||cnt;
end;
if (nc = 0) then;
do;
xo = con||x;
end;
byx = inv(xo`*xo)*xo`*y;
mse = sum((y-(xo*byx))##2)/(n-2-nc);
olscm = (mse*inv(xo`*xo));
byxse = sqrt(olscm[2,2]);
byx = byx[2,1];
end;
xzo
= con||xz;
byzx
= inv(xzo`*xzo)*xzo`*y;
byzx2
= byzx[3:(nv-nc),1];
if
(d = 1) then;
do;
mse = sum((y-(xzo*byzx))##2)/(n-nv);
covmat = mse*inv(xzo`*xzo);
olscm = vecdiag(covmat);
sse = mse*(n-nv);
sst = sum((y-(sum(y)/n))##2);
r2 = 1-(sse/sst);
ar2 = 1-(mse/(sst/(n-1)));
fr = ((n-nv)*r2)/((1-r2)*ncol(xz));
pfr = 1-probf(fr,ncol(xz),(n-nv));
if (nc > 0) then;
do;
bcon = byzx[(nv-nc+1):nv,1];
bconse = sqrt(olscm[(nv-nc+1):nv,1]);
end;
byzx2se = sqrt(olscm[3:(nv-nc),1]);
cprime = byzx[2,1];
cprimese=sqrt(olscm[2,1]);
end;
indeff2
= (bzx#byzx2);
zs
= (bzx/bzxse)#(byzx2/byzx2se);
temp
= t(sum(indeff2)//indeff2);
indeff[d,]=temp;
if
(d = 1) then;
do;
vs = nm[1:(nv-nc),1];
rn = {"DV = " "IV = " "MEDS
= "};
print "Dependent, Independent, and
Proposed Mediator Variables";
print vs [rowname = rn];
if (nc > 0) then;
do;
vs = nm[(nv-nc+1):nv,1];
print "Statistical Controls";
rn = {"CONTROLS="};
print vs [rowname = rn];
end;
nms = nm[3:(nv-nc),1];
te = bzx/bzxse;
df = n-2-nc;
p = 2*(1-probt(abs(te),df));
bzxmat = bzx||bzxse||te||p;
b[2:(nv-1-nc),1]=bzx;
se2 = bzxse#bzxse;
cnm = {"Coeff" "se" "t"
"p"};
print "IV to Mediators (a paths)";
print bzxmat [rowname = nms colname = cnm
format = 9.4];
te = byzx2/byzx2se;
df = n-nv;
p = 2*(1-probt(abs(te),df));
byzx2mat=byzx2||byzx2se||te||p;
print "Direct Effects of Mediators on DV
(b paths)";
print byzx2mat [rowname = nms colname = cnm
format = 9.4];
te=byx/byxse;
df = n-2-nc;
p=2*(1-probt(abs(te),df));
byxmat = byx||byxse||te||p;
xnm=nm[2,1];
print "Total effect of IV on DV (c
path)";
print byxmat [rowname = xnm colname = cnm
format = 9.4];
te=cprime/cprimese;
df = n-nv;
p=2*(1-probt(abs(te),df));
cprimmat = cprime||cprimese||te||p;
print "Direct Effect of IV on DV (c'
path)";
print cprimmat [rowname = xnm colname = cnm
format = 9.4];
if (nc > 0) then;
do;
df = n-nv;
nms = nm[(nv-nc+1):nv,1];
te=bcon/bconse;
p=2*(1-probt(abs(te),df));
bconmat = bcon||bconse||te||p;
print "Partial Effect of Control
Variables on DV";
print bconmat [rowname = nms colname =
cnm format = 9.4];
end;
dvms=r2||ar2||fr||ncol(xz)||(n-nv)||pfr;
print "Fit Statistics for DV Model";
cnm = {"R-sq" "adj R-sq" "F"
"df1" "df2" "p"};
print dvms [colname = cnm format = 9.4];
end;
end;
lvout
= indeff[2:(n+1),];
tdotm
= lvout[+,]/n;
tm
= j(n,ncol(lvout),1)*diag(tdotm);
topa=(((n-1)/n)*(tm-lvout))##3;
topa=topa[+,];
bota
=((((n-1)/n)*(tm-lvout))##2);
bota=bota[+,];
bota=6*sqrt(bota##3);
ahat
= topa/bota;
nms
= {"TOTAL"}//nm[3:(nv-nc),1];
indsam
= indeff[1,]`;
boot
= indeff[(n+2):nrow(indeff),];
mnboot
= (boot[+,]/btn)`;
xt=boot-j(btn,1)*boot[:,];
cv=(xt`*xt)/btn;
se=sqrt(vecdiag(cv));
if
(btn > 1) then;
do;
create bootstp from boot [colname='indirect'];
append from boot;
nnn = j(1,(nv-1-nc),-999);
boot = nnn//boot;
do e = 1 to (nv-1-nc);
do i = 2 to (btn+1);
ix = boot[i,e];
do k = i to 2 by -1;
k2 = k;
if (boot[(k-1),e] > ix) then;
do;
boot[k,e]=boot[(k-1),e];
end;
else;
if (boot[(k-1),e] <= ix) then;
do;
goto stpit;
end;
end;
stpit:
boot[k2,e]=ix;
end;
end;
boot = boot[2:(btn+1),];
end;
xp=j((nrow(mnboot)+2),1,0);
do
i = 1 to (nrow(mnboot)+2);
if (i <= nrow(mnboot)) then;
do;
pv = (boot[,i] < indsam[i,1]);
pv = pv[+,]/btn;
end;
else;
pv = zbca[(i-nrow(mnboot)),1];
p=pv;
if (pv > 0.5) then;
do;
p = 1-pv;
end;
y5 = sqrt(-2*log(p));
xp[i,1]=y5+((((y5*p4+p3)*y5+p2)*y5+p1)*y5+p0)/((((y5*q4+q3)*y5+q2)*y5+q1)*y5+q0);
if (pv <= 0.5) then;
do;
xp[i,1]=-xp[i,1];
end;
end;
bbb
= nrow(mnboot);
zz
= xp[1:bbb,1];
zlo
= zz + ((zz+xp[(bbb+1),1])/(1-ahat`#(zz+xp[(bbb+1),1])));
zup
= zz + ((zz+xp[(bbb+2),1])/(1-ahat`#(zz+xp[(bbb+2),1])));
ahat
= 0;
zlobc
= zz + ((zz+xp[(bbb+1),1])/(1-ahat`#(zz+xp[(bbb+1),1])));
zupbc
= zz + ((zz+xp[(bbb+2),1])/(1-ahat`#(zz+xp[(bbb+2),1])));
zlo
= probnorm(zlo);
zup
= probnorm(zup);
zlobc
= probnorm(zlobc);
zupbc
= probnorm(zupbc);
blow
= int(zlo*(btn+1));
bhigh
= int(zup*(btn+1))+1;
blowbc
= int(zlobc*(btn+1));
bhighbc
= int(zupbc*(btn+1))+1;
lowbca
= j(nrow(blow),1,0);
upbca
= lowbca;
do
i = 1 to nrow(blow);
if (blow[i,1] < 1) then;
do;
blow[i,1]=1;
end;
lowbca[i,1]=boot[blow[i,1],i];
if (bhigh[i,1] > btn) then;
do;
bhigh[i,1]=btn;
end;
upbca[i,1]=boot[bhigh[i,1],i];
end;
lowbc
= j(nrow(blow),1,0);
upbc
= lowbca;
do
i = 1 to nrow(blowbc);
if (blowbc[i,1] < 1) then;
do;
blowbc[i,1]=1;
end;
lowbc[i,1]=boot[blowbc[i,1],i];
if (bhighbc[i,1] > btn) then;
do;
bhighbc[i,1]=btn;
end;
upbc[i,1]=boot[bhighbc[i,1],i];
end;
print
"*****************************************************";
print
"BOOTSTRAP RESULTS FOR INDIRECT EFFECTS";
res
= indsam||mnboot||(mnboot-indsam)||se;
cn
= {"Data" "Boot" "Bias" "SE"};
print
"Total and Specific Indirect Effect(s) of IV on DV via MEDS";
print
res [rowname = nms colname = cn format = 9.4];
lowperc
= boot[blowp,];
upperc
= boot[bhighp,];
ci
= lowbca||upbca;
cn
= {"Lower" "Upper"};
if
(&bca ^= 0) then;
do;
print "Bias Corrected and Accelerated
Confidence Intervals";
print ci [rowname = nms colname = cn format =
9.4];
end;
if
(&bc ^= 0) then;
do;
ci = lowbc||upbc;
print "Bias Corrected Confidence
Intervals";
print ci [rowname = nms colname = cn format =
9.4];
end;
if (&percent ^= 0) then;
do;
ci = lowperc`||upperc`;
print "Percentile Confidence
Intervals";
print ci [rowname = nms colname = cn format =
9.4];
end;
print
"Level of Confidence for Confidence Intervals";
print
conf;
print
"Number of Bootstrap Resamples";
print
btn;
quit;
%mend;
