/* This macro estimates total indirect and specific indirect effects, */.
/* and bootstrap confidence intervals in single-step mediator models.  The syntax is */.
/* */.
/* INDIRECT y = dv/x = iv/m = mlist covlist/c = cov/boot = z/conf = ci/normal = n/contrast = t/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), */.
/* n is set to 1 to print normal theory standard errors for indirect effects, */.
/* t is set to 1 to do all possible pairwise contrasts between indirect effects, */.
/* 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 */.
/* cov defaults to 0 (meaning covlist is a null string), t defaults to 0, z defaults to 1000, ci default to 95, n defaults to 0 */.
/* 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 */.
/* version 2 Dec 13 2005 */.

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)/normal = !charend ('/') !default(0)/contrast = !charend ('/') !default(0)).
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 resid = make(n,(nv-nc),0).
compute info = make((2*(nv-nc-2)+1),(2*(nv-nc-2)+1),0).
compute imat = make(ncol(info),4,1).
compute imat(1:(nv-nc-2),1)=t({2:(nv-nc-1):1}).
compute imat(1:(nv-nc-2),3)=t({2:(nv-nc-1):1}).
compute imat((nv-nc-1):(ncol(info)-1),2)=t({2:(nv-nc-1):1}).
compute imat((nv-nc-1):(ncol(info)-1),4)=t({2:(nv-nc-1):1}).
compute imat((nv-nc-1):(ncol(info)-1),1)=make((nv-nc-2),1,(nv-nc)).
compute imat((nv-nc-1):(ncol(info)-1),3)=make((nv-nc-2),1,(nv-nc)).
compute imat(ncol(info),:)={(nv-nc),1,(nv-nc),1}.
compute cname={"C1";"C2";"C3";"C4";"C5";"C6";"C7";"C8";"C9";"C10";"C11";"C12";"C13";"C14";"C15";"C16";"C17"}.
compute cname={cname;"C18";"C19";"C20";"C21";"C22";"C23";"C24";"C25";"C26";"C27";"C28";"C29";"C30";"C31"}.
compute cname={cname;"C32";"C33";"C34";"C35";"C36";"C37";"C38";"C39";"C40";"C41";"C42";"C43";"C44";"C45"}.
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 = 1.
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 resid(:,#k-1) = ytmp-(xo*bzxt).
       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 resid(:,ncol(resid))=y-(xzo*byzx).
    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.
    print n/title = "Sample size"/format F10.0.
    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.
    do if (!normal <> 0 and !c = 0).
    compute bmat = make((nv-nc),(nv-nc),0).
    compute bmat(2:(nv-nc-1),1) = bzx.
    compute bmat((nv-nc),2:(nv-nc-1))=t(byzx2).
    compute bmat((nv-nc),1) = cprime.
    compute imbinv = inv(ident(ncol(bmat))-bmat).
    compute imbtinv=inv(ident(ncol(bmat))-t(bmat)).
    compute resid(:,1)=x-(csum(x)/(n)).
    compute psi = sscp(resid)/(n).
    compute invpsi = inv(psi).
    compute ibpsiib = imbinv*psi*imbtinv.
    loop ic = 1 to ncol(info).
    loop ic2 = 1 to ncol(info).
    compute info(ic,ic2)=(n-1)*((imbinv(imat(ic2,4),imat(ic,1))*imbinv(imat(ic,2),imat(ic2,3)))+(ibpsiib(imat(ic2,4),imat(ic,2))*invpsi(imat(ic,1),imat(ic2,3)))).
    end loop.
    end loop.
    compute varcov = inv(info). 
    compute varcov = varcov(1:(2*(nv-nc-2)),1:(2*(nv-nc-2))).
    compute ses = diag(varcov).
    compute avar = ses(1:nrow(bzxse),1).
    compute bvar = ses((nrow(bzxse)+1):nrow(ses),1).
    do if ((nv-nc-2) > 1 and (!contrast = 1)).
      compute prws=make(((nv-nc-2)*(nv-nc-3)/2),1,0).
      compute prwse=prws.
      compute kk=1.
      loop ic = 1 to (nv-nc-3).
        loop ic2 = (ic+1) to (nv-nc-2).
          compute vf2 = ((byzx2(ic,1)**2)*varcov(ic,ic))-(2*byzx2(ic,1)*byzx2(ic2,1)*(varcov(ic,ic2))).
          compute vf2=vf2+((byzx2(ic2,1)**2)*varcov(ic2,ic2))+((bzx(ic,1)**2)*(bvar(ic,1))).
          compute vf2=vf2-(2*bzx(ic,1)*bzx(ic2,1)*covmat((2+ic),(2+ic2)))+((bzx(ic2,1)**2)*(bvar(ic2,1))).
          compute cnt = indeff2(ic,1)-indeff2(ic2,1).
          compute prws(kk,1)=cnt.
          compute prwse(kk,1)=sqrt(vf2).
          compute kk=kk+1.
        end loop.
      end loop.
    compute cnam2 = cname(1:(kk-1),1).
    end if.
    compute dermat = {byzx2;bzx}.
    compute totse = sqrt(t(dermat)*varcov*dermat).
    compute specse = sqrt((byzx2&*byzx2)&*(avar)+(bzx&*bzx)&*(bvar)).
    compute specse = {totse; specse}.       
    compute specz = {csum(indeff2);indeff2}&/specse.
    compute ind22 = {csum(indeff2);indeff2}.
    compute nms = {"TOTAL";nm(3:(nv-nc),1)}.
    do if ((nv-nc-2) > 1 and (!contrast = 1)).
       compute ind22 = {ind22;prws}.
       compute specse = {specse;prwse}.
       compute specz = {specz;(prws&/prwse)}.
       compute nms = {nms;cnam2}.
    end if.
    compute pspec= 2*(1-cdfnorm(abs(specz))).
    compute spec = {ind22, specse, specz, pspec}.

    
    print/title = "******************************************************************".
    print/title = "           NORMAL THEORY TESTS FOR INDIRECT EFFECTS".
    print spec/title = "Indirect Effects of IV on DV through Proposed Mediators (ab paths)"/rnames = nms/clabels "Effect" "se" "Z" "p"/format = f9.4.
    end if.
  end if.
end loop.

do if (btn > 1).
  compute nms = {"TOTAL"; nm(3:(nv-nc),1)}.
  do if ((nv-nc-2) > 1 and (!contrast = 1)).
    compute crst = make((n+1+btn),((nv-nc-2)*(nv-nc-3)/2),0).
    compute kk=1.
    loop ic = 2 to (nv-nc-2).
      loop ic2 = (ic+1) to (nv-nc-1).
        compute crst(:,kk)=indeff(:,ic)-indeff(:,ic2).
        compute kk=kk+1.
      end loop.
    end loop.
    compute indeff = {indeff,crst}.
    compute cnam2 = cname(1:(kk-1),1).
    compute nms = {nms;cnam2}.
  end if.
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 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))).
save boot/outfile = indirect.sav/names = nms.
compute nnn = make(1,ncol(indeff),-999).
compute boot = {nnn;boot}.
loop #e = 1 to ncol(indeff).
  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),:).
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 = "Indirect Effects of IV on DV through Proposed Mediators (ab paths)"/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/title = "*****************************************************************".
print conf/title = "Level of Confidence for Confidence Intervals:".
print btn/title = "Number of Bootstrap Resamples:".
end if.
do if ((nv-nc-2) > 1 and (!contrast = 1) and ((!normal = 1 and !c = 0) OR btn > 999))).
print/title = "*****************************************************************".
print/title = "  INDIRECT EFFECT CONTRAST DEFINITIONS: Ind_Eff1 MINUS Ind_Eff2".
compute kk=1.
compute prwsv=make(((nv-nc-2)*(nv-nc-3)/2),2,0).
 loop ic = 1 to (nv-nc-3).
        loop ic2 = (ic+1) to (nv-nc-2).
          compute prwsv(kk,1)=nm(ic+2,1).
          compute prwsv(kk,2)=nm(ic2+2,1).
          compute kk=kk+1.
       end loop.
end loop.
compute prwsv = {cnam2, prwsv}.
print prwsv/title = " "/clabels = "Contrast" "IndEff_1" "IndEff_2"/format A9.
end if.
END MATRIX.
RESTORE.
!END DEFINE.