%MACRO INDTEST(p1=2,p2=4,p3=5,p4=8,draws=100000,seed=-1,alpha=.05);
*Macro requires that this ods statement be added to proc mixed with 
 the designated file names:
    ods output covb=acovfix asycov=acovrand solutionf=estfix 
               covparms=estrand convergencestatus=converge;

*Macro arguments;
  *p1 is location of a in vector of fixed effects estimates, 
    that is, the row in Solution for Fixed Effects table of Mixed output;
  *p2 is location of b in vector of fixed effects estimates (as above);
  *p3 is location of c prime in vector of fixed effects estimates (as above);
  *p4 is location of COV(aj,bj) within the covariance parameter output, 
    that is, the row in Covariance Parameter Estimates table of Mixed output;
  *draws is the number of random draws to create the simulated sampling
   distribution for the Monte Carlo confidence intervals;
  *seed is the random number need to generate simulated sampling 
   distribution for the Monte Carlo confidence intervals;
  *alpha is the desired Type I error rate;
  *NOTE draws*(alpha/2) must equal a whole number;

PROC IML;

USE CONVERGE;
READ ALL VAR{STATUS} INTO STATUS;
IF STATUS = 0 THEN DO;
*Collecting necessary model estimates from fixed effects and G matrix;
USE ESTFIX;
READ ALL VAR {EFFECT} INTO PARAMS;
READ ALL VAR {ESTIMATE} INTO EST;
a = EST[&p1.];
b = EST[&p2.];
c = EST[&p3.];

USE ESTRAND;
READ ALL VAR {CovParm} INTO PARAMS2;
READ ALL VAR {ESTIMATE} INTO EST;
sigma_ajbj = EST[&p4.];

Estimates = a//b//c//sigma_ajbj;

*Collecting necessary elements from asymptotic covariance matrices; 
USE acovfix;
READ ALL VAR _NUM_ INTO COVEST;
COVEST = COVEST[,2:NCOL(COVEST)];
COV_ABC = COVEST[&p1.||&p2.||&p3., &p1.||&p2.||&p3.];

USE acovrand;
READ ALL VAR _NUM_ INTO COVEST;
COVEST = COVEST[,2:NCOL(COVEST)];
COV_Sigma = COVEST[&p4.,&p4.]; *covariance of Sigma_ajbj;

Covariances = BLOCK(COV_ABC,COV_Sigma);

names = {"a" "b" "c prime" "cov(aj,bj)"};
PRINT Estimates[label="Estimates of a, b, c prime, and cov(aj,bj)" rownames=names];
PRINT Covariances[label="Asymptotic covariance matrix of a, b, c prime, and cov(aj,bj)"
                  rownames=names colnames=names];

/****************************************************************************************/
/*Computing analytical estimate of the average indirect effect and its sampling variance*/
/****************************************************************************************/

Eab = Estimates[1]*Estimates[2] + Estimates[4];
Va = Covariances[1,1];
Vb = Covariances[2,2];
Cab = Covariances[2,1];
VSajbj = Covariances[4,4];
Vab = (a**2)*Vb + (b**2)*Va + Va*Vb + 2*a*b*Cab + Cab**2 + VSajbj;
CR = Eab/SQRT(Vab);
P = (1-(PROBNORM(ABS(CR))))*2; *two-tailed p-value based on standard normal dist;
LCL = Eab - probit(1 - &alpha./2)*SQRT(Vab); *95% Confidence limits;
UCL = Eab + probit(1 - &alpha./2)*SQRT(Vab);
Results = Eab||(SQRT(Vab))||CR||P||&alpha.||LCL||UCL;
labels={"Estimated average indirect effect",
        "Standard error",
		"Critical Ratio",
		"P-Value", "alpha level",
		"LCL", "UCL"};
PRINT (Results`)[label="Analytical Estimate of Average Indirect Effect and Standard Error" rowname=labels]; 

*computing average total effect;
c = estimates[3];
Etot = Eab + c;
Vc = Covariances[3,3];
Cac = Covariances[3,1];
Cbc = Covariances[3,2];
Vtot = Vab + Vc + 2*b*Cac + 2*a*Cbc;
CR = Etot/SQRT(Vtot);
P = (1-(PROBNORM(ABS(CR))))*2;
LCL = Etot - probit(1-&alpha./2)*SQRT(Vtot);
UCL = Etot + probit(1-&alpha./2)*SQRT(Vtot);
Results = Etot||(SQRT(Vtot))||CR||P||&alpha.||LCL||UCL;
labels={"Estimated average total effect",
        "Standard error",
		"Critical Ratio",
		"P-Value", "alpha level",
		"LCL", "UCL"};
PRINT (Results`)[label="Analytical Estimate of Average Total Effect and Standard Error" rowname=labels]; 


/***********************************************************************/
/*Simulating distribution of average indirect and total effect for     */
/*asymmetric CI                                                        */
/***********************************************************************/

*drawing random normal variates a, b, c, sigma(aj,bj);
Z = rannor(J(&draws.,4,&seed.)); *starts with standard normal;
L = root(covariances);
Sim_Est = Z*L + J(&draws.,1,1)*Estimates`;
Sim_Ind = Sim_Est[,1]#Sim_Est[,2]+Sim_Est[,4];
Sim_Tot = Sim_Ind + Sim_Est[,3];

Sim = Sim_Ind||Sim_Tot;

CREATE Sim FROM Sim[COLNAME={Ind Tot}];
APPEND FROM Sim;

END;

*If failed to converge, then...;
ELSE DO;
hold={-999 -999 -999 -999 -999 -999};
CREATE Results from hold[COLNAME={Est SE CR p LCL UCL}];
APPEND from hold;
hold = {-999};
CREATE SIM from hold[COLNAME={Eab}];
APPEND FROM hold;
END;
QUIT;

proc sort data=Sim; by Ind; run;
data rank; set Sim;
  if _N_ = (&alpha./2)*&draws. then output;
  if _N_ = (1 - &alpha./2)*&draws. then output;
run;
proc transpose data=rank out=ind_centiles prefix=CL; var Ind; run;

proc sort data=Sim; by Tot; run;
data rank; set Sim;
  if _N_ = (&alpha./2)*&draws. then output;
  if _N_ = (1 - &alpha./2)*&draws. then output;
run;
proc transpose data=rank out=tot_centiles prefix=CL; var Tot; run;

data centiles; set ind_centiles tot_centiles;
 if _N_ = 1 then effect = "average indirect effect";
 if _N_ = 2 then effect = "average total effect";
 rename CL1=LCL_sim CL2=UCL_sim;
 alpha = &alpha.;
 draws = &draws.;
 drop=_NAME_;
run;


proc print data=centiles noobs; 
 title 'Monte Carlo confidence interval';
 var effect alpha draws LCL_sim UCL_sim; 
run;

title;

%MEND;