Hello I wish to call the following macro and run the analysis below. The macro runs well but the last part to simulate a data and runs does not work. Any help?
The problem is from the data test (marked bold below) , no problem with the macro. The error is at the tail end
%macro indirect(data=,y=,x=,m=,c=0,boot=1000,conf=95,percent=0,bc=1,bca=0, normal=0,contrast=0);
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;
resid = j(n,(nv-nc),0);
info = j((2*(nv-nc-2)+1),(2*(nv-nc-2)+1),0);
imat = j(ncol(info),4,1);
imat[1:(nv-nc-2),1]=(2:(nv-nc-1))`;
imat[1:(nv-nc-2),3]=(2:(nv-nc-1))`;
imat[(nv-nc-1):(ncol(info)-1),2]=(2:(nv-nc-1))`;
imat[(nv-nc-1):(ncol(info)-1),4]=(2:(nv-nc-1))`;
imat[(nv-nc-1):(ncol(info)-1),1]=j((nv-nc-2),1,(nv-nc));
imat[(nv-nc-1):(ncol(info)-1),3]=j((nv-nc-2),1,(nv-nc));
imat[ncol(info),1:4]={1 1 1 1};
imat[ncol(info),1]=nv-nc;
imat[ncol(info),3]=nv-nc;
bzx = j(nv-2-nc,1,0);
bzxse = j(nv-2-nc,1,0);
b=j((nv-1-nc),(nv-1-nc),0);
cname={"C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8", "C9", "C10", "C11", "C12", "C13", "C14", "C15", "C16"};
cname=cname//{"C17", "C18", "C19", "C20", "C21", "C22", "C23", "C24", "C25", "C26", "C27", "C28", "C29"};
cname=cname//{"C30", "C31", "C32", "C33", "C34", "C35", "C36", "C37", "C38", "C39", "C40", "C41", "C42"};
cname=cname//{"C43", "C44", "C45"};
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 = 1;
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;
resid[,(k-1)]=ytmp-(xo*bzxt);
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[,(nv-(nc-1)):nv];
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;
resid[,ncol(resid)]=y-(xzo*byzx);
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;
print "Sample size";
print n;
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];
if (&normal =^ 0) then;if(&c = 0) then;do;
do;
bmat=j((nv-nc),(nv-nc),0);
bmat[2:(nv-nc-1),1]=bzx;
bmat[(nv-nc),2:(nv-nc-1)]=byzx2`;
bmat[(nv-nc),1]=cprime;
imbinv = inv(i(ncol(bmat))-bmat);
imbtinv = inv(i(ncol(bmat))-bmat`);
resid[,1]=x-(x[+]/n);
psi = (resid`*resid)/(n-1);
invpsi = inv(psi);
ibpsiib = imbinv*psi*imbtinv;
do ic = 1 to ncol(info);
do ic2 = 1 to ncol(info);
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;
end;
varcov=inv(info);
varcov=varcov[1:(2*(nv-nc-2)),1:(2*(nv-nc-2))];
ses = vecdiag(varcov);
avar = ses[1:nrow(bzxse),1];
bvar = ses[(nrow(bzxse)+1):nrow(ses),1];
if ((nv-nc-2) > 1) then;do;if(&contrast = 1) then;do;
prws=j(((nv-nc-2)*(nv-nc-3)/2),1,0);
prwse=prws;
kk=1;
do ic = 1 to (nv-nc-3);
do ic2 = (ic+1) to (nv-nc-2);
vf2=((byzx2[ic,1]##2)*varcov[ic,ic])-(2*byzx2[ic,1]*byzx2[ic2,1]*(varcov[ic,ic2]));
vf2=vf2+((byzx2[ic2,1]##2)*varcov[ic2,ic2])+((bzx[ic,1]##2)*(bvar[ic,1]));
vf2=vf2-(2*bzx[ic,1]*bzx[ic2,1]*covmat[(2+ic),(2+ic2)])+((bzx[ic2,1]##2)*(bvar[ic2,1]));
cnt=indeff2[ic,1]-indeff2[ic2,1];
prws[kk,1]=cnt;
prwse[kk,1]=sqrt(vf2);
kk=kk+1;
end;
end;
cnam2=cname[1:(kk-1),1];
end;end;
dermat=byzx2//bzx;
totse=sqrt((dermat)`*varcov*dermat);
specse = sqrt((byzx2#byzx2)#(avar)+(bzx#bzx)#(bvar));
specse = totse//specse;
indsum=indeff2[+];
specz = (indsum//indeff2)/specse;
ind22 = indsum//indeff2;
nms = {"TOTAL"}//nm[3:(nv-nc),1];
if ((nv-nc-2) > 1) then;do;if(&contrast = 1) then;do;
ind22 = ind22//prws;
specse = specse//prwse;
specz2 = prws/prwse;
specz = specz//specz2;
nms = nms//cnam2;
end;end;
pspec = 2*(1-probnorm(abs(specz)));
spec = ind22||specse||specz||pspec;
cnm = {"Effect" "se" "Z" "p"};
print "*****************************************************";
print "NORMAL THEORY TESTS FOR INDIRECT EFFECTS";
print "Indirect Effects of IV on DV through Mediators (ab paths)";
print spec [rowname = nms colname = cnm format = 9.4];
end;end;
end;
end;
if (btn > 1) then;do;
nms = {"TOTAL"}//nm[3:(nv-nc),1];
if (nv-nc-2) > 1 then do;if (&contrast = 1) then do;
crst=j((n+1+btn),((nv-nc-2)*(nv-nc-3)/2),0);
kk=1;
do ic = 2 to (nv-nc-2);
do ic2 = (ic+1) to (nv-nc-1);
crst[,kk]=indeff[,ic]-indeff[,ic2];
kk=kk+1;
end;
end;
indeff = indeff||crst;
cnam2=cname[1:(kk-1),1];
nms = nms//cnam2;
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;
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));
create bootstp from boot [colname='indirect'];
append from boot;
nnn = j(1,ncol(indeff),-999);
boot = nnn//boot;
do e = 1 to (ncol(indeff));
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),];
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 "Indirect Effects of IV on DV through Mediators (ab paths)";
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 "*****************************************************";
print "Level of Confidence for Confidence Intervals";
print conf;
print "Number of Bootstrap Resamples";
print btn;
end;
if (&normal = 1) then;do;if (&c = 0) then;do;prt =1;end;end;
if (btn > 999) then;do;prt = 1;end;
if ((nv-nc-2) > 1) then;do;if (&contrast = 1) then;do;
if (prt = 1) then;do;
print "*****************************************************";
print "Indirect Effect Contrast Definitions: IndEff_1 minus IndEff2";
kk=1;
prwsv = j(((nv-nc-2)*(nv-nc-3)/2),2,"XXXXXXXX");
do ic = 1 to (nv-nc-3);
do ic2 = (ic+1) to (nv-nc-2);
prwsv[kk,1]=nm[ic+2,1];
prwsv[kk,2]=nm[ic2+2,1];
kk=kk+1;
end;
end;
prwsv = cnam2||prwsv;
cn = {"Contrast" "IndEff_1" "IndEff_2"};
print prwsv [colname = cn];
end;end;end;
quit;
%mend;
/* Macro ends*/;
/*Call macro into data set*/;
data test;
seed=-1;
do i = 1 to 100;
Y=10+(.5)*rannor(1);
X=2+(0.3)*rannor(1);
/*** These are the mediators *****/;
M1=1+(0.3)*rannor(1);
M2=3+(0.4)*rannor(1);
M2=5+(0.8)*rannor(1);
M4=0+(.9)*rannor(1);
/***** These are the covariates *******/;
gender=ranbin(0,1,.5);
age=60+(5)*rannor(1);
edu=12+(1)*rannor(1);
output;
end;
run;
/*So i want to call the macro indirect and fit this model, can any one help with the last step*/
%INDIRECT [DATA = test, Y = yvar, X = xvar, M = M1 M2 M3 M4{gender age edu}]; /* Note M=Mvlist are M1, M2, M3 and M4 and covlist are age gender and edu */
{,C = {cov}(0**)}
{,BOOT = {z}(1000**)}
{,CONF = {ci}(95**)}
{,NORMAL = {t}(0**)}
{,CONTRAST = {n}(0**)}
{,PERCENT = {p}(0**)}
{,BC = {b}(1**)}
{,BCA = {d}(0)};
ERROR:
ERROR
982
983 /*So i want to call the macro indirect and fit this model, can any one help with the last
983! step*/
984 %INDIRECT [DATA = test, Y = y, X = x, M = M1 M2 M3 M4{gender age edu}]; /* Note M=Mvlist
984! are M1, M2, M3 and M4 and covlist are age gender and edu */
NOTE: IML Ready
NOTE: Line generated by the invoked macro "INDIRECT".
1 use &data;
-
