/*
Copy an MCMC Model

Description:

	Copy all datasets associated with an MCMC model.
	
Parameters:

	- multivar_mcmc_lib: The SAS library that contains the output datasets of an MCMC model. (default = WORK)
	- multivar_mcmc_model: Prefix for the output datasets of the MCMC model (max 23 characters).
												 See nci_multivar_mcmc macro for more details.
	- outlib: The SAS library to copy the MCMC model datasets into. (default = WORK)
	- outname: The prefix used to name the copied MCMC model datasets (max 23 characters).
						 Replaces the prefix of the input model.
	
Output:

	The MCMC model datasets given in the input are copied to outlib and prefixed with outname.
	
Details:

	The MCMC model output from %nci_multivar_mcmc() is a collection of SAS datasets.
	This macro copies all of the associated datasets and ensures that they are named correctly.
	This is useful when making modifications to a model while preserving a copy of the original.
*/
%macro copy_mcmc_model(multivar_mcmc_lib=WORK,
											 multivar_mcmc_model=,
											 outlib=WORK,
											 outname=);
													
	**Extract macro variables;
	data _NULL_;
		set &multivar_mcmc_lib..&multivar_mcmc_model._vars;
		
		num_variables = 2*num_episodic + num_daily;
		call symputx("num_variables", num_variables, "L");
		
		call symputx("has_never_consumers", has_never_consumers, "L");
	run;
	
	data _NULL_;
		set &multivar_mcmc_lib..&multivar_mcmc_model._iters;
		
		call symputx("has_mcmc_u", has_saved_mcmc_u, "L");
		
		has_post_u = (num_post > 0);
		call symputx("has_post_u", has_post_u, "L");
	run;
	
	**Copy datasets;
	
	**beta;
	%do var_j = 1 %to &num_variables.;
	
		data &outlib..&outname._beta&var_j._mean;
			set &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._mean;
		run;
		
		data &outlib..&outname._beta&var_j._var;
			set &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._var;
		run;
	
		data &outlib..&outname._beta&var_j.;
			set &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j.;
		run;
	%end;
	
	**sigma-u;
	data &outlib..&outname._sigma_u_mean;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_mean;
	run;
	
	data &outlib..&outname._sigma_u_var;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_var;
	run;
	
	data &outlib..&outname._sigma_u;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u;
	run;
	
	**sigma-e;
	data &outlib..&outname._sigma_e_mean;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_mean;
	run;
	
	data &outlib..&outname._sigma_e_var;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_var;
	run;
	
	data &outlib..&outname._sigma_e;
		set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e;
	run;
	
	**MCMC U matrices;
	%if &has_mcmc_u. = 1 %then %do;
	
		data &outlib..&outname._u_mcmc;
			set &multivar_mcmc_lib..&multivar_mcmc_model._u_mcmc;
		run;
	%end;
	
	**post-MCMC conditional U matrices;
	%if &has_post_u. = 1 %then %do;
	
		data &outlib..&outname._u_post;
			set &multivar_mcmc_lib..&multivar_mcmc_model._u_post;
		run;
	%end;
	
	**variables;
	data &outlib..&outname._vars;
		set &multivar_mcmc_lib..&multivar_mcmc_model._vars;
	run;
	
	**covariates;
	data &outlib..&outname._covars;
		set &multivar_mcmc_lib..&multivar_mcmc_model._covars;
	run;
	
	**subjects and weights;
	data &outlib..&outname._subjects;
		set &multivar_mcmc_lib..&multivar_mcmc_model._subjects;
	run;
	
	**MCMC iterations, burn-in, thinning, and post-MCMC iterations;
	data &outlib..&outname._iters;
		set &multivar_mcmc_lib..&multivar_mcmc_model._iters;
	run;
	
	**backtransformation data;
	data &outlib..&outname._backtran;
		set &multivar_mcmc_lib..&multivar_mcmc_model._backtran;
	run;
	
	**random seed;
	data &outlib..&outname._seed;
		set &multivar_mcmc_lib..&multivar_mcmc_model._seed;
	run;
	
	%if &has_never_consumers. = 1 %then %do;
	
		**alpha1;
		data &outlib..&outname._alpha1_mean;
			set &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_mean;
		run;
		
		data &outlib..&outname._alpha1_var;
			set &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_var;
		run;
	
		data &outlib..&outname._alpha1;
			set &multivar_mcmc_lib..&multivar_mcmc_model._alpha1;
		run;
		
		**consumer probabilities;
		data &outlib..&outname._conprob1_mean;
			set &multivar_mcmc_lib..&multivar_mcmc_model._conprob1_mean;
		run;
		
		data &outlib..&outname._conprob1_var;
			set &multivar_mcmc_lib..&multivar_mcmc_model._conprob1_var;
		run;
		
		data &outlib..&outname._conprob1;
			set &multivar_mcmc_lib..&multivar_mcmc_model._conprob1;
		run;
		
		**never-consumer covariates;
		data &outlib..&outname._covarsnc;
			set &multivar_mcmc_lib..&multivar_mcmc_model._covarsnc;
		run;
	%end;
	
	**Output listing;
	data &outlib..&outname._obj;
		set &multivar_mcmc_lib..&multivar_mcmc_model._obj;
		
		dataset = cats("&outname._", name);
	run;
%mend copy_mcmc_model;

/*
Extract NCI Multivar MCMC Parameters

Description:
	
	Extract posterior means and covariances from an NCI multivar MCMC model.
	
Parameters:

	- multivar_mcmc_lib: The SAS library that contains the output datasets of an MCMC model. (default = WORK)
	- multivar_mcmc_model: Prefix for the output datasets of the MCMC model (max 23 characters).
												 See nci_multivar_mcmc macro for more details.
	- do_mean: Y/N indicator. If Y, calculate posterior means of MCMC parameters
	- do_variance: Y/N indicator. If Y, calculate posterior variances of MCMC parameters
	- do_covariance: Y/N indicator. If Y, calculate posterior covariance matrix of MCMC parameters.
	- outlib: The SAS library to store output datasets. (default = WORK)
	- outname: The prefix used to name output datasets (max 26 characters).
	
Output:

	The following SAS data sets are created in outlib and prefixed with outname:
	
		If do_mean is Y:
		- _means: A SAS data set of posterior means of MCMC parameters.
		
		If do_variance is Y:
		- _var: A SAS data set of poster variances of MCMC parameters.
		
		If do_covariance is Y:
		- _cov: A SAS data set of posterior covariances of MCMC parameters.
								 
Details:

	This function calculates the posterior means and variances of
	parameters using the MCMC samples in an NCI multivar MCMC model. Means
	and covariances are calculated for the fixed effect coefficients (beta),
	random effect covariance (Sigma-u), and residual error covariance
	(Sigma-e) for all models. For models allowing never-consumers, the
	never-consumer model coefficients (alpha1) and mean consumer probability
	are also included. The calculations use samples after the burn-in period,
	spaced apart by the thinning number. It is important to note that the
	variances and covariaces returned by this function represent the variance
	of the Monte Carlo process, and it is distinct from the covariances that
	would be found from replication techniques such as bootstrap or BRR.
*/
%macro extract_parameters(multivar_mcmc_lib=WORK,
													multivar_mcmc_model=,
													do_mean=Y,
													do_variance=Y,
													do_covariance=Y,
													outlib=WORK,
													outname=);
													
	%local var_j;
	
	**number of variables in the model;
	data _NULL_;
		set &multivar_mcmc_lib..&multivar_mcmc_model._vars;
		
		call symputx("num_variables", num_variables, "L");
		call symputx("has_never_consumers", has_never_consumers, "L");
	run;
	
	**Calculate statistics;
	
	**Means;
	%if %upcase(&do_mean.) = Y %then %do;
	
		**beta;
		%do var_j = 1 %to &num_variables.;
		
			proc iml;
				
				use &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._mean;
					read all var _NUM_ into beta_j_mean;
				close &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._mean;
				
				beta_j_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._beta&var_j.");
				
				create _beta&var_j._mean from beta_j_mean[colname=beta_j_names];
					append from beta_j_mean;
				close _beta&var_j._mean;
			quit;
		%end;
		
		**sigma-u;
		proc iml;
		
			use &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_mean;
				read all var _NUM_ into sigma_u_mean;
				sigma_u_mean = symsqr(sigma_u_mean)`;
			close &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_mean;
			
			sigma_u_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._sigma_u");
			
			create _sigma_u_mean from sigma_u_mean[colname=sigma_u_names];
				append from sigma_u_mean;
			close _sigma_u_mean;
		quit;
		
		**sigma-e;
		proc iml;
		
			use &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_mean;
				read all var _NUM_ into sigma_e_mean;
				sigma_e_mean = symsqr(sigma_e_mean)`;
			close &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_mean;
			
			sigma_e_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._sigma_e");
			
			create _sigma_e_mean from sigma_e_mean[colname=sigma_e_names];
				append from sigma_e_mean;
			close _sigma_e_mean;
		quit;
		
		**never-consumer parameters;
		%if &has_never_consumers. = 1 %then %do;
		
			**alpha1;
			proc iml;
			
				use &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_mean;
				 	read all var _NUM_ into alpha1_mean;
				close &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_mean;
				 
				alpha1_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._alpha1");
				 
				create _alpha1_mean from alpha1_mean[colname=alpha1_names];
				 	append from alpha1_mean;
				close _alpha1_mean;
			quit;
			
			**consumer probabilities;
			proc iml;
			
				use &multivar_mcmc_lib..&multivar_mcmc_model._conprob1_mean;
					read all var _NUM_ into conprob1_mean;
				close &multivar_mcmc_lib..&multivar_mcmc_model._conprob1_mean;
				
				conprob1_name = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._conprob1");
				
				create _conprob1_mean from conprob1_mean[colname=conprob1_name];
					append from conprob1_mean;
				close _conprob1_mean;
			quit;
		%end;
		
		data &outlib..&outname._means;
			
			%do var_j = 1 %to &num_variables.;
				
				set _beta&var_j._mean;
			%end;
			
			set _sigma_u_mean;
			set _sigma_e_mean;
			
			%if &has_never_consumers. = 1 %then %do;
			
				set _alpha1_mean;
				set _conprob1_mean;
			%end;
		run;
	%end;
	
	**Variances;
	%if %upcase(&do_variance.) = Y %then %do;
	
		**beta;
		%do var_j = 1 %to &num_variables.;
		
			proc iml;
			
				use &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._var;
					read all var _NUM_ into beta_j_var;
				close &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j._var;
				
				beta_j_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._beta&var_j.");
				
				create _beta&var_j._var from beta_j_var[colname=beta_j_names];
					append from beta_j_var;
				close _beta&var_j._var;
			quit;
		%end;
		
		**sigma-u;
		proc iml;
		
			use &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_var;
				read all var _NUM_ into sigma_u_var;
				sigma_u_var = symsqr(sigma_u_var)`;
			close &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u_var;
			
			sigma_u_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._sigma_u");
			
			create _sigma_u_var from sigma_u_var[colname=sigma_u_names];
				append from sigma_u_var;
			close _sigma_u_var;
		quit;
		
		**sigma-e;
		proc iml;
		
			use &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_var;
				read all var _NUM_ into sigma_e_var;
				sigma_e_var = symsqr(sigma_e_var)`;
			close &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e_var;
			
			sigma_e_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._sigma_e");
			
			create _sigma_e_var from sigma_e_var[colname=sigma_e_names];
				append from sigma_e_var;
			close _sigma_e_var;
		quit;
		
		**never-consumer parameters;
		%if &has_never_consumers. = 1 %then %do;
		
			**alpha1;
			proc iml;
			
				use &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_var;
					read all var _NUM_ into alpha1_var;
				close &multivar_mcmc_lib..&multivar_mcmc_model._alpha1_var;
				
				alpha1_names = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._alpha1");
				
				create _alpha1_var from alpha1_var[colname=alpha1_names];
					append from alpha1_var;
				close _alpha1_var;
			quit;
			
			**consumer probabilities;
			proc iml;
			
				use &multivar_mcmc_lib..&multivar_mcmc_model._conprob1_var;
					read all var _NUM_ into conprob1_var;
				close &multivar_mcmclib..&multivar_mcmc_model._conprob1_var;
				
				conprob1_name = contents("&multivar_mcmc_lib.", "&multivar_mcmc_model._conprob1");
				
				create _conprob1_var from conprob1_var[colname=conprob1_name];
					append from conprob1_var;
				close _conprob1_var;
			quit;
		%end;
		
		data &outlib..&outname._var;
		
			%do var_j = 1 %to &num_variables.;
			
				set _beta&var_j._var;
			%end;
			
			set _sigma_u_var;
			set _sigma_e_var;
			
			%if &has_never_consumers. = 1 %then %do;
			
				set _alpha1_var;
				set _conprob1_var;
			%end;
		run;
	%end;
	
	**Covariance matrix;
	%if %upcase(&do_covariance.) = Y %then %do;
	
		**number of iterations and burn-in;
		data _NULL_;
			set &multivar_mcmc_lib..&multivar_mcmc_model._iters;
			
			call symputx("num_mcmc_iterations", num_mcmc_iterations, "L");
			call symputx("num_burn", num_burn, "L");
			call symputx("num_thin", num_thin, "L");
		run;
	
		**MCMC parameters at thinned iterations after burn-in;
		data _mcmc_samples;
		
			%do var_j = 1 %to &num_variables.;
			
				set &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j.;
			%end;
			
			set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u;
			set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e;
			
			%if &has_never_consumers. = 1 %then %do;
			
				set &multivar_mcmc_lib..&multivar_mcmc_model._alpha1;
				set &multivar_mcmc_lib..&multivar_mcmc_model._conprob1;
			%end;
			
			if _N_ > &num_burn. and
				 _N_ <= &num_mcmc_iterations. and
				 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
		run;
		
		proc corr data=_mcmc_samples cov outp=_mcmc_cov noprint;
		
			var _ALL_;
		run;
		
		data &outlib..&outname._cov (keep = _NAME_
		
																				beta: sigma_e: sigma_u:
																				
																				%if &has_never_consumers. = 1 %then %do;
																				
																					alpha1: conprob1
																				%end;
																);
			set _mcmc_cov;
			if _TYPE_ = "COV";
		run;
	%end;
%mend extract_parameters;

/*
Modify Burn-In and Thinning for MCMC Model

Description:

	This macro modifies the burn-in and thinning values of a
	multivar MCMC model that are used to calculate posterior means and
	covariances. This can be useful for determining the impact of different
	values and finding optimal numbers that balance accuracy and computation
	speed.
	
Parameters:

	- multivar_mcmc_lib: The SAS library that contains the output datasets of an MCMC model. (default = WORK)
	- multivar_mcmc_model: Prefix for the output datasets of the MCMC model (max 23 characters).
												 See nci_multivar_mcmc macro for more details.
	- num_mcmc_iterations: An integer specifying the new total number of iterations in the MCMC chain, including burn-in.
	- num_burn: An integer specifying the new number of burn-in iterations in the MCMC chain.
	- num_thin: An integer specifying the new number of iterations between MCMC samples used for calculating posterior means.
	- outlib: The SAS library that contains the datasets of the new MCMC model. If blank, uses the multivar_mcmc_lib. (default = WORK)
	- outname: Prefix for the datasets of the new MCMC model. If blank, uses multivar_mcmc_model. (default = blank)
	
Output:

	New MCMC model datasets are created in outlib and prefixed without outname.
	
	For more information, see %nci_multivar_mcmc().
	
Details:

	Since multivar MCMC models save MCMC samples from every
	iteration, the number of iterations, burn-in, and thinning values used to
	calculate posterior means and covariances can be changed. This utility can
	be used to update those values while performing error-checking.

	For models with stored random effect (u) matrices from the main MCMC
	chain, not all u matrices are stored by default. As a result, the burn-in
	for models with main MCMC chain u matrices cannot be decreased, and the
	new thinning value must be a multiple of the original thinning value. An
	error will be thrown if u matrices are not found for the new thinned
	iterations.

	Since changing the number of iterations, burn-in, and thinning can alter
	the posterior means, post-MCMC 'u' matrix draws conditional on the
	posterior mean may not be valid and will throw a warning.
*/
%macro burn_and_thin(multivar_mcmc_lib=WORK,
										 multivar_mcmc_model=,
										 num_mcmc_iterations=,
										 num_burn=,
										 num_thin=,
										 outlib=WORK,
										 outname=);	
										 
	%if &outlib. = %str() %then %do;
	
		%let outlib = &multivar_mcmc_lib.;
	%end;					
	
	%if &outname. = %str() %then %do;
	
		%let outname = &multivar_mcmc_model.;
	%end;
	
	%copy_mcmc_model(multivar_mcmc_lib=&multivar_mcmc_lib.,
									 multivar_mcmc_model=&multivar_mcmc_model.,
									 outlib=&outlib.,
									 outname=&outname.);
									 
	data &outlib..&outname._iters;
		set &multivar_mcmc_lib..&multivar_mcmc_model._iters;
		
		%if &num_mcmc_iterations. ^= %str() %then %do;
		
			num_mcmc_iterations = &num_mcmc_iterations.;
		%end;
		
		%if &num_burn. ^= %str() %then %do;
		
			num_burn = &num_burn.;
		%end;
		
		%if &num_thin. ^= %str() %then %do;
		
			num_thin = &num_thin.;
		%end;
		
		call symputx("num_traces", num_traces, "L");
		
		call symputx("num_mcmc_iterations", num_mcmc_iterations, "L");
		call symputx("num_burn", num_burn, "L");
		call symputx("num_thin", num_thin, "L");
		
		call symputx("has_saved_mcmc_u", has_saved_mcmc_u, "L");
		
		call symputx("num_post", num_post, "L");
	run;
	
	%if &num_mcmc_iterations. > &num_traces. %then %do;
	
		data _NULL_;
		
			put "ERROR: Number of MCMC iterations cannot be higher than the number of traces in the model";
			abort abend 2;
		run;
	%end;
	
	%if &num_burn. >= &num_mcmc_iterations. %then %do;
	
		data _NULL_;
		
			put "ERROR: Number of MCMC iterations must be greater than the number of burn-in iterations";
			abort abend 2;
		run;
	%end;
	
	%if &has_saved_mcmc_u. = 1 %then %do;
	
		data _NULL_;
			set &multivar_mcmc_lib..&multivar_mcmc_model._u_mcmc end=last_obs;
			by iteration;
			if iteration > &num_burn. and
				 iteration <= &num_mcmc_iterations. and
				 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
				 
			retain num_found;
			
			if _N_ = 1 then num_found = 0;
			
			if first.iteration then num_found = num_found + 1;
			
			if last_obs = 1 then do;
			
				num_expected = int((&num_mcmc_iterations. - (&num_burn. + 1))/&num_thin.) + 1;
				if num_found ^= num_expected then do;
				
					put "ERROR: MCMC chain U matrices not saved for requested iterations";
					abort abend 2;
				end;
			end;
		run;
	%end;
	
	%if &num_post. > 0 %then %do;
	
		data _NULL_;
		
			put "WARNING: Post-MCMC U matrices may not be valid because posterior means may have changed";
		run;
	%end;
	
	%if %sysfunc(exist(&outlib..&outname._ll)) = 1 %then %do;
	
		data _NULL_;
		
			put "WARNING: Log-likelihood may not be valid because posterior means may have changed";
		run;
	%end;
	
	**re-calculate posterior means and variances;
	data _NULL_;
		set &outlib..&outname._vars;
		
		call symputx("num_variables", num_variables, "L");
		
		call symputx("has_never_consumers", has_never_consumers, "L");
	run;
	
	**beta;
	%do var_j = 1 %to &num_variables.;
	
		data _beta&var_j._samples;
			set &outlib..&outname._beta&var_j.;
			if _N_ > &num_burn. and
				 _N_ <= &num_mcmc_iterations. and
				 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
		run;
		
		proc means data=_beta&var_j._samples noprint;
			
			var _ALL_;
			
			output out=_beta&var_j._mean_new (keep = beta&var_j.:) mean=;
		run;
		
		proc means data=_beta&var_j._samples noprint;
		
			var _ALL_;
			
			output out=_beta&var_j._var_new (keep = beta&var_j.:) var=;
		run;
		
		proc iml;
		
			use _beta&var_j._mean_new;
				read all var _ALL_ into beta_j_mean_new;
			close _beta&var_j._mean_new;
			
			use &outlib..&outname._beta&var_j._mean;
				read all var _NUM_ into beta_j_mean[rowname=_VARIABLE_ colname=coefficient_names_j];
			close &outlib..&outname._beta&var_j._mean;
			
			beta_j_mean = beta_j_mean_new;
			
			create &outlib..&outname._beta&var_j._mean from beta_j_mean[rowname=_VARIABLE_ colname=coefficient_names_j];
				append from beta_j_mean[rowname=_VARIABLE_];
			close &outlib..&outname._beta&var_j._mean;

		
			use _beta&var_j._var_new;
				read all var _ALL_ into beta_j_var_new;
			close _beta&var_j._var_new;
			
			use &outlib..&outname._beta&var_j._var;
				read all var _NUM_ into beta_j_var[rowname=_VARIABLE_ colname=coefficient_names_j];
			close &outlib..&outname._beta&var_j._var;
			
			beta_j_var = beta_j_var_new;
			
			create &outlib..&outname._beta&var_j._var from beta_j_var[rowname=_VARIABLE_ colname=coefficient_names_j];
				append from beta_j_var[rowname=_VARIABLE_];
			close &outlib..&outname._beta&var_j._var;
		quit;
	%end;
	
	**sigma-u;
	data _sigma_u_samples;
		set &outlib..&outname._sigma_u;
		if _N_ > &num_burn. and
			 _N_ <= &num_mcmc_iterations. and
			 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
	run;
	
	proc means data=_sigma_u_samples noprint;
	
		var _ALL_;
		
		output out=_sigma_u_mean_new (keep = sigma_u:) mean=;
	run;
	
	proc means data=_sigma_u_samples noprint;
		
		var _ALL_;
		
		output out=_sigma_u_var_new (keep = sigma_u:) var=;
	run;
	
	proc iml;
	
		use _sigma_u_mean_new;
			read all var _ALL_ into sigma_u_mean_new;
		close _sigma_u_mean_new;
		
		use &outlib..&outname._sigma_u_mean;
			read all var _NUM_ into sigma_u_mean[rowname=_VARIABLE_ colname=variable_names];
		close &outlib..&outname._sigma_u_mean;
		
		sigma_u_mean = sqrsym(sigma_u_mean_new`);
		
		create &outlib..&outname._sigma_u_mean from sigma_u_mean[rowname=_VARIABLE_ colname=variable_names];
			append from sigma_u_mean[rowname=_VARIABLE_];
		close &outlib..&outname._sigma_u_mean;
		
		
		use _sigma_u_var_new;
			read all var _ALL_ into sigma_u_var_new;
		close _sigma_u_var_new;
		
		use &outlib..&outname._sigma_u_var;
			read all var _NUM_ into sigma_u_var[rowname=_VARIABLE_ colname=variable_names];
		close &outlib..&outname._sigma_u_var;
		
		sigma_u_var = sqrsym(sigma_u_var_new`);
		
		create &outlib..&outname._sigma_u_var from sigma_u_var[rowname=_VARIABLE_ colname=variable_names];
			append from sigma_u_var[rowname=_VARIABLE_];
		close &outlib..&outname._sigma_u_var;
	quit;
	
	**sigma-e;
	data _sigma_e_samples;
		set &outlib..&outname._sigma_e;
		if _N_ > &num_burn. and
			 _N_ <= &num_mcmc_iterations. and
			 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
	run;
	
	proc means data=_sigma_e_samples noprint;
	
		var _ALL_;
		
		output out=_sigma_e_mean_new (keep = sigma_e:) mean=;
	run;
	
	proc means data=_sigma_e_samples noprint;
	
		var _ALL_;
		
		output out=_sigma_e_var_new (keep = sigma_e:) var=;
	run;
	
	proc iml;
	
		use _sigma_e_mean_new;
			read all var _ALL_ into sigma_e_mean_new;
		close _sigma_e_mean_new;
		
		use &outlib..&outname._sigma_e_mean;
			read all var _NUM_ into sigma_e_mean[rowname=_VARIABLE_ colname=variable_names];
		close &outlib..&outname._sigma_e_mean;
		
		sigma_e_mean = sqrsym(sigma_e_mean_new`);
		
		create &outlib..&outname._sigma_e_mean from sigma_e_mean[rowname=_VARIABLE_ colname=variable_names];
			append from sigma_e_mean[rowname=_VARIABLE_];
		close &outlib..&outname._sigma_e_mean;
		
		
		use _sigma_e_var_new;
			read all var _ALL_ into sigma_e_var_new;
		close _sigma_e_var_new;
		
		use &outlib..&outname._sigma_e_var;
			read all var _NUM_ into sigma_e_var[rowname=_VARIABLE_ colname=variable_names];
		close &outlib..&outname._sigma_e_var;
		
		sigma_e_var = sqrsym(sigma_e_var_new`);
		
		create &outlib..&outname._sigma_e_var from sigma_e_var[rowname=_VARIABLE_ colname=variable_names];
			append from sigma_e_var[rowname=_VARIABLE_];
		close &outlib..&outname._sigma_e_var;
	quit;
	
	**never-consumer parameters;
	%if &has_never_consumers. = 1 %then %do;
	
		**alpha1;
		data _alpha1_samples;
			set &outlib..&outname._alpha1;
			if _N_ > &num_burn. and
				 _N_ <= &num_mcmc_iterations. and
				 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
		run;
		
		proc means data=_alpha1_samples noprint;
		
			var _ALL_;
			
			output out=_alpha1_mean_new (keep = alpha1:) mean=;
		run;
		
		proc means data=_alpha1_samples noprint;
		
			var _ALL_;
			
			output out=_alpha1_var_new (keep = alpha1:) var=;
		run;
		
		proc iml;
		
			use _alpha1_mean_new;
				read all var _ALL_ into alpha1_mean_new;
			close _alpha1_mean_new;
			
			use &outlib..&outname._alpha1_mean;
				read all var _NUM_ into alpha1_mean[rowname=_VARIABLE_ colname=never_consumer_coefficient_names];
			close &outlib..&outname._alpha1_mean;
			
			alpha1_mean = alpha1_mean_new;
			
			create &outlib..&outname._alpha1_mean from alpha1_mean[rowname=_VARIABLE_ colname=never_consumer_coefficient_names];
				append from alpha1_mean[rowname=_VARIABLE_];
			close &outlib..&outname._alpha1_mean;
			
			
			use _alpha1_var_new;
				read all var _ALL_ into alpha1_var_new;
			close _alpha1_var_new;
			
			use &outlib..&outname._alpha1_var;
				read all var _NUM_ into alpha1_var[rowname=_VARIABLE_ colname=never_consumer_coefficient_names];
			close &outlib..&outname._alpha1_var;
			
			alpha1_var = alpha1_var_new;
			
			create &outlib..&outname._alpha1_var from alpha1_var[rowname=_VARIABLE_ colname=never_consumer_coefficient_names];
				append from alpha1_var[rowname=_VARIABLE_];
			close &outlib..&outname._alpha1_var;
		quit;
		
		**consumer probabilities;
		data _conprob1_samples;
			set &outlib..&outname._conprob1;
			if _N_ > &num_burn. and
				 _N_ <= &num_mcmc_iterations. and
				 mod(_N_, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
		run;
		
		proc means data=_conprob1_samples noprint;
		
			var _ALL_;
			
			output out=_conprob1_mean_new (keep = conprob1) mean=;
		run;
		
		proc means data=_conprob1_samples noprint;
		
			var _ALL_;
			
			output out=_conprob1_var_new (keep = conprob1) var=;
		run;
		
		proc iml;
		
			use _conprob1_mean_new;
				read all var _ALL_ into conprob1_mean_new;
			close _conprob1_mean_new;
			
			use &outlib..&outname._conprob1_mean;
				read all var _NUM_ into conprob1_mean[rowname=_VARIABLE_ colname={"prob"}];
			close &outlib..&outname._conprob1_mean;
			
			conprob1_mean = conprob1_mean_new;
			
			create &outlib..&outname._conprob1_mean from conprob1_mean[rowname=_VARIABLE_ colname={"prob"}];
				append from conprob1_mean[rowname=_VARIABLE_];
			close &outlib..&outname._conprob1_mean;
			
			
			use _conprob1_var_new;
				read all var _ALL_ into conprob1_var_new;
			close _conprob1_var_new;
			
			use &outlib..&outname._conprob1_var;
				read all var _NUM_ into conprob1_var[rowname=_VARIABLE_ colname={"prob"}];
			close &outlib..&outname._conprob1_var;
			
			conprob1_var = conprob1_var_new;
			
			create &outlib..&outname._conprob1_var from conprob1_var[rowname=_VARIABLE_ colname={"prob"}];
				append from conprob1_var[rowname=_VARIABLE_];
			close &outlib..&outname._conprob1_var;
		quit;
	%end;
%mend burn_and_thin;

/*
Draw Parameter Samples from an MCMC Chain

Description:

	This macro extracts samples of model parameters from an MCMC chain at a specified interval.
	
Parameters:

	- multivar_mcmc_lib: The SAS library that contains the output datasets of an MCMC model. (default = WORK)
	- multivar_mcmc_model: Prefix for the output datasets of the MCMC model (max 23 characters).
												 See nci_multivar_mcmc macro for more details.
	- num_draws:  integer specifying the number of parameter samples to draw from the MCMC chain. (default = 10)
	- draw_spacing: An integer that specifies that number of thinned samples between parameter draws. (default = 5)
	- outlib: The SAS library to store output datasets. (default = WORK)
	- outname: The prefix used to name output datasets (max 21 characters, default = DRAW).
	
Output:

	The following SAS data sets are created in outlib:
	
		outname1-outname&num_draws.: Multivar MCMC datasets for each parameter set that was drawn from multivar_mcmc_model.
																 For a complete listing of multivar MCMC datasets, see Output section for nci_multivar_mcmc.
	
Details:

	Drawing parameter values from an MCMC chain after the burn-in period
	is equivalent to drawing from the posterior distribution of the parameters. The
	variation in the MCMC chain is used to sample different parameter values as
	opposed to finding a covariance matrix. Using the MCMC chain has the
	advantage of using the actual posterior distribution of the parameters
	instead of making an assumption. Accurately sampling parameter values is
	required for measurement error correction methods such as multiple
	imputation.

	The function starts by drawing the parameter set from the first iteration
	after burn-in, and then draws parameter sets every draw_spacing
	iterations until num_draws sets have been drawn. Each set of parameters is 
	used to construct a new MCMC model that has only one iteration.
	
	If random effects (U) matrices from the MCMC chain are saved, the corresponding 
	U matrix for each drawn parameter set will be saved as a post-MCMC U matrix draw 
	in the new model.
*/
%macro draw_parameters(multivar_mcmc_lib=WORK,
											 multivar_mcmc_model=,
											 num_draws=10,
											 draw_spacing=5,
											 outlib=WORK,
											 outname=DRAW);
											 
	%local i;
										 
	data _NULL_;
		set &multivar_mcmc_lib..&multivar_mcmc_model._iters;
		
		call symputx("num_mcmc_iterations", num_mcmc_iterations, "L");
		call symputx("num_burn", num_burn, "L");
		call symputx("num_thin", num_thin, "L");
		
		call symputx("has_saved_mcmc_u", has_saved_mcmc_u, "L");
	run;
	
	data _NULL_;
		set &multivar_mcmc_lib..&multivar_mcmc_model._vars;
		
		call symputx("num_variables", num_variables, "L");
		call symputx("has_never_consumers", has_never_consumers, "L");
	run;
	
	**Check if enough iterations exist for parameter draws;
	data _NULL_;
	
		first_draw = &num_burn. + 1;
		last_draw = first_draw + &draw_spacing. * &num_thin. * (&num_draws. - 1);
		
		if last_draw > &num_mcmc_iterations. then do;
		
			put "ERROR: Not enough iterations after burn-in to make the requested number of parameter draws";
		end;
		
		call symputx("first_draw", first_draw, "L");
		call symputx("last_draw", last_draw, "L");
	run;
	
	%do i = 1 %to &num_draws.;
	
		data _NULL_;
		
			draw = &first_draw. + &draw_spacing. * &num_thin. * (&i. - 1);
			call symputx("draw", draw, "L");
		run;
		
		**draw parameters from specified iteration;
		%copy_mcmc_model(multivar_mcmc_lib=&multivar_mcmc_lib.,
										 multivar_mcmc_model=&multivar_mcmc_model.,
										 outlib=&outlib.,
										 outname=&outname.&i.);
		
		**beta;
		%do var_j = 1 %to &num_variables.;
		
			data &outlib..&outname.&i._beta&var_j.;
				set &multivar_mcmc_lib..&multivar_mcmc_model._beta&var_j.;
				if _N_ = &draw.;
			run;
			
			proc iml;
			
				use &outlib..&outname.&i._beta&var_j.;
					read all var _ALL_ into beta_j;
				close &outlib..&outname.&i._beta&var_j.;
				
				use &outlib..&outname.&i._beta&var_j._mean;
					read all var _NUM_ into beta_j_mean[rowname=_VARIABLE_ colname=beta_j_names];
				close &outlib..&outname.&i._beta&var_j._mean;
				
				beta_j_mean = beta_j;
				
				create &outlib..&outname.&i._beta&var_j._mean from beta_j_mean[rowname=_VARIABLE_ colname=beta_j_names];
					append from beta_j_mean[rowname=_VARIABLE_];
				close &outlib..&outname.&i._beta&var_j._mean;
			quit;
			
			data &outlib..&outname.&i._beta&var_j._var;
				set &outlib..&outname.&i._beta&var_j._var;
				
				call missing(of _NUMERIC_);
			run;
		%end;
		
		**sigma-u;
		data &outlib..&outname.&i._sigma_u;
			set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_u;
			if _N_ = &draw.;
		run;
		
		proc iml;
			
			use &outlib..&outname.&i._sigma_u;
				read all var _ALL_ into sigma_u;
			close &outlib..&outname.&i._sigma_u;
			
			use &outlib..&outname.&i._sigma_u_mean;
				read all var _NUM_ into sigma_u_mean[rowname=_VARIABLE_ colname=sigma_u_names]; 
			close &outlib..&outname.&i._sigma_u_mean;
			
			sigma_u_mean = sqrsym(sigma_u`);
			
			create &outlib..&outname.&i._sigma_u_mean from sigma_u_mean[rowname=_VARIABLE_ colname=sigma_u_names];
				append from sigma_u_mean[rowname=_VARIABLE_];
			close &outlib..&outname.&i._sigma_u_mean;
		quit;
		
		data &outlib..&outname.&i._sigma_u_var;
			set &outlib..&outname.&i._sigma_u_var;
			
			call missing(of _NUMERIC_);
		run;
		
		**sigma-e;
		data &outlib..&outname.&i._sigma_e;
			set &multivar_mcmc_lib..&multivar_mcmc_model._sigma_e;
			if _N_ = &draw.;
		run;
		
		proc iml;
		
			use &outlib..&outname.&i._sigma_e;
				read all var _ALL_ into sigma_e;
			close &outlib..&outname.&i._sigma_e;
			
			use &outlib..&outname.&i._sigma_e_mean;
				read all var _NUM_ into sigma_e_mean[rowname=_VARIABLE_ colname=sigma_e_names]; 
			close &outlib..&outname.&i._sigma_e_mean;
			
			sigma_e_mean = sqrsym(sigma_e`);
			
			create &outlib..&outname.&i._sigma_e_mean from sigma_e_mean[rowname=_VARIABLE_ colname=sigma_e_names];
				append from sigma_e_mean[rowname=_VARIABLE_];
			close &outlib..&outname.&i._sigma_e_mean;
		quit;
		
		data &outlib..&outname.&i._sigma_e_var;
			set &outlib..&outname.&i._sigma_e_var;
			
			call missing(of _NUMERIC_);
		run;
		
		**U matrix;
		%if &has_saved_mcmc_u. = 1 %then %do;
		
			data &outlib..&outname.&i._u_post (keep = post_mcmc_iteration u_col:);
				set &outlib..&outname.&i._u_mcmc;
				if iteration = &draw.;
				
				post_mcmc_iteration = 1;
			run;
			
			proc datasets library=&outlib.;
			
				delete &outname.&i._u_mcmc;
			run;
			
			data &outlib..&outname.&i._obj;
				set &outlib..&outname.&i._obj;
				
				if name = "u_mcmc" then delete;
			run;
		%end;
		%else %do;
		
			proc datasets library=&outlib.;
			
				delete &outname.&i._u_post;
			run;
			
			data &outlib..&outname.&i._obj;
				set &outlib..&outname.&i._obj;
				
				if name = "u_post" then delete;
			run;
		%end;
		
		**never-consumer parameters;
		%if &has_never_consumers. = 1 %then %do;
		
			**alpha1;
			data &outlib..&outname.&i._alpha1;
				set &multivar_mcmc_lib..&multivar_mcmc_model._alpha1;
				if _N_ = &draw.;
			run;
			
			proc iml;
			
				use &outlib..&outname.&i._alpha1;
					read all var _ALL_ into alpha1;
				close &outlib..&outname.&i._alpha1;
				
				use &outlib..&outname.&i._alpha1_mean;
					read all var _NUM_ into alpha1_mean[rowname=_VARIABLE_ colname=alpha1_names];
				close &outlib..&outname.&i._alpha1_mean;
				
				alpha1_mean = alpha1;
				
				create &outlib..&outname.&i._alpha1_mean from alpha1_mean[rowname=_VARIABLE_ colname=alpha1_names];
					append from alpha1_mean[rowname=_VARIABLE_];
				close &outlib..&outname.&i._alpha1_mean;
			quit;
			
			data &outlib..&outname.&i._alpha1_var;
				set &outlib..&outname.&i._alpha1_var;
				
				call missing(of _NUMERIC_);
			run;
			
			**consumer probabilities;
			data &outlib..&outname.&i._conprob1;
				set &multivar_mcmc_lib..&multivar_mcmc_model._conprob1;
				if _N_ = &draw.;
			run;
			
			proc iml;
			
				use &outlib..&outname.&i._conprob1;
					read all var _ALL_ into conprob1;
				close &outlib..&outname.&i._conprob1;
				
				use &outlib..&outname.&i._conprob1_mean;
					read all var _NUM_ into conprob1_mean[rowname=_VARIABLE_ colname=conprob1_names];
				close &outlib..&outname.&i._conprob1_mean;
				
				conprob1_mean = conprob1;
				
				create &outlib..&outname.&i._conprob1_mean from conprob1_mean[rowname=_VARIABLE_ colname=conprob1_names];
					append from conprob1_mean[rowname=_VARIABLE_];
				close &outlib..&outname.&i._conprob1_mean;
			quit;
			
			data &outlib..&outname.&i._conprob1_var;
				set &outlib..&outname.&i._conprob1_var;
				
				call missing(of _NUMERIC_);
			run;
		%end;
		
		data &outlib..&outname.&i._iters;
		
			num_traces = 1;
			
			num_mcmc_iterations = 1;
			num_burn = 0;
			num_thin = 1;
			num_post = &has_saved_mcmc_u.;
			
			has_saved_mcmc_u = 0;
		run;
	%end;
%mend draw_parameters;