/* This macro estimates total indirect and specific indirect effects, */.
/* and bootstrap confidence intervals in single-step multiple mediator models.  The syntax is */.
/* */.
/* INDIRECT y = dv/x = iv/m = mlist covlist/c = cov/boot = z/conf = ci/percent = p/bc = b/bca = d */.
/* */.
/* where dv is the dependent variable, iv is the independent variable, and mlist is a */.
/* list of one or more mediator variables through with the IV's effect is presumably transmitted to the */.
/* dv, covlist is a list of covariates, cov is the number of covariates in covlist, z is the number of bootstrap */.
/* resamples desired in increments of 1000, ci is the desired confidence for confidence intervals (1 to 99), */.
/* p is set to 1 to print percentile confidence intervals, b is set to 1 for bias-corrected confidence */.
/* intervals, and d is set to 0 to disable printing of bias-corrected and adjusted confidence intervals */.
/* c defaults to 0 (meaning covlist is a null string), z defaults to 1000, ci default to 95. Bias-corrected and */.
/* accelerated confidence intervals are printed by default.  bc and bca confidence intervals are not printed by default */.
/* Anything after mlist in the syntax is optional */.
/* The macro assumes missing values are represented in the data with a period.  User defined missing */.
/* values will be treated as valid data, so user defined missing data should be purged from the */.
/* data prior to running the macro */.
/* */.
/* Macro written by Andrew F. Hayes, School of Communication, The Ohio State University, hayes.338@osu.edu */.

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(:,(nv-(nc-1)):nv)).
      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((nv-!c+1):nv,1).
      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.