options ls=80 ps=60 formdlim='-';
ods rtf file='Insert correct path\Results_mediation_ch9.rtf'; 

/*Include macro for mediation analysis*/
%include 'Insert correct path\Mediation_macro.sas';

/*Create the modafinil dataset from Section 9.4.1. 
PctNeg is the outcome. 
N3_1 is duration of N3 sleep at baseline. 
N3_2 is duration of N3 sleep at time 2. 
Trt is 1 for modafinil and 0 for placebo. */
data sleep;
input id $ Trt PctNeg N3_1 N3_2;
cards;
T001	0	0.00	64.5	46.5
T002	1	56.25	55.0	79.0
T003	0	0.00	28.5	19.5
T004	0	14.29	103.0	80.0
T005	1	90.00	5.5	21.5
T006	0	.	1.5	.
T007	0	11.11	8.5	5.5
T008	0	.	48.5	.
T009	1	.	6.5	8.5
T010	1	33.33	7.0	45.5
T011	1	72.22	13.0	43.0
T012	1	0.00	63.5	70.5
T013	0	33.33	9.0	20.0
T014	0	36.36	86.0	94.0
T015	1	5.88	59.5	87.5
T016	1	.	18.0	36.5
T017	0	0.00	1.0	0.5
T018	0	56.25	90.5	130.5
T020	0	71.43	2.5	4.5
T021	0	18.18	83.0	89.5
T022	1	.	2.5	22.0
T025	1	30.00	15.5	52.5
T027	0	28.57	16.5	23.5
T028	0	15.38	7.5	25.5
T029	0	100.00	25.0	53.5
T030	1	.	3.0	.
T031	0	0.00	16.5	26.5
T032	0	0.00	0.5	0.5
T033	1	88.89	3.0	28.0
T034	1	88.24	88.0	109.5
T035	0	37.50	82.0	70.0
T036	0	0.00	39.0	7.5
T037	0	100.00	56.0	64.0
T039	1	52.63	55.0	49.5
T043	0	0.00	28.5	22.0
T044	1	64.71	54.5	70.0
T045	1	12.50	41.0	46.0
T046	1	33.33	0.5	0.0
T047	1	.	65.0	61.5
T049	1	100.00	41.5	61.0
T051	0	0.00	40.0	19.0
T053	1	10.53	74.0	55.0
T054	0	33.33	80.0	24.0
T055	0	11.11	12.0	5.5
T056	0	.	17.0	8.0
T057	1	.	13.5	29.5
T058	1	95.00	36.0	83.5
;

/*Create change in n3 sleep duration variable*/
data sleep;set sleep;
    n3sleep=n3_2-n3_1; 
run;

title 'Mediation analysis without an interaction in modafinil dataset';
%mediation(data=sleep,yvar=pctneg,avar=trt,mvar=n3sleep,cvar=,a0=0,a1=1,m=0,nc=0,yreg=linear,mreg=linear,
interaction=false,casecontrol=,output=full,c=,boot=false);

title 'Mediation analysis with an interaction in modafinil dataset';
%mediation(data=sleep,yvar=pctneg,avar=trt,mvar=n3sleep,cvar=,a0=0,a1=1,m=0,nc=0,yreg=linear,mreg=linear,
interaction=true,casecontrol=,output=,c=,boot=false);


/*Read in data for the intent to smoke example from Section 9.4.2.*/
data IntentLong;
   infile 'Insert correct path\IntentToSmoke_Ch9.txt';
   input intent treatment cig_use;
   id = _n_;
run;

title 'Table 9.3. Effects of a school-based drug prevention program on smoking';
proc freq data=IntentLong;
   table intent*treatment*cig_use;
run;

title 'Mediation analysis based on logistic model for the outcome in drug prevention program study';
%mediation(data=IntentLong,yvar=cig_use,avar=treatment,mvar=intent,cvar=,a0=0,a1=1,m=0,nc=0,yreg=logistic,mreg=logistic,
interaction=false,casecontrol=,output=full,c=,boot=false);

title 'Mediation analysis based on log-linear model for the outcome in drug prevention program study';
%mediation(data=IntentLong,yvar=cig_use,avar=treatment,mvar=intent,cvar=,a0=0,a1=1,m=0,nc=0,yreg=loglinear,mreg=logistic,
interaction=false,casecontrol=,output=full,c=,boot=false);


/*Interaction is not statistically significant so this code is not neededed.
Mediation analysis based on logistic model for the outcome without interaction
%mediation(data=IntentLong,yvar=cig_use,avar=treatment,mvar=intent,cvar=,a0=0,a1=1,m=0,nc=0,yreg=logistic,mreg=logistic,
interaction=true,casecontrol=,output=full,c=,boot=false);

Mediation analysis based on logistic model for the outcome with interaction
%mediation(data=IntentLong,yvar=cig_use,avar=treatment,mvar=intent,cvar=,a0=0,a1=1,m=0,nc=0,yreg=loglinear,mreg=logistic,
interaction=true,casecontrol=,output=full,c=,boot=false);*/


/*Read in the simulated dataset from Section 9.4.3. 
y is the output, m is the mediator, x the predictor, id the subject id. */
data simdata;
   infile 'Insert correct path\simdataCh9.txt';
   input id x m y;
run;

/*Reorganize the data so that analysis can be performed*/
data long;
    set simdata;
    Z = y; /*Combine the outcome and the mediator in one variable Z. Make Z equal to y*/
    Sy = 1; /*Create a selection variable Sy and set to 1 to indicate Z is Y*/
    Sm = 0; /*Create a selection variable Sm and set to 0 to indicate Z is not M*/
    dv = 'Y'; /*Create a variable dv to differentiate Y from M*/
    output; /*Output the first record for the observation*/
    Z = m; /*Make Z equal to M*/
    Sy = 0; /*Set the selection variable Sy to 0 to indicate Z is not Y*/
    Sm = 1; /*Set the selection variable Sm to 1 to indicate Z is M*/
    dv = 'M'; /*Create variable dv also to differentiate Y from M*/
    output; /*Output the second record for the observation*/
run;

title 'Joint mixed model analysis of the outcome on the mediator and predictor, and of the mediator on the predictor in simulated data set';
proc mixed data=long asycov covtest noclprint method=reml maxiter=1000;
    class id dv;
    model Z = Sm Sm*x Sy Sy*x Sy*m / noint solution ddfm=kr covb;
    random Sm Sm*x Sy Sy*x Sy*m /g gcorr type=un subject=id;
    repeated / group=dv subject=id;
    ods output covb=acovfix asycov=acovrand solutionf=estfix covparms=estrand convergencestatus=converge;
run;

/*Reorganize the data to calculate CDE, NDE, NIE and TIE*/
data results (keep=Effect Estimate); 
   set estfix;
   if Effect="Sm*x" or Effect="x*Sy" or Effect="Sy*m";
run;

proc sort data=results; by Estimate; run;
proc transpose data=results out=res (drop= _NAME_); 
   id Effect;
   var Estimate;
run;

data results2 (keep=CovParm Estimate); 
    set estrand;
    if CovParm="UN(5,2)";
run;

proc transpose data=results2 out=resbis (drop= _NAME_ ); 
   id CovParm;
   var Estimate;
run;

data res;
   merge res resbis;
   nde=x_Sy;
   nie=Sm_x*Sy_m + UN_5_2_;
   tce = nde + nie;
   pctmed = nie /tce;
run;

proc print data=res; 
   var nde nie tce pctmed;
run;

ods rtf close;

/*Use bootstrap to calculate standard errors*/
proc sort data=long nodupkey out=boots; by id; run;
proc sort data=boots; by id; run;

/* Bootstrap macro of Bind et al (2016)*/
%macro bootstrap;
%do i=1 %to 100;
proc surveyselect data=boots out=bootsample (keep=id)
method=urs samprate=1 outhits rep=1;
run;
ods listing close;

proc sql;
    create table bootstraped as select * from bootsample,long
    where bootsample.id=long.id;
quit;

proc mixed data=bootstraped asycov covtest noclprint method=reml maxiter=1000;
    class id dv;
    model Z = Sm Sm*x Sy Sy*x Sy*m / noint solution ddfm=kr covb;
    random Sm Sm*x Sy Sy*x Sy*m /g gcorr type=un subject=id;
    repeated / group=dv subject=id;
    ods output covb=acovfix asycov=acovrand solutionf=estfix covparms=estrand convergencestatus=converge;
run;

data results (keep=Effect Estimate); 
   set estfix;
   if Effect="Sm*x" or Effect="x*Sy" or Effect="Sy*m";
run;

proc sort data=results; by Estimate; run;
proc transpose data=results out=res2 (drop= _NAME_); 
   id Effect;
   var Estimate;
run;

data results2 (keep=CovParm Estimate); 
    set estrand;
    if CovParm="UN(5,2)";
run;

proc transpose data=results2 out=resbis (drop= _NAME_ ); 
   id CovParm;
   var Estimate;
run;

data res2;
   merge res2 resbis;
   nde = x_Sy;
   nie = Sm_x*Sy_m + UN_5_2_;
   tce = nde + nie;
   pctmed = nie /tce;
run;

proc append base=res data=res2; run;

%end;
%mend bootstrap;

/*Run the bootstrap macro*/
%bootstrap;

data res; set res;
if _N_ = 1 then delete; run;

title '95% bootstrap confidence interval for the NDE';
proc univariate data=res ;
output out=nde pctlpts=2.5 97.5 pctlpre=P;
var nde;
run;

title '95% bootstrap confidence interval for the NIE';
proc univariate data=res ;
output out=nie pctlpts=2.5 97.5 pctlpre=P;
var nie;
run;

title '95% bootstrap confidence interval for the TCE';
proc univariate data=res ;
output out=tce pctlpts=2.5 97.5 pctlpre=P;
var tce;
run;

title '95% bootstrap confidence interval for the percent mediated effect';
proc univariate data=res ;
output out=pctmed pctlpts=2.5 97.5 pctlpre=P;
var pctmed;
run;

title '95% bootstrap confidence interval for the NDE';
proc print data=nde; run;
title '95% bootstrap confidence interval for the NIE';
proc print data=nie; run;
title '95% bootstrap confidence interval for the TCE';
proc print data=tce; run;
title '95% bootstrap confidence interval for the percent mediated effect';
proc print data=pctmed; run;








