/***************************************************************************************************************************************************************/
/***************************************************************************************************************************************************************/

/* 						Marginal structural model approach for mediation for a continuous mediator and outcome, 												*/
/*											and a multicategorical exposure (treatment) with four levels 																	*/



/****************************************************************************************************************************************************************/
/* 																	DATA MANIPULATION 																			*/
/*																																								*/
/* The four exposure levels a0, a1, a2, a3 (with a0 as reference level) have to be coded as 0,1,2,3, respectively.												*/
/* Categorical adjustment variables need to be coded as a series of dummy variables before being entered as covariates in the exposure and mediator models.		*/										
/*																																								*/							
/*																	MACRO VARIABLES																				*/
/*																																								*/
/* The main macro is %MSM_mediation; to use it, the user needs to specify the values for the following macro variables:											*/
/*																																								*/							
/* input_data:  	input data that include the multicategorical exposure variable, the continuous outcome and mediator variables, as well as the covariates	*/
/*					to be adjusted for in the models; input data also has to include the CLUSTER VARIABLE if the CLUSTER BOOTSTRAP is planned to be performed; 	*/
/* A: 				the name of the multicategorical exposure with four levels a0, a1, a2, a3 (coded as 0,1,2,3, respectively); it is not necessary to create 	*/
/*					the dummy variables for the exposure levels a1, a2, a3;																						*/											
/* M:				the name of the continuous mediator variable;																								*/
/* Y:				the name of the continuous outcome variable;																								*/
/* interaction: 	the user needs to specify INTERACTION=0 or INTERACTION=1 for the outcome model without or with interaction between the exposure and the 	*/
/*					mediator, respectively;																														*/	
/* adjusted: 		the user needs to specify ADJUSTED=0 or ADJUSTED=1 to obtain unadjusted or adjusted NIE, NDE and TE, respectively;							*/
/* cvar_A:			adjustment variables (covariates) for the exposure model; categorical variables need to be coded in input data as a series of dummy 		*/
/*					variables before being entered as covariates; 																								*/															
/* cvar_M:			adjustment variables (covariates) for the mediator model; categorical variables need to be coded in input data as a series of dummy 		*/
/*					variables before being entered as covariates;																								*/
/* ridging: 		the user has to specify the value of this macro variable ONLY for adjusted analyses (ADJUSTED=1) to define the ridging method to be used 	*/
/*              	in the adjusted multinomial logistic regression; the possible values are ABSOLUTE, RELATIVE, NONE (Remark: RIDGING=RELATIVE is used by 		*/
/*              	default in proc logistic); in other words, the user has to specify RIDGING=RELATIVE if they would like to use the SAS default option;		*/
/* truncation:  	the user has to specify the value of this macro variable (TRUNCATION=1) ONLY if right (largest) weight truncation is planned;               */ 
/* cut_off: 		if weight truncation is planned (TRUNCATION=1 was specified), the user has to specify the cut-off value (between 0 and 100); 				*/
/*	            	for example, the user has to set CUT_OFF=95 for a 5% right (largest) weight truncation; 													*/								
/* boot:	    	to calculate confidence intervals using robust standard errors, the user needs to specify BOOT=0; to calculate confidence intervals by      */
/*					bootstrapping, the user needs to specify BOOT=1;																							*/																
/* nboot:			the user needs to specify a value for the number of bootstrap samples to consider if BOOT=1 was specified; this value is an upper bound and */
/*					can exceed the actual number of bootstrap samples used (see variable n below); 																*/
/* n:				the user needs to specify the desired number (<= nboot) of bootstrap samples without convergence failure in the exposure model; 			*/
/*					these n non problematic bootstrap samples will be randomly selected from the nboot bootstrap samples and will be used to construct 95% CI   */
/*					(e.g., in our adjusted analyses, the estimated percentage of convergence failure was close to 5% and we took nboot=6000 and n=5000); 		*/
/*					for crude analyses, the user typically specifies the same value as for nboot;																*/
/* SAMPLINGUNIT: 	the user needs to specify the name of the cluster variable (for cluster sampling) if BOOT=1 was specified for CLUSTER BOOTSTRAP.			*/										
/*																																								*/
/* The user needs to specify the macro variable values ONLY for the main macro %%MSM_mediation. 																*/
/*																																								*/														
/*																		SAS OUTPUT EXPLANATION 																	*/																			 
/* After %MSM_mediation macro execution, the following information is displayed in the SAS results viewer window:												*/		
/* type of analysis (crude or adjusted, with or without interaction);																							*/
/* method used to calculate the confidence intervals for NIE, NDE and TE (CI based on robust standard errors, bootsrap CI);										*/																			
/* mediation effect point estimates (e.g., contrast_1 displays the mediation effects if the exposure changed from a0 to a1);									*/																
/* mediation effect interval estimates (e.g., CI_1_low and CI_1_upp are the upper and lower bounds of the 95% confidence interval for the contrast_1).  		*/													
/* 																																								*/
/****************************************************************************************************************************************************************/


%macro point_estimates(data);

%global conv_status;

data mydata;                                                             
set &data;  
IDtemp=_n_; 
run; 

%if &adjusted=0 %then %do;

	%let conv_status=0;

	data mydata; 
	set mydata; 
	Atemp = &A; 
	A_duplic = 0; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
	A_duplic = 1; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
	A_duplic = 2; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
	A_duplic = 3; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
	run;

	ods output  ModelFit = Model_Fit (keep=ValueDF);
	proc genmod data = mydata; 
	class Atemp (param=ref ref='0');
	model Mtemp = Atemp; 
	output out = mydata P = predM;
	run;

	data rmse; set Model_Fit(obs=1); _RMSE_=SQRT(ValueDF); keep _RMSE_; run;

	data mydata;
	set mydata;
	if _n_=1 then set rmse;
	Atemp = A_duplic;
	weightDIR=PDF('normal', &M, predM, _RMSE_);
	run;

	proc genmod data=mydata; 
	class Atemp (param=ref ref='0');
	model Mtemp = Atemp;
	output out = mydata P = predM_duplic;
	run;

	data mydata;
	set mydata;
	weightINDIR = PDF('normal', &M, predM_duplic, _RMSE_);
	w = weightINDIR/weightDIR;
	run;

%end;

%if &adjusted=1 %then %do;

	proc logistic data = mydata noprint;
	model &A = / link = glogit;
	output out = mydata predprobs = i;
	run;

	data mydata;
	set mydata;
	if &A = 0 then prA=IP_0;
	if &A = 1 then prA=IP_1;
	if &A = 2 then prA=IP_2;
	if &A = 3 then prA=IP_3;
	drop _FROM_ _INTO_ IP_0-IP_3;
	run;

	ods output  ConvergenceStatus=convergence (keep=Status);
	proc logistic data = mydata;
	model &A = &cvar_A /link=glogit ridging=&ridging;
	output out=mydata predprobs=i;
	run; 

	proc sql noprint;
	select Status
 	into :conv_status trimmed
	from convergence;
	quit;

	/* %put status=&conv_status ; */

	%if &conv_status=0 %then %do;

		data mydata;
		set mydata;
		if &A = 0 then do; prAcondC=IP_0; end;
		if &A = 1 then do; prAcondC=IP_1; end;
		if &A = 2 then do; prAcondC=IP_2; end;
		if &A = 3 then do; prAcondC=IP_3; end;
		drop _FROM_ _INTO_ IP_0-IP_3;
		run;

		data mydata; 
		set mydata;
		Atemp= &A;
		A_duplic = 0; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
		A_duplic = 1; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
		A_duplic = 2; Mtemp=.; if &A = A_duplic then Mtemp = &M; output;
		A_duplic = 3; Mtemp=.; if &A = A_duplic then Mtemp = &M ; output;
		run;

		

		ods output  ModelFit = Model_Fit (keep=ValueDF);
		proc genmod data = mydata; 
		class Atemp (param=ref ref='0');
		model Mtemp = Atemp &cvar_M; 
		output out=mydata P=predM;
		run;

		data rmse; set Model_Fit(obs=1); _RMSE_=SQRT(ValueDF); keep _RMSE_; run;
 
		data mydata;
		set mydata;
		if _n_=1 then set rmse;
		weightDIR=PDF('normal', &M, predM, _RMSE_);
		Atemp = A_duplic;
		run;

		proc genmod data = mydata; 
		class Atemp (param=ref ref='0');
		model Mtemp= Atemp &cvar_M; 
		output out=mydata P=predM_duplic;
		run;

		data mydata;
		set mydata;
		weightINDIR=PDF('normal', &M, predM_duplic, _RMSE_);
		w = (prA/prAcondC)*(weightINDIR/weightDIR);
		run;

	%end;

%end;



%if &conv_status=0 %then %do;

	%if &truncation = 1 %then %do;

		proc univariate data=mydata noprint;
		var w;
		output out=Percentile_w pctlpts = &cut_off pctlpre =P_;
		run;

		data mydata;
		set mydata;
		if _n_=1 then set Percentile_w;
		if w > P_&cut_off then w = P_&cut_off; 
		run;
				
	%end;

	%if &interaction = 0 %then %do;
	
		ods output  GEEEmpPEst=Model_ParamEst (keep=Estimate);
		proc genmod data = mydata; 
		class IDtemp &A (param=ref ref='0') A_duplic (param=ref ref='0');
		model &Y = &A A_duplic;
		weight w;
		repeated subject=IDtemp / type=ind;
		run;

		data param; set Model_ParamEst (firstobs=2); run;
		proc transpose data=param out=param; run;

		data point_estimates; 
		set param;
		rename COL1=NDE_a1 COL2=NDE_a2 COL3=NDE_a3 COL4=NIE_a1 COL5=NIE_a2 COL6=NIE_a3;
		TE_a1=COL1+COL4; TE_a2=COL2+COL5; TE_a3=COL3+COL6;
		keep COL1-COL6 TE_a1-TE_a3;
		run;


	%end;

	%if &interaction = 1 %then %do;
				
		ods output  GEEEmpPEst = Model_ParamEst (keep=Estimate);
		proc genmod data = mydata; 
		class IDtemp &A (param=ref ref='0') A_duplic (param=ref ref='0');
		model &Y = &A | A_duplic;
		weight w;
		repeated subject=IDtemp / type=ind;
		run;

		data param; set Model_ParamEst (firstobs=2); run;
		proc transpose data=param out=param; run;
	
		data point_estimates ; 
		set  param;
		rename COL1=NDE_a1 COL2=NDE_a2 COL3=NDE_a3;
		NIE_a1= COL4+COL7; NIE_a2= COL5+COL11; NIE_a3= COL6+COL15;
		TE_a1=COL1+COL4+COL7; TE_a2=COL2+COL5+COL11; TE_a3=COL3+COL6+COL15;
		keep COL1-COL3 NIE_a1-NIE_a3 TE_a1-TE_a3;
		run;

	%end;

	
%end;

%mend point_estimates;





%macro CI(effect);

proc univariate data= boot noprint;
var &effect._a1 &effect._a2 &effect._a3;
output out=CI_&effect pctlpts=2.5 97.5 pctlpre=&effect._1 &effect._2 &effect._3 pctlname= _low _upp ;
run;

data CI_&effect; 
retain &effect._1_low &effect._1_upp &effect._2_low &effect._2_upp &effect._3_low &effect._3_upp; set CI_&effect;
run;

%mend;




%macro MSM_mediation(input_data,A,M,Y,interaction,adjusted,cvar_A,cvar_M,ridging,truncation,cut_off,boot,nboot,SAMPLINGUNIT,n);

%if &adjusted=0 %then %do; 	
	%let pref=crude;
%end;
%if &adjusted=1 %then %do; 	
	%let pref=adjusted;
%end;

%point_estimates(data=&input_data);
data contrasts; set point_estimates; run;


/* 95% CI using robust standard errors */
%if &boot=0 %then %do; 

ods html close;
options nonotes nosource nosource2 errors=0; ods results off;

	%if &interaction=0 %then %do;	

		ods output  GEEEmpPEst=Model_ParamEst (keep=LowerCL UpperCL)  GEERCov=cov_matrix;
		proc genmod data = mydata; 
		class IDtemp &A (param=ref ref='0') A_duplic (param=ref ref='0');
		model &Y = &A A_duplic;
		weight w;
		repeated subject=IDtemp / type=ind covb;
		run;

		proc iml;

		use Model_ParamEst; read point (2:7) into CI_NDE_NIE;
	
		use contrasts; read all into contrasts; 
	
		use cov_matrix; read all var _NUM_ into cov[colname=varNames]; 
		cov=cov[2:nrow(cov),2:nrow(cov)]; 

		se_TE_a1=sqrt(cov[1,1]+cov[4,4]+2*cov[1,4]); 
		se_TE_a2=sqrt(cov[2,2]+cov[5,5]+2*cov[2,5]); 
		se_TE_a3=sqrt(cov[3,3]+cov[6,6]+2*cov[3,6]); 

	
		TE_a1_low=contrasts[1,7]-quantile("Normal",.975)*se_TE_a1;
		TE_a1_upp=contrasts[1,7]+quantile("Normal",.975)*se_TE_a1;

		TE_a2_low=contrasts[1,8]-quantile("Normal",.975)*se_TE_a2;
		TE_a2_upp=contrasts[1,8]+quantile("Normal",.975)*se_TE_a2;

		TE_a3_low=contrasts[1,9]-quantile("Normal",.975)*se_TE_a3;
		TE_a3_upp=contrasts[1,9]+quantile("Normal",.975)*se_TE_a3;
	
	
		NDE_GEE_CI = contrasts[1,1]||CI_NDE_NIE[1,]      ||contrasts[1,2]||CI_NDE_NIE[2,]      ||contrasts[1,3]||CI_NDE_NIE[3,];
		NIE_GEE_CI = contrasts[1,4]||CI_NDE_NIE[4,]      ||contrasts[1,5]||CI_NDE_NIE[5,]      ||contrasts[1,6]||CI_NDE_NIE[6,]; 
		TE_GEE_CI  = contrasts[1,7]||TE_a1_low||TE_a1_upp||contrasts[1,8]||TE_a2_low||TE_a2_upp||contrasts[1,9]||TE_a3_low||TE_a3_upp;  

		all = NDE_GEE_CI// NIE_GEE_CI//TE_GEE_CI;

		create  results from all 
		[colname={"contrast_1" "CI1_low" "CI1_upp" "contrast_2" "CI2_low" "CI2_upp" "contrast_3" "CI3_low" "CI3_upp"}];
		append from all; close results;

		quit;

	%end;

	%if &interaction=1 %then %do;	
	
		ods output  GEEEmpPEst=Model_ParamEst (keep=LowerCL UpperCL)  GEERCov=cov_matrix;
		proc genmod data = mydata; 
		class IDtemp &A (param=ref ref='0') A_duplic (param=ref ref='0');
		model &Y = &A|A_duplic;
		weight w;
		repeated subject=IDtemp / type=ind covb;
		run;

		proc iml;

		use Model_ParamEst; read point (2:4) into CI_NDE;
	
		use contrasts; read all into contrasts; 

		use cov_matrix; read all var _NUM_ into cov[colname=varNames];
		cov=cov[2:nrow(cov),2:nrow(cov)];

		
		se_NIE_a1=sqrt(cov[4,4]+cov[7,7]  +2*cov[4,7]);  
		se_NIE_a2=sqrt(cov[5,5]+cov[11,11]+2*cov[5,11]); 
		se_NIE_a3=sqrt(cov[6,6]+cov[15,15]+2*cov[6,15]);



		se_TE_a1=sqrt(cov[1,1]+cov[4,4]+cov[7,7]  +2*cov[1,4]+2*cov[4,7] +2*cov[1,7]);  
		se_TE_a2=sqrt(cov[2,2]+cov[5,5]+cov[11,11]+2*cov[2,5]+2*cov[5,11]+2*cov[2,11]); 
		se_TE_a3=sqrt(cov[3,3]+cov[6,6]+cov[15,15]+2*cov[3,6]+2*cov[6,15]+2*cov[3,15]); 


		
		NIE_a1_low=contrasts[1,4]-quantile("Normal",.975)*se_NIE_a1;
		NIE_a1_upp=contrasts[1,4]+quantile("Normal",.975)*se_NIE_a1; 

		NIE_a2_low=contrasts[1,5]-quantile("Normal",.975)*se_NIE_a2; 
		NIE_a2_upp=contrasts[1,5]+quantile("Normal",.975)*se_NIE_a2;

		NIE_a3_low=contrasts[1,6]-quantile("Normal",.975)*se_NIE_a3;
		NIE_a3_upp=contrasts[1,6]+quantile("Normal",.975)*se_NIE_a3;
	


		TE_a1_low=contrasts[1,7]-quantile("Normal",.975)*se_TE_a1; 
		TE_a1_upp=contrasts[1,7]+quantile("Normal", .975)*se_TE_a1;

		TE_a2_low=contrasts[1,8]-quantile("Normal",.975)*se_TE_a2; 
		TE_a2_upp=contrasts[1,8]+quantile("Normal",.975)*se_TE_a2;

		TE_a3_low=contrasts[1,9]-quantile("Normal",.975)*se_TE_a3; 
		TE_a3_upp=contrasts[1,9]+quantile("Normal",.975)*se_TE_a3; 
	
	
		NDE_GEE_CI = contrasts[1,1]||CI_NDE[1,]            ||contrasts[1,2]||CI_NDE[2,]            ||contrasts[1,3]||CI_NDE[3,]; 
		NIE_GEE_CI = contrasts[1,4]||NIE_a1_low||NIE_a1_upp||contrasts[1,5]||NIE_a2_low||NIE_a2_upp||contrasts[1,6]||NIE_a3_low||NIE_a3_upp; 
		TE_GEE_CI  = contrasts[1,7]||TE_a1_low ||TE_a1_upp ||contrasts[1,8]||TE_a2_low ||TE_a2_upp ||contrasts[1,9]||TE_a3_low ||TE_a3_upp;  

		all = NDE_GEE_CI// NIE_GEE_CI//TE_GEE_CI;

		create  results from all 
		[colname={"contrast_1" "CI_1_low" "CI_1_upp" "contrast_2" "CI_2_low" "CI_2_upp" "contrast_3" "CI_3_low" "CI_3_upp"}];
		append from all;close results;

		quit;

		%end;
	
	options notes source source2 errors=20; ods results on;
	ods html;

	data &pref._results_GEE; retain effect; set results;
	if _n_=1 then effect="NDE";	if _n_=2 then effect="NIE";	if _n_=3 then effect="TE";
	array _nums {*} _numeric_;
	do i = 1 to dim(_nums);
 		_nums{i} = round(_nums{i},.0001);
	end;
	drop i;
	run;
	title "&pref analysis with 95% CI based on robust standard errors, interactions = &interaction";
	proc print data= &pref._results_GEE; run;
	title;	

%end; /* end boot=0 */


/* 95% CI by bootstrapping */

%if &boot=1 %then %do;

	proc datasets library = WORK noprint; delete  boot_0 convergency_count; run; quit;

	proc surveyselect data = &input_data
	out=bootdata seed=1983 method=urs samprate=1 outhits rep=&nboot noprint;
	SAMPLINGUNIT &SAMPLINGUNIT; 
	run;

	ods html close;
	options nonotes nosource nosource2 errors=0; ods results off;

	%do i=1 %to &nboot;
	%put iter=&i;

		data data_for_boot; set bootdata (where=(Replicate=&i));  run; 

		/*
		proc sql noprint; select count(*) into :OBSCOUNT from data_for_boot; quit;
		%put Count=&OBSCOUNT.;
		*/

		%point_estimates(data=data_for_boot);

		
		data convergency; conv_status=&conv_status; run; 
		proc append base = convergency_count data = convergency force; run;

		%if &conv_status = 0 %then %do;
			proc append base = boot_0 data = point_estimates force; run;
		%end;

	%end;

	proc sql noprint; select count(*) into :OBSCOUNT from boot_0; quit;
	%put Count_boot_0=&OBSCOUNT.;


	/* to select bootstrap iterations without convergence problems */
	proc surveyselect data = boot_0
	out=boot seed=123 method=srs n=&n noprint;
	run;


	%CI(effect=NDE);
	%CI(effect=NIE);
	%CI(effect=TE);


	proc iml;
	use contrasts; read all into contrasts; 

	use CI_NDE; read all into NDE_CI;
	use CI_NIE; read all into NIE_CI; 
	use CI_TE;  read all into TE_CI; 

	NDE_with_boot_CI = contrasts[1,1]|| NDE_CI[1,1:2]||contrasts[1,2]||NDE_CI[1,3:4]||contrasts[1,3]||NDE_CI[1,5:6]; 
	NIE_with_boot_CI = contrasts[1,4]|| NIE_CI[1,1:2]||contrasts[1,5]||NIE_CI[1,3:4]||contrasts[1,6]||NIE_CI[1,5:6]; 
	TE_with_boot_CI  = contrasts[1,7]|| TE_CI[1,1:2] ||contrasts[1,8]||TE_CI[1,3:4] ||contrasts[1,9]||TE_CI[1,5:6];  

	all = NDE_with_boot_CI// NIE_with_boot_CI//TE_with_boot_CI;

	create  results from all 
	[colname={"contrast_1" "CI_1_low" "CI_1_upp" "contrast_2" "CI_2_low" "CI_2_upp" "contrast_3" "CI_3_low" "CI_3_upp"}];
	append from all;close results;

	quit;

	options notes source source2 errors=20; ods results on;
	ods html;


	data &pref._results_boot; retain effect; set results;
	if _n_=1 then effect="NDE";	if _n_=2 then effect="NIE";	if _n_=3 then effect="TE";
	array _nums {*} _numeric_;
	do i = 1 to dim(_nums);
 	 _nums{i} = round(_nums{i},.0001);
	end;
	drop i;
	run;
	title "&pref analysis with bootsrap 95% CI, interactions = &interaction";
	proc print data= &pref._results_boot; run;
	title;	

%end; /* end boot =1 */

%mend;