/*******************************************************************************
Program name: MI_SLOGIT.sas   
 
Programmers: Qixuan Chen, Ph.D (qc2138@columbia.edu)
            Department of Biostatistics
			Mailman School of Public Health
			Columbia University 

			Xiaohui Jiang, M.S. (ljcindy@umich.edu)
			Department of Environmental Health Sciences
			University of Michigan School of Public Health

Date: September 3rd, 2009   
 
Purpose: (1) Use MIANALYZE to combine the multiple regression coefficients
        estimates obtained using PROC SURVEYLOGISTIC.In the data set "&dset", 
		there should be one variable called "_imputation_".

		(2) The user need to specify the following macro variables:
		"dset"		-- source data
		"outcome"	-- outcome of the model
		"var"	    -- the list of covariates
		"catVar"-- the string of categorical variables in "varlist"
		"strata"	-- strata in a survey sample
		"cluster"	-- cluster in a survey sample
		"weight"	-- sample weight
		"output1"   -- the output of the combined beta coefficient estimates
		"output2"   -- the output of the combined p-values (useful when 
                       there are categorical variables)
********************************************************************************/


%macro MI_SLOGIT(dset, outcome, var, catVar, strata, cluster, weight, output1, output2);
	proc surveylogistic data=&dset NAMELEN=40;
		weight &weight;
		strata &strata;
		cluster &cluster;
		class &catVar; 
		model &outcome(event='1')= &var/RSQUARE covb; 
		by _imputation_;
		ods output ParameterEstimates=out ModelInfo=design covb=covariance;
	run;

	data _null_(where=(_imputation_=1)); set design;
		if Description="Number of Strata" then call symput('nstrata',Value);
		if Description="Number of Clusters" then call symput('ncluster',Value);
	run;
	%let df=%eval(&ncluster-&nstrata);

	%if &catVar=   %then %do;
	    data out1; set out;
			Parameter=trim(Variable);
			keep _imputation_ Parameter Estimate StdErr;
		run;
	%end;
	%else %do;
		data out1; set out;
			Parameter=trim(Variable)||ClassVal0;
			keep _imputation_ Parameter Estimate StdErr;
		run;
	%end;

	data testvar; set out; where _imputation_=1 and Variable ne 'Intercept'; run;

	proc freq data=testvar; tables Variable; ods output OneWayFreqs =freq; run;
	data freq; set freq; keep Variable frequency; run;
	proc sort data=freq;  by frequency; run;
	data _null_; set freq;
       if _N_=1 then call symput('min',frequency);
	run;
    proc sort data=freq;  by descending frequency; run;
	   data _null_; set freq;
       if _N_=1 then call symput('max',frequency);
	run;

	%if &max=1 %then %do;
	/* remember here use data=out */
		proc mianalyze parms=out covb(effectvar=stacking)=covariance edf=&df;
			modeleffects Intercept &var;
			ods output ParameterEstimates=result;
		run;
		data result; set result;
			if Probt=. then Probt=(1-probt(abs(Estimate/sqrt(StdErr)),&df))*2;
		run;
		data result; set result;
			if Probt=<0.0001 then Probt=0;
		run;
		data test; set result;
			format var $40.;
			if Parm="Intercept" then delete;
			var=Parm;
			Prob=round(Probt,0.001);
			keep var Prob;
	    run;
		proc datasets library=work;
			delete out out1 testvar freq covariance;
		run;
		quit;
	%end;

	%else %do;
	   %if &min=1 %then %do;
			data continuous; set freq; if frequency=1; drop frequency; run;
			proc sql noprint;
				select count(*) into :ncont from continuous;
				%let ncont=&ncont;
				select variable into : cont1-:cont&ncont from continuous;
			quit;
		%end;
		data categorical; set freq; if frequency>1; run;
		proc sql noprint;
			select count(*) into :nCate from categorical;
		  	%let nCate=&nCate; 
		  	select variable into :cate1-:cate&nCate from categorical;
		  	select frequency into :catFreq1-:catFreq&ncate from categorical;
		quit;

		data MIvar; set out1; where _imputation_=1; run;
		proc sql noprint;
			select Parameter into :MIvarlist separated by ' ' from MIvar;	
		quit;
		%put &MIvarlist;
		/* here use data=out1 */
		proc mianalyze parms=out1 covb(effectvar=stacking)=covariance edf=&df;
			modeleffects &MIvarlist;
			%if &min=1 %then %do;
				%do i=1 %to &nCont;
					&&cont&i : test &&cont&i;
				%end;
			%end;
			%do k=1 %to &nCate;
				&&cate&k : test %do j=1 %to %eval(&&catFreq&k-1);
								  	%str(&&cate&k)&j.=
								%end;
								%str(&&cate&k)%eval(&&catFreq&k)=0 /mult;
			%end;
			ods output TestMultStat=test ParameterEstimates=result;
		run;

		data test; set test;
			rename Test=Parm;
			if ProbF=. then ProbF=1-probf(FValue,NumDF,&df);
		run;
		data test(keep=Parm Prob); set test;
			rename ProbF=Prob;
		run;

		%if &min=1 %then %do;
			data continuous; set continuous; 
				format Parm $40.;
				index=1;
				Parm=variable; 
				drop variable;
			run;
			data result2(keep=Parm Prob); set result; 
				if tValue=. then Probt=(1-Probt(abs(Estimate/sqrt(StdErr)),&df))*2; 
				rename Probt=Prob; 
			run;
			proc sort data=continuous; by Parm; run;
			proc sort data=result2; by Parm; run;
			data result2; merge continuous result2; by Parm; run;
			data result2; set result2; if index=1; run;      
			data test; set test result2;run;
		%end;

		data test; set test;
			format var $40.;
			if Prob=<.0001 then Prob=0;
			var=Parm;
			keep var Prob;
		run;

	    proc datasets library=work;
	    	delete CATEGORICAL CONTINUOUS COVARIANCE DESIGN FREQ MIVAR OUT OUT1 TESTVAR result2;
		run; quit;
	%end;

	data &output1; set result; run;
	data &output2; set test; run;

%mend MI_SLOGIT;