libname RS "C:\NCME 2016\RMPW\";
options nofmterr; 
*---------------------------------------------------------------;
* Estimate Propensity Scores; 
* --------------------------------------------------------------;
%let cov=emp_prior pqtrunc50 pqtrunc51 pqtrunc52 pqtrunc53 pqtrunc30 hispanic pqtrunc49	nevmar; 
                                                                          /*specify covariates*/
%macro PS (indata, tr, mv, tmv, outdata);
  proc sql noprint;                  
    select distinct min(&tr), max(&tr), min(&mv),max(&mv)
    into :ltr separated by '',:htr separated by '',:lmv separated by '',:hmv separated by ''
    from &indata;
  quit; 
  data &outdata;
    set &indata(keep=&tr);
	id=_N_;
  run;
  %do i= &ltr %to &htr; 
    %if &tmv=1 %then %let link=logit; 
	%if &tmv=2 %then %let link=clogit;
	%if &tmv=3 %then %let link=glogit;
    proc logistic data=&indata outest=est outmodel=pout;
      model &mv=&cov/link=&link;
      where &tr=&i; 
    run;
    %if &tmv=1 or &tmv=3 %then %do; 
	  data rs; 
	    set &indata; 
		&tr=&i; 
	  run;
	%end;    
    %if &tmv=2 %then %do; 
      proc sql; 
        create table rs as
        select A.*,B.intercept_0 as int&i,&i as &tr
        from &indata (drop=&tr) A, est B;
      quit;
    %end;
    proc logistic inmodel=pout;
      score data=rs out=pred&i ;
    run;  
    data &outdata;
      merge pred&i(drop=&tr) &outdata;
        %do j=&lmv %to &hmv; 
          %if &tmv=1 or &tmv=3 %then %do; 
            lp&i&j=log(p_&j/(1-p_&j));
		  %end;
          %if &tmv=2 %then %do; 
            lp&i&j=int&i-log(p_&lmv/(1-p_&lmv)); 
		  %end;
	      rename p_&j=p&i&j;
        %end; 
    run;
  %end;    
%mend PS;
		
*** binary mediator; 
%PS  (indata=RS.rside,         /*specify input data*/
      tr=treat,                /*specify treatment*/
      mv=emp,                  /*specify mediator*/
	  tmv=1,                   /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=bpred);          /*specify output data*/

*** ordinal mediator; 
%PS  (indata=RS.rside,         /*specify input data*/
      tr=treat,                /*specify treatment*/
      mv=empcat,               /*specify mediator*/
	  tmv=2,                   /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=opred);          /*specify output data*/

*** multi-category mediator; 
%PS  (indata=RS.rside,         /*specify input data*/
      tr=treat,                /*specify treatment*/
      mv=empcat,               /*specify mediator*/
	  tmv=3,                   /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=cpred);          /*specify output data*/

*** moderated mediator ; 
* apply the macro above after recoding the treatment x moderator combination
  and constructing a new &tr variable; 
data mdata ; 
  set RS.rside;
  mtreat=.;
  if teen_parent=0 and treat=0 then mtreat=0; /*group first by moderator and then by treatment*/
  if teen_parent=0 and treat=1 then mtreat=1;
  if teen_parent=1 and treat=0 then mtreat=2; 
  if teen_parent=1 and treat=1 then mtreat=3;
run;
%PS  (indata=mdata,                  /*specify input data*/
      tr=mtreat,                     /*specify moderated treatment*/
      mv=emp,                        /*specify mediator*/
	  tmv=1,                         /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=mpred);                /*specify output data*/

*---------------------------------------------------------------;
* Identify Common Support; 
* --------------------------------------------------------------;
%macro CS (indata, tr, mv, tmv, outdata);
  proc sql noprint;                  
    select distinct min(&tr), max(&tr), min(&mv),max(&mv)
    into :ltr,:htr,:lmv,:hmv from &indata;
  quit; 
  %do i=&ltr %to &htr; 
    %do j=&lmv %to &hmv;
      proc sql noprint;
        create table temp&i&j as
        select *, min(lp&i&j)as min&i&j, max(lp&i&j)as max&i&j
        from &indata
        group by &tr, &mv;
      quit;
      proc sql noprint;
        select max(min&i&j), min(max&i&j), std(lp&i&j)
        into :lb&i&j,:ub&i&j,:st&i&j
        from temp&i&j;
      quit;
      data &indata;
        set &indata;
        if lp&i&j<(&&lb&i&j)-0.2*(&&st&i&j)or lp&i&j >(&&ub&i&j)+0.2*(&&st&i&j)then NV=1;
      run;
    %end; 
  %end;
  data &outdata (drop=NV);
    set &indata;
    if NV=1 then delete;
  run;
%mend CS;

*** binary mediator; 
%CS (indata=bpred,                 /*specify input data*/
     tr=treat,                     /*specify treatment*/
     mv=emp,                       /*specify mediator*/
	 tmv=1,                        /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
     outdata=bvalid);              /*specify output data*/

*** ordinal mediator; 
%CS  (indata=opred,                 /*specify input data*/
      tr=treat,                     /*specify treatment*/
      mv=empcat,                    /*specify mediator*/
	  tmv=2,                        /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=ovalid);              /*specify output data*/

*** same macro applies to multi-category mediator; 

*** moderated mediator ; 
* apply the macro above to each subsample defined by the moderator; 
data mpred0 mpred1;
  set mpred; 
  if teen_parent=0 then output mpred0;
  else output mpred1;
run;
%CS  (indata=mpred0,                 /*specify input data*/
      tr=mtreat,                     /*specify treatment*/
      mv=emp,                        /*specify mediator*/
	  tmv=1,                         /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=mvalid0);              /*specify output data*/
%CS  (indata=mpred1,                 /*specify input data*/
      tr=mtreat,                     /*specify treatment*/
      mv=emp,                        /*specify mediator*/
	  tmv=1,                         /*specify type of mediator: 1=binary, 2=ordinal, 3=multi-category*/
      outdata=mvalid1);              /*specify output data*/
*combine subsamples; 
data mvalid; 
  set mvalid0 mvalid1;
run;

*---------------------------------------------------------------;
* Generate Parametric RMPW weight; 
* --------------------------------------------------------------;
%macro WT (indata, tr, mv, outdata);
  proc sql noprint;                  
    select distinct min(&tr), max(&tr), min(&mv),max(&mv)
    into :ltr separated by '',:htr separated by '',:lmv separated by '',:hmv separated by ''
    from &indata;
  quit; 
  data &outdata;
    set &indata;
    %do j=&lmv %to &hmv;	
      if &tr=&ltr then rmpw=1;
      if &tr=&htr and &mv=&j then rmpw=p&ltr.&j/p&htr.&j; 
    %end; 	 
  run;
%mend WT;
	  
*** binary mediator; 
%WT  (indata=bvalid,                 /*specify input data*/
      tr=treat,                      /*specify treatment*/
      mv=emp,                        /*specify mediator*/
	  outdata=bweight);              /*specify output data*/

*** ordinal mediator; 
%WT  (indata=ovalid,                 /*specify input data*/
      tr=treat,                      /*specify treatment*/
      mv=empcat,                     /*specify mediator*/
	  outdata=oweight);              /*specify output data*/

*** moderated mediator; 
* apply the macro above to each subsample defined by the moderator (in this case mvalid0 mvalid1); 

%WT  (indata=mvalid0,                 /*specify input data*/
      tr=mtreat,                       /*specify treatment*/
      mv=emp,                         /*specify mediator*/
	  outdata=mweight0);              /*specify output data*/
%WT  (indata=mvalid1,                 /*specify input data*/
      tr=mtreat,                       /*specify treatment*/
      mv=emp,                         /*specify mediator*/
	  outdata=mweight1);              /*specify output data*/
data mweight; 
  set mweight0 mweight1;
run;

*---------------------------------------------------------------;
* Generate Nonparametric RMPW weight; 
* --------------------------------------------------------------;
%macro ST (indata, tr, mv, nst, outdata);
  proc sql noprint;                  
    select distinct min(&tr), max(&tr), min(&mv),max(&mv)
    into :ltr separated by '', :htr separated by '', :lmv separated by '',:hmv separated by ''
    from &indata;
  quit; 
  proc rank data=&indata out=temp0 descending groups=&nst;
    var lp&ltr&hmv; 
    ranks rk&ltr;
  run;
  proc sort data=temp0;
    by rk&ltr; 
  run;  
  proc rank data=temp0 out=temp1 descending groups=&nst;
    by rk&ltr;
    var lp&htr&hmv; 
    ranks rk&htr;
  run;
  proc summary data=temp1;
    class rk&ltr rk&htr;
    ways 2;
    output out=temp2(drop=_type_ _freq_);
  run;  
  data temp2; 
    set temp2;
	id=_n_;
	rename id=strata;
  run;
  proc sql;
    create table stra as 
    select A.*, B.strata from 
    temp1 A, temp2 B
    where A.rk&ltr=B.rk&ltr and A.rk&htr=B.rk&htr;
  quit;
  proc sql;
    create table st1 as 
    select *, mean(&mv)as pr
	from stra 
    group by &tr, strata; 
  quit; 
  %do i=&ltr %to &htr; 
    %do j=&lmv %to &hmv;
      %do k=1 %to %eval(&nst*&nst); 
        data _null_;
          set st1;
		  if &tr=&i and &mv=&j and strata=&k then do; 
          call symput("pr&i&hmv&k",pr);
		  end; 
	    run; 
	  %end; 
    %end; 
  %end; 
  data &outdata; 
    set st1; 
    %do k=1 %to %eval(&nst*&nst);
      if strata=&k then do; 
        np&ltr=&&pr&ltr&hmv&k; 
	    np&htr=&&pr&htr&hmv&k; 
	  end;  
    %end;   
	if &tr=&ltr then nrmpw=1; 
	if &tr=&htr and &mv=&hmv then nrmpw=np0/np1;
	if &tr=&htr and &mv=&lmv then nrmpw=(1-np0)/(1-np1); 
  run;  
%mend ST; 

*** non-parametric rmpw for binary mediator 3x3; 
%ST (indata=bvalid,                    /*specify input data*/
      tr=treat,                         /*specify treatment*/
	  mv=emp,                           /*specify mediator*/
	  nst=3,                            /*specify number of strata for each dimension*/
	  outdata=bnpw3);                   /*specify output data*/

* non-parametric rmpw for binary mediator 4x4; 
%ST (indata=bvalid,                     /*specify input data*/
      tr=treat,                         /*specify treatment*/
	  mv=emp,                           /*specify mediator*/
	  nst=4,                            /*specify number of strata for each dimension*/
	  outdata=bnpw4);                   /*specify output data*/
 
*---------------------------------------------------------------;
* Generate Duplicate Observations; 
* --------------------------------------------------------------;
%macro DP (indata, tr, wt, outdata);
  data &outdata;
    set &indata (drop=D);
    do D = 0 to &tr;
    output;
    end; 
  run; 
  data &outdata;
    set &outdata;
	if D=1 then &wt=1;
  run;
%mend DP; 	

* binary mediator;   
  *parametic weight; 
%DP (indata=bweight,                    /*specify input data*/
     tr=treat,                          /*specify treatment*/
	 wt=rmpw,                          /*specify weight*/
	 outdata=bndata1);                  /*specify output data*/

*nonparametric weight 3x3; 
%DP (indata=bnpw3,                     /*specify input data*/
     tr=treat,                         /*specify treatment*/
	 wt=nrmpw,                         /*specify weight*/
	 outdata=bndata3);                 /*specify output data*/

*nonparametric weight 4x4; 
%DP (indata=bnpw4,                     /*specify input data*/
     tr=treat,                         /*specify treatment*/
	 wt=nrmpw,                         /*specify weight*/
	 outdata=bndata4);                 /*specify output data*/

* ordinal mediator;   
%DP (indata=oweight,                   /*specify input data*/
     tr=treat,                         /*specify treatment*/
	 wt=rmpw,                          /*specify weight*/
	 outdata=ondata);                  /*specify output data*/

* moderated mediator; 
%DP (indata=mweight,                  /*specify input data*/
     tr=treat,                        /*specify treatment*/
	 wt=rmpw,                         /*specify weight*/
	 outdata=mndata);                 /*specify output data*/

*---------------------------------------------------------------;
* Run Outcome models; 
* --------------------------------------------------------------;
proc surveyreg data=bndata1; 
  cluster id; 
  model trunc_dep12sm2 =treat D pqtrunc50 pqtrunc51	pqtrunc52 pqtrunc53 pqtrunc30 
        hispanic pqtrunc49 nevmar; 
  weight rmpw; 
run;

%let con=treat D pqtrunc50 pqtrunc51 pqtrunc52 pqtrunc53 pqtrunc30 hispanic pqtrunc49 nevmar; 
                                                           /*specify covariates*/
%macro OM (mydata,oc,tr,md,wt); 
  %if &md=N %then %do;
    proc surveyreg data = &mydata;
      cluster id;
      model &oc= &tr D &cov; 
      weight &wt;
    run;
  %end;
  %if &md^=N %then %do;
    proc surveyreg data = &mydata;
      cluster id;
      model &oc= &tr D &md &md*&tr &md*D &cov; 
      weight &wt;
    run;
  %end;
%mend OM; 

/** trunc_dep12sm12 as outcome **/

*** binary mediator, parametric rmpw ; 
%OM (mydata=bndata1,                   /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=rmpw);                         /*specify weight*/

*** binary mediator, nonparametric rmpw3x3; 
%OM (mydata=bndata3,                   /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=nrmpw);                        /*specify weight*/

*** binary mediator, nonparametric weight 4x4; 
%OM (mydata=bndata4,                   /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=nrmpw);                        /*specify weight*/

*** ordinal mediator, parametric weight; 
%OM (mydata=ondata,                    /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=rmpw);                         /*specify weight*/

*** moderated mediator, parametric weight; 
data mndata;
  set mndata;
  not_teen=1-teen_parent; 
run;

%OM (mydata=mndata,                    /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=teen_parent,                   /*specify moderator-if none, put N*/
	 wt=rmpw);                        /*specify weight*/

%OM (mydata=mndata,                    /*specify data*/
     oc=trunc_dep12sm2,                /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=not_teen,                      /*specify moderator-if none, put N*/
	 wt=rmpw);                        /*specify weight*/

/** depsmrk2 as outcome **/

*** binary mediator, parametric rmpw ; 
%OM (mydata=bndata1,                   /*specify data*/
     oc=depsmrk2,                      /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=rmpw);                         /*specify weight*/

*** binary mediator, nonparametric weight 4x4; 
%OM (mydata=bndata4,                   /*specify data*/
     oc=depsmrk2,                      /*specify outcome*/
     tr=treat,                         /*specify treatment*/
	 md=N,                             /*specify moderator-if none, put N*/
	 wt=nrmpw);                        /*specify weight*/

/** SEM/Path Analysis Approach **/

proc logistic data=RS.rside;
  model emp=treat;
  ods output ParameterEstimates=out1;
run;
data _null_;
  set out1 (where =(variable='treat')) ;
  call symput("mediator",estimate);
  call symput ("mediator_se", stderr);
run;
proc reg data=RS.rside; 
  model trunc_dep12sm2=treat emp; 
  ods output ParameterEstimates=out2;
run; 
data _null_;
  set out2;
  if variable='treat' then do;
  call symput("treat_b", estimate);
  call symput ("treat_se", stderr);
  end;
  if variable ='emp' then do; 
  call symput("emp_b", estimate);
  call symput ("emp_se", stderr);
  end;
run;
%put &mediator &mediator_se &treat_b &treat_se &emp_b &emp_se; 
data sem; 
  direct=&treat_b;
  direct_se=&treat_se;
  indirect=&emp_b*&treat_b;
  indirect_se=sqrt((&mediator*&mediator)*(&emp_se*&emp_se)+(&mediator_se*&mediator_se)*(&emp_b*&emp_b)); 
run;
proc print data=sem;
run;