/*
History:
2026-02-04 Update to allow for two-way interactions. 
2026-01-05 Update macro to allow for categorical covariates.
           The options cat_covariates, cat_reference, and cat_maxn have been added.

Significance Testing for MCMC Covariates

Description:

	Tests the significance of covariates in an MCMC model using replication techniques.
	
Parameters:

	- test_covariates: Space-delimited list of covariates to be tested. 
	- base_weight: Variable with weighting for each subject.
	- rep_weights: Space-delimited list of variables with the replicate weights for each subject.
								 Can be specified using DATA step variable list syntax.
	- brr: Flag specifying whether replicate weights are for balanced repeated replication (BRR). (Y/N, default = N)
	- fay_factor: Numeric between 0 and 1 specifying the Fay factor used to generate the BRR weights used in the function.
								Applies only when brr is Y. (default = 1)
	- df: Integer specifying the number of degrees of freedom for Hotelling's T^2 test.
	- initial_mcmc_seed: Numeric starting seed for the random number generator. 
											 If specified, the seed will be incremented for each MCMC replicate. 
											 If blank, uses a randomly generated integer from -10^7 to 10^7, exclusive for each replicate. (default = blank)
	- outlib: The SAS library to store output datasets. (default = WORK)
	- outname: The name of the output dataset.
	- cat_covariates: Space-delimited list of names of categorical covariates. Any character variable
	                  in covariates will be treated as categorical. For each categorical covariate,
					  dummy variables will be created. See cat_reference below for specifying
					  the reference categories.
	- cat_reference:  Space-delimited list of reference categories for the categorical variables
                      in cat_covariates. The length and order of cat_reference must match the 
                      length and order of cat_covariates.					  
	- cat_maxn: The maximum number of categories for a categorical variable. The default is 25.
	- DEBUG:  0 or 1 for debugging purposes. The default is 0.
	- v143code: Flag to run alternative code for IML version 14.3. (0/1, default is determined at run time)
	
	See %nci_multivar_mcmc() for other parameters.
	
Output:

	The following SAS data set is created in outlib:
	
		outname: A SAS data set that contains the joint p-value for all of the covariates in test_covariates.
	
Details:

	The variance/covariance matrix of the beta parameters in the MCMC
	model is found by fitting MCMC model replicates with bootstrap or balanced
	repeated replication (BRR). The original weighting variable for the dataset
	is given in base_weight and the replicate weighting variables are given
	in rep_weights. For survey datasets that BRR is applicable to, BRR
	weights can be generated by brr_weights. When using BRR, the value of
	fay_factor must match the Fay factor used to generate the weights.

	The covariates to test are specified in test_covariates. Elements of
	test_covariates that are not present in the model are ignored.
	The test_covariates can include covariates from the episodic indicator,
	episodic amount, and daily amount parts of the model. This function does
	not test never-consumer covariates.

	The df parameter is the number of degrees of freedom to use for the
	Hotelling's t-squared test. For survey datasets, this is the total number
	of PSUs across all strata minus the number of strata.
*/
%macro covariate_test(test_covariates=,
											base_weight=,
											rep_weights=,
											brr=N,
											fay_factor=1,
											df=,
											initial_mcmc_seed=,
											pre_mcmc_lib=WORK,
											pre_mcmc_data=,
											id=,
											repeat_obs=,
											episodic_variables=,
											episodic_indicators=,
											episodic_amounts=,
											daily_variables=,
											daily_amounts=,
											never_consumer_variable=,
											never_consumer_indicator=,
											never_consumer_amount=,
											default_covariates=,
											default_intercept=Y,
											episodic_indicator_covariates=,
											episodic_indicator_intercept=,
											episodic_amount_covariates=,
											episodic_amount_intercept=,
											daily_amount_covariates=,
											daily_amount_intercept=,
											individual_covariates=,
											individual_intercept=,
											never_consumer_covariates=,
											never_consumer_intercept=Y,
											num_mcmc_iterations=12000,
											num_burn=2000,
											num_thin=25,
											sigma_u_prior=,
											sigma_u_constant=N,
											low_memory_mode=N,
											outlib=WORK,
											outname=,
											cat_covariates=,
											cat_reference=,
											cat_maxn=25,
											DEBUG=0,
											v143code=);
											
	%local mcmc_seed rep var_j dat allv catvdata intvdata;
				 
	**Set the value for v143code if needed;
	%let v143code = %left(%trim(&v143code));
	%if %length(&v143code) = 0 %then %do;
		%store_IML_version(_tmp_dat_);
		%get_v143code(_tmp_dat_);
		%let v143code = 0;
	%end;	
	
	** Check for errors;
	%let rep_weights     = %normStr(&rep_weights);
	%let test_covariates = %normStr(&test_covariates);
	%let df              = %normStr(&df);
	%if %eval(%length(&rep_weights) = 0) %then %errorExit(9, rep_weights must be specified, exit=1);
	%if %eval(%length(&test_covariates) = 0) %then %errorExit(9, test_covariates must be specified, exit=1);
	%if %eval(%length(&df) = 0) %then %errorExit(9, df must be specified, exit=1);
	
	** Create dummy variables for categorical covariates;
	%local dat allv catvdata;
	%let dat      = &pre_mcmc_lib..&pre_mcmc_data._mcmc_in; /* input data */
	%let catvdata = _catv_info_tmp; /* output data of categorical variables */
	%let allv     = &test_covariates &default_covariates &episodic_indicator_covariates &episodic_amount_covariates;
	%let allv     = &allv &daily_amount_covariates &never_consumer_covariates;
	%let allv     = %normStr(&allv, upcase=1, removeMultBlanks=1);
	%createDummyVariables(&catvdata, &dat, &cat_covariates, refcats=&cat_reference, 
	                      allvars=&allv, maxncats=&cat_maxn, suffix=_CAT, DEBUG=&DEBUG);		
	**Update covariate lists for categorical covariates;	
	%let test_covariates               = %updateListWithDummyVars(&test_covariates,               &catvdata, suffix=_CAT);	
	%let default_covariates            = %updateListWithDummyVars(&default_covariates,            &catvdata, suffix=_CAT);	
	%let episodic_indicator_covariates = %updateListWithDummyVars(&episodic_indicator_covariates, &catvdata, suffix=_CAT);		
	%let episodic_amount_covariates    = %updateListWithDummyVars(&episodic_amount_covariates,    &catvdata, suffix=_CAT);	
	%let daily_amount_covariates       = %updateListWithDummyVars(&daily_amount_covariates,       &catvdata, suffix=_CAT);	
	%let never_consumer_covariates     = %updateListWithDummyVars(&never_consumer_covariates,     &catvdata, suffix=_CAT);	
	%updateDataWithDummyVars(&individual_covariates, &catvdata, var="covariates", suffix=_CAT);	
	
	**Update covariates list for interactions;
	%let intvdata = _intv_info_tmp; /* output data of interaction variables */
	%createInteractions(&intvdata, &dat, &allv, &catvdata, sep=X);
	%let test_covariates               = %updateListWithInteractions(&test_covariates,               &intvdata);
	%let default_covariates            = %updateListWithInteractions(&default_covariates,            &intvdata);
	%let episodic_indicator_covariates = %updateListWithInteractions(&episodic_indicator_covariates, &intvdata);
	%let episodic_amount_covariates    = %updateListWithInteractions(&episodic_amount_covariates,    &intvdata);
	%let daily_amount_covariates       = %updateListWithInteractions(&daily_amount_covariates,       &intvdata);
	%let never_consumer_covariates     = %updateListWithInteractions(&never_consumer_covariates,     &intvdata);

	**process variable and covariate lists;
	%check_and_process_input(pre_mcmc_data=&pre_mcmc_lib..&pre_mcmc_data.,
													 id=&id.,
													 weight=&base_weight.,
													 repeat_obs=&repeat_obs.,
													 episodic_variables=&episodic_variables.,
													 episodic_indicators=&episodic_indicators.,
													 episodic_amounts=&episodic_amounts.,
													 daily_variables=&daily_variables.,
													 daily_amounts=&daily_amounts.,
													 never_consumer_variable=&never_consumer_variable.,
													 never_consumer_indicator=&never_consumer_indicator.,
													 never_consumer_amount=&never_consumer_amount.,
													 default_covariates=&default_covariates.,
													 default_intercept=&default_intercept.,
													 episodic_indicator_covariates=&episodic_indicator_covariates.,
													 episodic_indicator_intercept=&episodic_indicator_intercept.,
													 episodic_amount_covariates=&episodic_amount_covariates.,
													 episodic_amount_intercept=&episodic_amount_intercept.,
													 daily_amount_covariates=&daily_amount_covariates.,
													 daily_amount_intercept=&daily_amount_intercept.,
													 individual_covariates=&individual_covariates.,
													 individual_intercept=&individual_intercept.,
													 never_consumer_covariates=&never_consumer_covariates.,
													 never_consumer_intercept=&never_consumer_intercept.,
													 sigma_u_prior=&sigma_u_prior.,
													 num_mcmc_iterations=&num_mcmc_iterations.,
													 num_burn=&num_burn.,
													 num_thin=&num_thin.,
													 num_post=0,
													 save_mcmc_u=N,
													 save_all_u=N,
													 outname=_init);
													 
	data _NULL_;
		set _init_vars end=last;
				
		call symputx("num_variables", num_variables, "L");
				
		num_test = countw("&test_covariates.", " ");
		call symputx("num_test", num_test, "L");
	run;
			
	data _NULL_;
		set _init_covars end=last;
				
		retain max_num_covariates 0;
				
		if num_covariates > max_num_covariates then max_num_covariates = num_covariates;
				
		if last = 1 then call symputx("max_num_covariates", max_num_covariates, "L");
	run;
	
	**select beta coefficients to test;
	data _NULL_;
		set _init_covars end=last;
		
		length selected_j1-selected_j&num_test. $32.;
		
		retain num_selected 0;
					 
		array selected_j{&num_test.} selected_j1-selected_j&num_test.;
		
		array covariates{&max_num_covariates.} covariate1-covariate&max_num_covariates.;
		
		do test_j = 1 to &num_test.;
		
			test_coefficient_j = scan("&test_covariates.", test_j, " ");
			
			do covar_j = 1 to num_covariates;
			
				if compare(test_coefficient_j, covariates{covar_j}, "i") = 0 then do;
				
					num_selected = num_selected + 1;
					selected_j{test_j} = covariates{covar_j};
				end;
			end;
		end;
		
		call symputx(cats("selected", _N_), catx(" ", of selected_j{*}), "L");
		
		if last = 1 then call symputx("num_selected", num_selected, "L");
	run;
	
	**check that enough degrees of freedom are present for selected coefficients;
	%if &num_selected. >= &df. %then %do;
	
		data _NULL_;
			
			put "ERROR: Insufficient degrees of freedom to test selected coefficients";
			abort abend 2;
		run;
	%end;											 
	
	**expand list of replicate weights;
	data _NULL_;
		set &pre_mcmc_lib..&pre_mcmc_data._mcmc_in (keep = &rep_weights.);
		
		array rep_weights{*} &rep_weights.;
		
		call symputx("num_replicates", dim(rep_weights), "L");
	run;
	
	data _weights (keep = _base _rep1-_rep&num_replicates.);
		set &pre_mcmc_lib..&pre_mcmc_data._mcmc_in (keep = &base_weight. &rep_weights.);
		
		length _base _rep1-_rep&num_replicates. $32.;
		
		array rep_weights{&num_replicates.} &rep_weights.;
		array rep_names{&num_replicates.} _rep1-_rep&num_replicates.;
		
		_base = vname(&base_weight.);
		
		do rep = 1 to &num_replicates.;
		
			rep_names{rep} = vname(rep_weights{rep});
		end;
	run;
											
	**run MCMC with base weights and all replicate weights and extract chosen beta parameters;
	%let mcmc_seed = &initial_mcmc_seed.;
	
	%do rep = 0 %to &num_replicates.;
	
		data _NULL_;
			set _weights;
			
			%if &rep. = 0 %then %do;
			
				call symputx("weight", _base, "L");
			%end;
			%else %do;
			
				call symputx("weight", _rep&rep., "L");
			%end;
		run;
		
		**Run MCMC;
		%nci_multivar_mcmc(pre_mcmc_lib=&pre_mcmc_lib.,
											 pre_mcmc_data=&pre_mcmc_data.,
											 id=&id.,
											 repeat_obs=&repeat_obs.,
											 weight=&weight.,
											 episodic_variables=&episodic_variables.,
											 episodic_indicators=&episodic_indicators.,
											 episodic_amounts=&episodic_amounts.,
											 daily_variables=&daily_variables.,
											 daily_amounts=&daily_amounts.,
											 default_covariates=&default_covariates.,
											 episodic_indicator_covariates=&episodic_indicator_covariates.,
											 episodic_amount_covariates=&episodic_amount_covariates.,
											 daily_amount_covariates=&daily_amount_covariates.,
											 individual_covariates=&individual_covariates.,
											 default_intercept=&default_intercept.,
											 episodic_indicator_intercept=&episodic_indicator_intercept.,
											 episodic_amount_intercept=&episodic_amount_intercept.,
											 daily_amount_intercept=&daily_amount_intercept.,
											 individual_intercept=&individual_intercept.,
											 never_consumer_variable=&never_consumer_variable.,
											 never_consumer_indicator=&never_consumer_indicator.,
											 never_consumer_amount=&never_consumer_amount.,
											 never_consumer_covariates=&never_consumer_covariates.,
											 never_consumer_intercept=&never_consumer_intercept.,
											 mcmc_seed=&mcmc_seed.,
											 num_mcmc_iterations=&num_mcmc_iterations.,
											 num_burn=&num_burn.,
											 num_thin=&num_thin.,
											 sigma_u_prior=&sigma_u_prior.,
											 sigma_u_constant=&sigma_u_constant.,
											 low_memory_mode=&low_memory_mode.,
											 outname=_rep,
											 DEBUG=&DEBUG,
											 v143code=&v143code);
		
		**Extract selected beta parameter means;
		data _beta_selected&rep.;
			%do var_j = 1 %to &num_variables.;
			
				set _rep_beta&var_j._mean (keep = &&selected&var_j..);
			%end;
		run;
		
		**increment seed if it is specified;
		%if &mcmc_seed. ^= %str() %then %do;
		
			%let mcmc_seed = %eval(&mcmc_seed. + 1);
		%end;
	%end;
	
	data _beta_boot;
		set _beta_selected0-_beta_selected&num_replicates.;
	run;
	
	proc iml;
	
		num_replicates = &num_replicates.;
	
		use _beta_boot;
		
			read point(1) var _ALL_ into beta_point;
			beta_point = beta_point`;
			
			read point(2:(num_replicates + 1)) var _ALL_ into beta_reps;
			beta_reps = beta_reps`;
		close _beta_boot;
	
		**calculate variance-covariance matrix of beta estimates;
		if upcase("&brr.") = "N" then fay_factor = 1;
		else fay_factor = &fay_factor.;
		
		beta_diff = beta_reps - beta_point;
		beta_sum_sq_diff = beta_diff * beta_diff`;
		beta_covariance = beta_sum_sq_diff/(num_replicates * fay_factor**2);
		
		** Check for singular matrix;
		p_value = .;
		rank = round(trace(ginv(beta_covariance)*beta_covariance));
		
		**Hotelling T^2 test;
		if rank >= nrow(beta_covariance) then do;
			t2 = beta_point` * inv(beta_covariance) * beta_point;
		
			N = &df.;
			p = &num_selected.;
		
			F_statistic = t2*(N-p)/(p*(N-1));
			p_value = 1 - probf(F_statistic, p, N-p);
		end;
		
		create &outlib..&outname. from p_value[colname={"p"}];
			append from p_value;
		close &outlib..&outname.;
	quit;
%mend covariate_test;