**Produce trace plots of beta values and means;
%macro beta_trace_plots(beta=,
												variable_data=,
												covariate_data=,
												num_mcmc_iterations=,
												num_burn=,
												num_thin=,
												num_variables=,
												max_num_covariates=,
												max_num_coefficients=,
												plot_type=);
												
	%local var_j coef_j;
	
	%do var_j = 1 %to &num_variables.;
	
		**get variable;
		data _NULL_;
			set &variable_data.;
			
			call symputx("variable", variable&var_j., "L");
		run;
		
		**get coefficients;
		data _coefficients_j (keep = column1-column&max_num_coefficients. coefficient1-coefficient&max_num_coefficients.);
			set &covariate_data.;
			if _N_ = &var_j.;
			
			length column1-column&max_num_coefficients. $32.;
			length coefficient1-coefficient&max_num_coefficients. $32.;
			
			array columns{&max_num_coefficients.} column1-column&max_num_coefficients.;
			array coefficients{&max_num_coefficients.} coefficient1-coefficient&max_num_coefficients.;
			
			%if &max_num_covariates. > 0 %then %do;
			
				array covariates{&max_num_covariates.} covariate1-covariate&max_num_covariates.;
			%end;
			
			if intercept = 1 then do;
			
				columns{1} = "beta&var_j._intercept";
				coefficients{1} = "Intercept";
			end;
			
			do covar_j = 1 to num_covariates;
			
				columns{intercept + covar_j} = cats("beta&var_j._covariate", covar_j);
				coefficients{intercept + covar_j} = covariates{covar_j};
			end;
			
			call symputx("num_coefficients_j", num_coefficients);
		run;
		
		**make plots;
		data _value_plot_data;
			set &beta.&var_j.;
			
			iteration = _N_;
		run;
		
		%do coef_j = 1 %to &num_coefficients_j.;
		
			data _NULL_;
				set _coefficients_j;
				
				call symputx("column", column&coef_j., "L");
				call symputx("coefficient", coefficient&coef_j., "L");
			run;
			
			%if %sysfunc(findw(&plot_type., VALUE, , %str(s))) ^= 0 %then %do;
			
				proc sgplot data=_value_plot_data;
				
					series x=iteration y=&column. / lineattrs=(color=black);
					
					refline &num_burn. / axis=X lineattrs=(color=red);
					
					xaxis label="Number of Iterations";
					yaxis label="Beta Value";
					
					title "Beta Values, &variable., &coefficient.";
				run;
			%end;
			
			%if %sysfunc(findw(&plot_type., MEAN, , %str(s))) ^= 0 %then %do;
			
				data _mean_plot_data;
					set _value_plot_data;
					if iteration > &num_burn. and 
						 iteration <= &num_mcmc_iterations. and 
						 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
						 
					retain mean_value num_samples 0;
					
					num_samples = num_samples + 1;
					mean_value = mean_value + (&column. - mean_value)/num_samples;
				run;
				
				proc sgplot data=_mean_plot_data;
					
					series x=iteration y=mean_value / lineattrs=(color=black);
					
					xaxis label="Number of Iterations";
					yaxis label="Beta Mean";
					
					title "Beta Mean, &variable, &coefficient.";
				run;
			%end;
		%end;
	%end;
%mend beta_trace_plots;

**Produces trace plots of Sigma-u and Sigma-e values and means;
%macro sigma_trace_plots(sigma=,
												 component=,
												 variable_data=,
												 num_mcmc_iterations=,
												 num_burn=,
												 num_thin=,
												 num_variables=,
												 plot_type=);
												 
	%local row col;
												 
	data _value_plot_data;
		set &sigma.;
			
		iteration = _N_;
	run;
		
	%do row = 1 %to &num_variables.;
		
		%do col = 1 %to &row.;
			
			data _NULL_;
				set &variable_data.;
				
				call symputx("row_variable", variable&row., "L");
				call symputx("col_variable", variable&col., "L");
			run;
				
			%if %sysfunc(findw(&plot_type., VALUE, , %str(s))) ^= 0 %then %do;
				
				proc sgplot data=_value_plot_data;
					
					series x=iteration y=sigma_&component._row&row._col&col. / lineattrs=(color=black);
						
					refline &num_burn. / axis=X lineattrs=(color=red);
						
					xaxis label="Number of Iterations";
					yaxis label="Sigma-&component. Value";
						
					title "Sigma-&component. Values, &row_variable., &col_variable.";
				run;
			%end;
				
			%if %sysfunc(findw(&plot_type., MEAN, , %str(s))) ^= 0 %then %do;
			
				data _mean_plot_data;
					set _value_plot_data;
					if iteration > &num_burn. and 
						 iteration <= &num_mcmc_iterations. and
						 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
						 
					retain mean_value num_samples 0;
					
					num_samples = num_samples + 1;
					mean_value = mean_value + (sigma_&component._row&row._col&col. - mean_value)/num_samples;
				run;
					
				proc sgplot data=_mean_plot_data;
					
					series x=iteration y=mean_value / lineattrs=(color=black);
						
					xaxis label="Number of Iterations";
					yaxis label="Sigma-&component. Mean";
						
					title "Sigma-&component. Mean, &row_variable., &col_variable.";
				run;
			%end;
		%end;
	%end;
%mend sigma_trace_plots;

**Produces overlay of Sigma-u and Sigma-e trace plots;
%macro sigma_overlay_plots(sigma_u=,
													 sigma_e=,
													 variable_data=,
													 num_mcmc_iterations=,
													 num_burn=,
													 num_thin=,
													 num_variables=,
													 plot_type=);
													 
	%local row col;
													 
	**combine Sigma-u and Sigma-e data side-by-side;
	data _value_plot_data;
		set &sigma_u.;
		set &sigma_e.;
		
		iteration = _N_;
	run;
	
	%do row = 1 %to &num_variables.;
		
		%do col = 1 %to &row.;
		
			data _NULL_;
				set &variable_data.;
				
				call symputx("row_variable", variable&row., "L");
				call symputx("col_variable", variable&col., "L");
			run;
			
			%if %sysfunc(findw(&plot_type., VALUE, , %str(s))) ^= 0 %then %do;
			
				proc sgplot data=_value_plot_data;
					
					series x=iteration y=sigma_u_row&row._col&col. / lineattrs=(color=blue);
					series x=iteration y=sigma_e_row&row._col&col. / lineattrs=(color=black);
					
					refline &num_burn. / axis=X lineattrs=(color=red);
					
					xaxis label="Number of Iterations";
					yaxis label="Parameter Value";
					
					label sigma_u_row&row._col&col. = "Sigma-u"
								sigma_e_row&row._col&col. = "Sigma-e";
					
					title "Overlay of Sigma-u and Sigma-e Values, &row_variable., &col_variable.";
				run;
				title;
			%end;
			
			%if %sysfunc(findw(&plot_type., MEAN, , %str(s))) ^= 0 %then %do;
			
				data _mean_plot_data;
					set _value_plot_data;
					if iteration > &num_burn. and 
						 iteration <= &num_mcmc_iterations. and 
						 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
						 
					retain sigma_u_mean sigma_e_mean num_samples 0;
					
					num_samples = num_samples + 1;
					
					sigma_u_mean = sigma_u_mean + (sigma_u_row&row._col&col. - sigma_u_mean)/num_samples;
					sigma_e_mean = sigma_e_mean + (sigma_e_row&row._col&col. - sigma_e_mean)/num_samples;
				run;
				
				proc sgplot data=_mean_plot_data;
				
					series x=iteration y=sigma_u_mean / lineattrs=(color=blue);
					series x=iteration y=sigma_e_mean / lineattrs=(color=black);
					
					xaxis label="Number of Iterations";
					yaxis label="Parameter Mean";
					
					label sigma_u_mean = "Sigma-u"
								sigma_e_mean = "Sigma-e";
					
					title "Overlay of Sigma-u and Sigma-e Means, &row_variable., &col_variable.";
				run;
				title;
			%end;
		%end;
	%end;
%mend sigma_overlay_plots;

**Produces trace plots of alpha1 values and means;
%macro alpha1_trace_plots(alpha1=,
													variable_data=,
													never_consumer_covariate_data=,
													num_mcmc_iterations=,
													num_burn=,
													num_thin=,
													num_never_consumer_covariates=,
													num_never_consumer_coefficients=,
													plot_type=);
													
	%local coef_h;
	
	data _NULL_;
		set &variable_data.;
		
		call symputx("variable", variable1);
	run;
	
	data _never_consumer_coefficients (keep = column1-column&num_never_consumer_coefficients. coefficient1-coefficient&num_never_consumer_coefficients.);
		set &never_consumer_covariate_data.;
		
		length column1-column&num_never_consumer_coefficients. $32.;
		length coefficient1-coefficient&num_never_consumer_coefficients $32.;
		
		array columns{&num_never_consumer_coefficients.} column1-column&num_never_consumer_coefficients.;
		array coefficients{&num_never_consumer_coefficients.} coefficient1-coefficient&num_never_consumer_coefficients.;
		
		%if &num_never_consumer_covariates. > 0 %then %do;
		
			array covariates{&num_never_consumer_covariates.} covariate1-covariate&num_never_consumer_covariates.;
		%end;
		
		if intercept = 1 then do;
		
			columns{1} = "alpha1_intercept";
			coefficients{1} = "Intercept";
		end;
		
		do covar_h = 1 to num_covariates;
		
			columns{intercept + covar_h} = cats("alpha1_covariate", covar_h);
			coefficients{intercept + covar_h} = covariates{covar_h};
		end;
	run;
													
	data _value_plot_data;
		set &alpha1.;
		
		iteration = _N_;
	run;
	
	%do coef_h = 1 %to &num_never_consumer_coefficients.;
	
		data _NULL_;
			set _never_consumer_coefficients;
			
			call symputx("column", column&coef_h., "L");
			call symputx("coefficient", coefficient&coef_h., "L");
		run;
	
		%if %sysfunc(findw(&plot_type., VALUE, , %str(s))) ^= 0 %then %do;
			
			proc sgplot data=_value_plot_data;
				
				series x=iteration y=&column. / lineattrs=(color=black);
					
				refline &num_burn. / axis=X lineattrs=(color=red);
					
				xaxis label="Number of Iterations";
				yaxis label="Alpha Value";
					
				title "Alpha Values, &variable., &coefficient.";
			run;
			title;
		%end;
			
		%if %sysfunc(findw(&plot_type., MEAN, , %str(s))) ^= 0 %then %do;
			
				data _mean_plot_data;
					set _value_plot_data;
					if iteration > &num_burn. and 
						 iteration <= &num_mcmc_iterations. and 
						 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
						 
					retain mean_value num_samples 0;
					
					num_samples = num_samples + 1;
					mean_value = mean_value + (&column. - mean_value)/num_samples;
				run;
				
				proc sgplot data=_mean_plot_data;
				
					series x=iteration y=mean_value / lineattrs=(color=black);
					
					xaxis label="Number of Iterations";
					yaxis label="Alpha Mean";
					
					title "Alpha Means, &variable., &coefficient.";
				run;
				title;
		%end;
	%end;
%mend alpha1_trace_plots;

**Produces trace plot of average consumer probabilities for the first episodic food;
%macro consumer_probability_trace_plots(consumer_probabilities=,
																				variable_data=,
																				num_mcmc_iterations=,
																				num_burn=,
																				num_thin=,
																				plot_type=);
																				
	data _NULL_;
		set &variable_data.;
		
		call symputx("variable", variable1);
	run;
																				
	data _value_plot_data;
		set &consumer_probabilities.;
		
		iteration = _N_;
	run;
	
	%if %sysfunc(findw(&plot_type., VALUE, , %str(s))) ^= 0 %then %do;
	
		proc sgplot data=_value_plot_data;
		
			series x=iteration y=conprob1 / lineattrs=(color=black);
			
			refline &num_burn. / axis=X lineattrs=(color=red);
			
			xaxis label="Number of Iterations";
			yaxis label="Consumer Probability Value";
			
			title "Consumer Probability Values, &variable.";
		run;
		title;
	%end;
	
	%if %sysfunc(findw(&plot_type., MEAN, , %str(s))) ^= 0 %then %do;
	
		data _mean_plot_data;
			set _value_plot_data;
			if iteration > &num_burn. and 
				 iteration <= &num_mcmc_iterations. and 
				 mod(iteration, &num_thin.) = mod(&num_burn. + 1, &num_thin.);
				 
			retain mean_value num_samples 0;
			
			num_samples = num_samples + 1;
			mean_value = mean_value + (conprob1 - mean_value)/num_samples;
		run;
		
		proc sgplot data=_mean_plot_data;
			
			series x=iteration y=mean_value;
			
			xaxis label="Number of Iterations";
			yaxis label="Consumer Probability Mean";
			
			title "Consumer Probability Means, &variable.";
		run;
		title;
	%end;
%mend consumer_probability_trace_plots;