/******************************************************************************************************************************************/
/*********************************************            MACRO NPC          **************************************************************/
/******************************************************************************************************************************************/
/*

LIST OF MACRO PARAMETERS:

dati	=	name of the dataset;
var_byn	=	list of binary variables;
var_cat	=	list of categorical variables, non-binary;
var_con	=	list of continuous variables;
dom_byn	=	list of directional marginal sub-hypothesis for binary variables:
			if XAu< XBu specify LESS
			if XAu> XBu specify GREAT
			if XAu=/ XBu specify NOTEQ;
dom_con	=	list of directional marginal sub-hypothesis for continuous variables, see above;
weights	=	list of weights for the variables, specify first weights for binary variables, 
			then weights for categorical variables and for continuous variables;
clas	=	variable defining the two groups (character variable);
nsample	=	number of conditional resamplings;
strato	=	variable defining strata (character variable);
paired	=	paired data, (character variable)(specify yes/no);
unit	=	variable identifying paired observations;
missing	=	presence of missing values, (specify yes/no).

*/



%macro npc(dati=, var_byn=, var_cat=, var_con=, dom_byn=, dom_con=, weights=, clas=, nsample=, strato=, 
      paired=, unit=, missing=);

	%macro conta(lista,num);
		%global &num;
		%let &num=0;
		%if &lista ne %then %do;
			data _null_;
			array lis [*] &lista;
			n=dim(lis);
			call symput("&num",trim(left(put(n,8.))));
			run;
		%end;
	%mend conta;
	%conta(&var_byn,nvar_byn);
	%conta(&var_cat,nvar_cat);
	%conta(&var_con,nvar_con);
	%let lis_var=&var_byn &var_cat &var_con;	
	%let var_con=&var_byn &var_con;	
	%let dom_con=&dom_byn &dom_con;
	%let nvar=%eval(&nvar_byn+&nvar_cat+&nvar_con);
	%let nvar_con=%eval(&nvar_byn+&nvar_con);
	%if %upcase(&paired)=YES and &strato ne %then %do;
		proc sort data=&dati;
			by &strato;
		run;
	%end;
	proc multtest data=&dati perm nsample=&nsample outsamp=permut(drop=_obs_ 
			%if %upcase(&paired)=YES %then %do;
			rename=(_strata_=&unit %if &strato ne %then &strato=_strata_;) %end;) noprint;
		%if %upcase(&paired)=YES %then %do;
			%if &strato ne %then by &strato;;
			strata &unit; 
		%end;
		%else %if &strato ne %then strata &strato;;
		class &clas;
		test mean(&lis_var);
	run;
	data dati(keep=_sample_ %if &strato ne %then _strata_; _class_ &lis_var 
			%if %upcase(&paired)=YES %then &unit;);
		set &dati;
		_SAMPLE_=0;
		rename %if &strato ne %then &strato=_STRATA_; &clas=_CLASS_;
	run;
	data permut;
		set dati permut;
	run;
	%if %upcase(&paired)=YES %then %do;
		proc sort data=permut;
			by _sample_ %if &strato ne %then _strata_; &unit;
		run;
		data permut(drop=&unit _class_ %do i=1 %to &nvar; diff&i %end;);
			set permut(rename=(%do i=1 %to &nvar; %scan(&lis_var,&i)=diff&i %end;));
			by _sample_ %if &strato ne %then _strata_; &unit;
			if first.&unit then do;
				retain &lis_var 0;
				%do i=1 %to &nvar; %scan(&lis_var,&i)=diff&i; %end;
			end;
			else do; 
				%do i=1 %to &nvar; %scan(&lis_var,&i)=%scan(&lis_var,&i)-diff&i; %end;
			end;
			if last.&unit;
		run;
		%if &var_con ne %then %do;
			data statis(drop=%do i=1 %to &nvar_con; sum&i test&i %end;);
				set permut(keep=_sample_ %if &strato ne %then _strata_; &var_con 
					rename=(%do i=1 %to &nvar_con; %scan(&var_con,&i)=test&i %end;));
				by _sample_ %if &strato ne %then _strata_;;
				if first.%if &strato ne %then _strata_; %else _sample_; then do;
						%do i=1 %to &nvar_con; sum&i=0; %end;
				end;
				%do i=1 %to &nvar_con; 
					sum&i+test&i;
					%if %upcase(%scan(&dom_con,&i))=LESS %then %do; 
						%scan(&var_con,&i)= (-1)*sum&i; %end;
					%else %if %upcase(%scan(&dom_con,&i))=GREAT %then %do; 
						%scan(&var_con,&i)=sum&i; %end;								
					%else %if %upcase(%scan(&dom_con,&i))=NOTEQ %then %do;
						%scan(&var_con,&i)=abs(sum&i); %end;
				%end;
				if last.%if &strato ne %then _strata_; %else _sample_;;
			run;
		%end;
	%end;
	%else %if %upcase(&paired)=NO and &var_con ne %then %do;
		proc sort data=permut;
			by _sample_ %if &strato ne %then _strata_; _class_;
		run;
		data somme(drop=_class_ %do i=1 %to &nvar_con; test&i %end;);
			set permut(keep=_sample_ %if &strato ne %then _strata_; _class_ &var_con 
				rename=(%do i=1 %to &nvar_con; %scan(&var_con,&i)=test&i %end;));
			by _sample_ %if &strato ne %then _strata_; _class_;
			if first._class_ then do;
				%do i=1 %to &nvar_con; sum&i=0; n_mis&i=0; %end;
			end;
			%do i=1 %to &nvar_con; 
				sum&i+test&i; if test&i ne . then n_mis&i+1;
			%end;
			if last._class_;
		run;
		proc sort data=somme;
			by _sample_ %if &strato ne %then _strata_;;
		run;
		data statis(keep=_sample_ %if &strato ne %then _strata_; &var_con);
			set somme;
			by _sample_ %if &strato ne %then _strata_;;
			if first.%if &strato ne %then _strata_;%else _sample_; then do;
				retain %do i=1 %to &nvar_con; varr&i m_rap&i %end; 0;
				%do i=1 %to &nvar_con; varr&i=sum&i; m_rap&i=n_mis&i; %end;
			end;
			else do;
				%do i=1 %to &nvar_con; 
					m_rap&i=m_rap&i/n_mis&i;
					m2_rap&i=1/m_rap&i;
					varr&i=(varr&i*sqrt(m2_rap&i))-(sum&i*sqrt(m_rap&i));
				%end;
			end;
			%do i=1 %to &nvar_con;
				%if %upcase(%scan(&dom_con,&i))=LESS %then %do; 
					%scan(&var_con,&i)=(-1)*varr&i; %end;
				%else %if %upcase(%scan(&dom_con,&i))=GREAT %then %do; 
					%scan(&var_con,&i)=varr&i; %end;
				%else %if %upcase(%scan(&dom_con,&i))=NOTEQ %then %do; 
					%scan(&var_con,&i)=abs(varr&i); %end;
			%end;
			if last.%if &strato ne %then _strata_; %else _sample_;; 
		run;
	%end;
	%if &var_cat ne %then %do;		
		proc sort data=permut;
			by _sample_ %if &strato ne %then _strata_;;
		run;		
		proc freq data=permut noprint;
			by _sample_ %if &strato ne %then _strata_;;
			tables _class_ /missing out=c;
		run;
		proc sort data=c;
			by _sample_ %if &strato ne %then _strata_; _class_;
		run;
		proc freq data=permut noprint;
			%if &strato ne %then by _sample_;;
			tables %if &strato ne %then _strata_; %else _sample_; /missing out=d;
		run;
		proc sort data=d;
			by _sample_ %if &strato ne %then _strata_;;
		run;
		%do i=1 %to &nvar_cat;
			proc freq data=permut noprint;
				by _sample_ %if &strato ne %then _strata_;;
				tables %scan(&var_cat,&i)*_class_ /missing out=a&i outpct;
				tables %scan(&var_cat,&i)         /missing out=b&i;
			run;
			proc sort data=a&i;
				by _sample_ %if &strato ne %then _strata_; _class_;
			run;
			data ac&i;
				merge a&i(drop=percent count) c(drop=percent);
				by _sample_ %if &strato ne %then _strata_; _class_;
			run;
			proc sort data=ac&i;
				by _sample_ %if &strato ne %then _strata_; %scan(&var_cat,&i);
			run;
			proc sort data=b&i;
				by _sample_ %if &strato ne %then _strata_; %scan(&var_cat,&i);
			run;
			data abc&i;
				merge ac&i b&i(drop=count pct_col pct_row);
				by _sample_ %if &strato ne %then _strata_; %scan(&var_cat,&i);
			run;
			proc sort data=abc&i;
				by _sample_ %if &strato ne %then _strata_;;
			run;
			data abcd&i;
				merge abc&i d(drop=percent rename=(count=tot));
				by _sample_ %if &strato ne %then _strata_;;
			run;
			proc sort data=abcd&i;
				by _sample_ %if &strato ne %then _strata_; %scan(&var_cat,&i);
			run;
			data stat&i (keep=_sample_ %if &strato ne %then _strata_; %scan(&var_cat,&i));
				set abcd&i(rename=(%scan(&var_cat,&i)=var&i));
				pct_r=round(pct_row,1);
				by _sample_ %if &strato ne %then _strata_; var&i;
				retain p_prop 0;
				if first.var&i and pct_r=100 then p_prop=
					(tot-count)*percent+((count*(pct_col-percent)**2)/percent);
				else if first.var&i and pct_r ne 100 then p_prop=(count*(pct_col-
					percent)**2)/percent;
				else p_prop=p_prop+((count*(pct_col-percent)**2)/percent);	
			if first.%if &strato ne %then _strata_;%else _sample_; 
				then %scan(&var_cat,&i)=0;
				if last.var&i then %scan(&var_cat,&i)+p_prop;
				if last.%if &strato ne %then _strata_; %else _sample_;; 
			run;		
			proc sort data=stat&i;
				by _sample_ %if &strato ne %then _strata_;;
			run;
		%end;
		data statis;
			merge %do i=1 %to &nvar_cat; stat&i %end; %if &var_con ne %then statis;;
			by _sample_ %if &strato ne %then _strata_;;
		run;
	%end;
	%if &strato ne %then %do;
		proc sort data=statis;
			by _strata_;
		run;
	%end;
	proc rank data=statis descending ties=high out=pval_par;
		%if &strato ne %then by _strata_;;
 		var &lis_var;
  		ranks %do i=1 %to &nvar; test&i._r %end;;
	run;
	
	data pval_par(drop=%do i=1 %to &nvar; test&i._r %end;);
		set pval_par;
		%do i=1 %to &nvar; lambda&i=(0.5+test&i._r-1)/(&nsample+1); 
		loglam&i=log(lambda&i); %end;
		%if &weights ne %then %do; comb1=-2*(%do i=1 %to %eval(&nvar-1); 
			(%scan(&weights,&i)*loglam&i)+ %end; (%scan(&weights,&nvar)*loglam&nvar)); %end;
		%else %do; comb1=-2*(%do i=1 %to %eval(&nvar-1); loglam&i+ %end; loglam&nvar); %end;
	run;
	proc print data=pval_par noobs label;
		where _sample_=0;
		title "Partial Tests";
		var %if &strato ne %then _strata_; %do i=1 %to &nvar; %scan(&lis_var,&i) 
			lambda&i %end;;
		label %if &strato ne %then _strata_="Strata"; %do i=1 %to &nvar; 
			%scan(&lis_var,&i)="%scan(&lis_var,&i): value of test statistic"
			lambda&i="%scan(&lis_var,&i): p-value" %end;;
	run;
	proc rank data=pval_par(drop=&lis_var %do i=1 %to &nvar; lambda&i %end;) descending ties=high 
		out=pval_par(drop=comb1);
		%if &strato ne %then by _strata_;; 
 		var comb1;
  		ranks comb1_r;
	run;
	%if &strato ne %then %do;
		proc sort data=pval_par;
			by _sample_;
		run;
	%end;
	data pval_par(keep=_sample_ lambda %if &strato ne %then _strata_;) 
	 		%if &strato ne %then comb(keep=_sample_ comb_tot);;
		set pval_par;
		lambda=(0.5+comb1_r-1)/(&nsample+1); 
		loglam=log(lambda); 
		if _sample_=0 %if &strato ne %then then output pval_par;;
		%if &strato ne %then %do;
			by _sample_;
			if first._sample_ then comb_to=0; 
			comb_to+loglam;
			comb_tot=comb_to*(-2);
			if last._sample_ then output comb;
		%end;
	run;
	proc print data=pval_par noobs label;
		title %if &strato ne %then "Second-order Combined Tests"; %else "Overall Combined Test";;
		var %if &strato ne %then _strata_; lambda;
		label %if &strato ne %then _strata_="Strata"; lambda="p-value";
	run;
	%if &strato ne %then %do;
		proc rank data=comb descending ties=high out=comb(keep=_sample_ combtot);
 			var comb_tot;
  			ranks combtot;
		run;
		data comb;
			set comb;
			if _sample_=0;
			lambdat=(0.5+combtot-1)/(&nsample+1);
		run;
		proc print data=comb noobs label;
			title "Overall Combined Test";
			var lambdat;
			label lambdat="p-value";
		run;	
	%end;
	
%mend npc;