/* History:
2026-02-04 Update to allow for two-way interactions.
2025-12-17 Update macro to allow for categorical covariates.
           The options cat_covariates, cat_reference, and cat_maxn have been added.
2025-11-21 Load the v143_modules catalog.		   
2025-10-30 Add options DEBUG and v143code to nci_multivar_mcmc macro.
*/

%include "&sysincludefiledir./internal/mcmc_check_input.sas";
%include "&sysincludefiledir./internal/mcmc_process_variables.sas";

%include "&sysincludefiledir./internal/iml_modules/mcmc_check_iml.sas";
%include "&sysincludefiledir./internal/iml_modules/mcmc_initialize.sas";
%include "&sysincludefiledir./internal/iml_modules/mcmc_likelihood.sas";
%include "&sysincludefiledir./internal/iml_modules/mcmc_main_loop.sas";
%include "&sysincludefiledir./internal/iml_modules/mcmc_update_parameters.sas";
%include "&sysincludefiledir./internal/iml_modules/v143_modules.sas";

/*
Fit MCMC measurement error model

Description:

	Fits a measurement error model for one or more episodic and/or
	daily variables using an MCMC method.
	
Parameters:

	- pre_mcmc_lib: The SAS library that contains the pre-MCMC datasets. (default = WORK)
	- pre_mcmc_data: Prefix for the pre-MCMC datasets (max 23 characters). The following datasets with the prefix must be present:
									 	 - _mcmc_in: A SAS data set with indicators, standardized amounts, and covariates for the MCMC model.
									 	 - _backtran: A SAS data set with the following columns:
									 	 							  - variable: The name of the variable.
									 	 							  - tran_lambda: The value of the lambda used to transform the variable.
									 	 							  - minamount: Half of the minimum non-zero value of the variable.
									 	 							  - tran_center: The mean of the Box-Cox transformed variable before standardization.
									 	 							  - tran_scale: The standard deviation of the Box-Cox transformed variable before standardization divided by sqrt(2).
									 	 							  - biomarker: Binary flag indicating whether the variable is a biomarker assumed to be unbiased on the transformed scale.
	- id: Variable that identifies each subject.
	- repeat_obs: Variable that distinguishes repeat observations for each subject.
	- weight: Variable with weighting for each subject.
	- episodic_variables: Space-delimited list of episodic variables.
	- episodic_indicators: Space-delimited list of consumption indicator variables for episodic variables. 
												 Not needed if pre_mcmc_data was created by nci_multivar_preprocessor.
	- episodic_amounts: Space-delimited list of consumption amount variables for episodic variables.
											Not needed if pre_mcmc_data was created by nci_multivar_preprocessor.
	- daily_variables: Space-delimited list of daily variables.
	- daily_amounts: Space-delimited list of consumption amount variables for daily variables.
									 Not needed if pre_mcmc_data was created by nci_multivar_preprocessor.
	- never_consumer_variable: One episodic variable to allow never-consumers. 
														 Can be a variable already listed as episodic or a new variable.
	- never_consumer_indicator: Consumption indicator variable for the never-consumer variable.
															Not needed if pre_mcmc_data was created by nci_multivar_preprocessor.
	- never_consumer_amount: Consumption amount variable for the never-consumer variable.
													 Not needed if pre_mcmc_data was created by nci_multivar_preprocessor.
	- default_covariates: Space-delimited list of covariates to be used for episodic indicators, episodic amounts, and daily amounts. 
												Does not affect the never-consumer model if present.
	- default_intercept: Flag to include an intercept in the models for episodic indicators, episodic amounts, and daily amounts. (Y/N, default = Y)
											 Does not affect the never-consumer model if present.
	- episodic_indicator_covariates: Space-delimited list of covariates to be used for episodic indicators. 
																	 If specified, overwrites default_covariates for episodic indicators.
	- episodic_indicator_intercept: Flag to include an intercept in the model for episodic indicators. (Y/N).
																	If specified, overwrites default_intercept for episodic indicators.
	- episodic_amount_covariates: Space-delimited list of covariates to be used for episodic amounts. 
																If specified, overwrites default_covariates for episodic amounts.
	- episodic_amount_intercept: Flag to include an intercept in the model for episodic amounts. (Y/N)
															 If specified, overwrites default_intercept for episodic amounts.
	- daily_amount_covariates: Space-delimited list of covariates to be used for daily amounts. 
														 If specified, overwrites default_covariates for daily amounts.
	- daily_amount_intercept: Flag to include an intercept in the model for daily amounts. (Y/N)
														If specified, overwrites default_intercept for daily amounts.
	- individual_covariates: A SAS data set of covariate lists for individual indicators and amounts.
													 It must contain the following columns:
													 	 - variable: Name of the indicator or amount.
													 	 - covariates: A space-delimited string of covariates for the indicator or amount.
													 Overrides default_covariates, episodic_indicator_covariates, episodic_amount_covariates, and daily_amount_covariates for each named variable.
	- individual_intercept: A SAS data set of intercepts for individual indicators and amounts.
													It must contain the following columns:
														- variable: Name of the indicator or amount.
														- intercept: Flag to include an intercept in the model for the indicator or amount (Y/N).
													Overrides default_intercept, episodic_indicator_intercept, episodic_amount_intercept, and daily_amount_intercept for each named variable.
	- never_consumer_covariates: Space-delimited list of covariates to be used for the never-consumer model. 
															 Not affected by default_covariates and must be specified separately.
	- never_consumer_intercept: Flag to include intercept in the never-consumer model.
															Not affected by default_intercept and must be specified separately. (Y/N, default = Y)
	- mcmc_seed: Positive integer starting seed for the random number generator. 
							 If blank, uses a randomly generated integer from 1 to 2*10^7, inclusive. (default = blank)
	- num_mcmc_iterations: Integer specifying the total number (including burn-in) of iterations in the MCMC chain. (default = 12000)
	- num_burn: Integer specifying the number of burn-in iterations in the MCMC chain. Must be smaller than num_mcmc_iterations. (default = 2000)
	- num_thin: Integer specifying the number of iterations between MCMC samples used for calculating posterior means. (default = 25)
	- num_post: Integer specifying the number of random effect (u) matrix draws done after the MCMC chain is finished, conditional on the posterior mean of the parameters. 
							Used for regression calibration, see Regression Calibration section. (default = 0)
	- sigma_u_prior: A SAS data set specifying a prior covariance matrix for the random effects (u). 
									 Must be a valid covariance matrix with dimensions equal to the total number of episodic indicators, episodic amounts, and daily amounts.
	- sigma_u_constant: Flag specifying if the covariance matrix of the random effects (u) should remain the same for each iteration. (Y/N, default = N)
	- save_mcmc_u: Flag specifying if random effects (u) matrices should be saved from the main MCMC chain.
								 If Y, saves u matrix every num_thin iterations after num_burn. 
								 This is used for fitting models with techniques such as multiple imputation. 
								 Does not affect conditional u matrices drawn after the MCMC chain. (Y/N, default = N)
	- save_all_u: Flag specifying if random effects (u) matrices are saved for every iteration. 
								Greatly increases disk space needed to store output. 
								Has no effect if save_u_main is N. (Y/N, default = N)
	- do_log_likelhood: Flag specifying if an estimate of the log-likelihood should be calculated after the MCMC chain has finished.
											See the Log-likelihood section for details. (Y/N, default = N)
	- low_memory_mode: Y/N indicator. If Y, MCMC parameter traces and U matrices are written to the output datasets one iteration at a time to save memory. (default = N)
	- outlib: The SAS library to store output datasets. (default = WORK)
	- outname: The prefix used to name output datasets. (max 18 characters)
	- DEBUG: Flag to print debugging information. (0/1, default = 0)
	- v143code: Flag to run alternative code for IML version 14.3. (0/1, default is determined at run time)
	- 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.
	
Output:

	The following SAS datasets are created in outlib and prefixed with outname:
	
		- _beta1_mean-_betaJ_mean: Posterior means (point estimates) of fixed effect coefficients for each of the J indicators and amounts.
		- _beta1_var-_betaJ_var: Posterior variances of the fixed effect coefficients for each of the J indicators and amounts.
		- _beta1-_betaJ: Fixed effect coefficients at every MCMC iteration for each of the J indicators and amounts.
		- _sigma_u_mean: Posterior means (point estimates) of between-person covariance matrix.
		- _sigma_u_var: Posterior variances of between-person covariance matrix.
		- _sigma_u: Upper triangle of between-person covariance matrix at every MCMC iteration.
		- _sigma_e_mean: Posterior means (point estimates) of within-person covariance matrix.
		- _sigma_e_var: Posterior variances of within-person covariance matrix.
		- _sigma_e: Upper triangle of within-person covariance matrix at every MCMC iteration.
		- _u_mcmc: Random effect (u) matrices saved from the main MCMC chain.
							 Contains the following columns:
							 	 - iteration: The MCMC iteration that the u matrix was saved from.
							 	 - u_col1-u_colJ: The random intercepts for each of the J indicators and amounts.
		- _u_post: Random effect (u) matrices drawn after the MCMC chain conditional on the posterior means of the parameters.
							 Contains the following columns:
							 	 - post_mcmc_iteration: The number of iterations after the MCMC chain that the u matrix was drawn.
							 	 - u_col1-u_colJ: The random intercepts for each of the J indicators and amounts.
		- _subjects: The subject information for the MCMC model.
								 Contains the following columns:
								 	 - subject_id: The subject identifier for each unique subject used in the model.
								 	 - subject_weighting: The weighting for each unique subject used in the model.
		- _iters: The iteration parameters used in the MCMC.
							Contains the following columns:
						  	- num_traces: Number of MCMC iterations saved in the output. 
						  	              Used to keep track of iterations saved if the iteration numbers are changed with %burn_and_thin().
								- num_mcmc_iterations: Total number of iterations in the MCMC chain including burn-in, same as num_mcmc_iterations parameter.
								- num_burn: Number of burn-in iterations that will be discarded when calculating posterior means, same as num_burn parameter.
								- num_thin: Number of iterations between MCMC samples used for calculating posterior means, same as num_thin parameter.
								- num_post: Number of u matrices drawn after the MCMC chain conditional on the posterior means of the parameters, same as num_post parameter.
								- has_saved_mcmc_u: 1/0 flag for whether u matrices were saved from the main MCMC chain, same as the save_u_main parameter.
		- _vars: Variables, indicators, and amount used in the model.
						 Contains the following columns:
						 	 - episodic_variable1-episodic_variableM: The names of the M episodic variables used in the model.
							 - episodic_indicator1-episodic_indicatorM: The names of the indicators for the M episodic variables used in the model.
							 - episodic_amount1-episodic_amountM: The names of the amounts for the M episodic variables used in the model.
							 - daily_variable1-daily_variableN: The names of the N episodic variables used in the model.
							 - daily_amount1-daily_amountN: The names of the amounts for the N daily variables used in the model.
							 - num_episodic: The number of episodic variables in the model.
							 - num_daily: The number of daily variables in the model.
							 - num_variables: Total number of indicators and amounts in the model.
							 - has_never_consumers: 1/0 flag for whether never-consumers are allowed for episodic_variable1.
		- _covars: Covariates for each indicator and amount in the model.
							 Contains the following columns:
							 	 - variable: The name of the indicator or amount.
								 - covariate1-covariateP: The names of the P covariates used in the model for the variable.
								 - intercept: 1/0 flag for whether an intercept was used in the model for the variable.
								 - num_covariates: The number of covariates specified for each variable.
								 - num_coefficients: Number of beta coefficients for each variable (number of covariates plus the intercept).
		- _backtran: The backtransformation dataset passed through from pre_mcmc_data.
		- _seed: The random number generator seed used to generate the results, see the mcmc_seed parameter for details.
		
	If a never-consumer variable is modeled, the following SAS data sets are created in outlib and prefixed with outname:
	
		- _alpha1_mean: Posterior means (point estimates) of never-consumer model coefficients.
		- _alpha1_var: Posterior variances of never-consumer model coefficients.
		- _alpha1: Never-consumer model coefficients at every MCMC iteration.
		- _conprob1_mean: Posterior mean (point estimate) of the probability of a person being a consumer of the never-consumer variable.
		- _conprob1_var: Posterior variance of the probability of a person being a consumer of the never-consumer variable.
		- _conprob1: Probability of a person being a consumer of the never-consumer variable at every MCMC iteration.
		- _covarsnc: Covariates for the never-consumer model.
								 Contains the following columns:
								 	 - covariate1-covariateP: The names of the P covariates used in the never-consumer model.
								 	 - intercept: 1/0 flag for whether an intercept was included in the never-consumer model.		
									 - num_covariates: The number of covariates in the never-consumer model.
									 - num_coefficients: Number of alpha coefficients in the never-consumer model (number of covariates plus the intercept).
												
	If do_log_likelihood is Y, the following SAS data sets are created in outlib and prefixed outname:
	
		- _ll: The estimated marginal log-likelihood of the model.
		
	A listing of the output datasets is created in outlib and prefixed with outname:
				 										
		- _obj: A SAS data set with the following columns:
							- name: Name of the output component.
							- dataset: Name of the SAS dataset containing the output component.
												
About the Model:

	The model fit by this function is a multivariate
	mixed-effects model with fixed effects (beta) and a random intercept for
	each variable (u). Each daily variable has a single continuous outcome
	variable for the amount consumed. Each episodic variable has two model
	parts, a binary indicator of consumption and a continuous variable for
	amount consumed. A probit link function is used for modeling the binary
	indicator. The parts of the model (daily amounts, episodic indicators, and
	episodic amounts) can have separate covariate lists. The random intercepts
	are assumed to have multivariate Gaussian distribution with covariance
	matrix Sigma-u which represents the between-person variance component. 
	The residual error is also assumed to have a multivariate Gaussian 
	distribution with covariance matrix Sigma-e which represents the 
	within-person variance component.
	
MCMC Procedure:

	The model is fit using a Gibbs sampler MCMC
	procedure. At each iteration, each parameter is updated sequentially using
	its distribution conditional on the other parameters. The conditional
	distributions and priors are outlined in Zhang, et al. (2011). After a
	specified burn-in period, subsequent samples of each parameter are averaged
	to calculate a posterior mean. The posterior means are the point estimates
	for the parameters. To ensure that the samples used to calculated the
	posterior mean are independent, the samples can be 'thinned'. For example,
	a thinning number of 5 means that every 5th sample after the burn-in period
	is used to calculate the posterior mean. When doing regression calibration,
	additional draws of the random intercept (u) matrix can be made
	conditional on the posterior means of beta, Sigma-u, and Sigma-e. This is
	useful when performing measurement error correction with regression
	calibration.
	
Never-consumers:

	The model can allow some subjects to be
	never-consumers of one episodic variable. If an episodic variable with
	never-consumers is specified, then the probability that each subject has
	any consumption is included as an additional part of the model for that
	variable. As with the other model parts, the never-consumer part of the
	model has its own covariate matrix and coefficients. To distiguish the
	coefficients from the regular fixed effect coefficients, they are referred
	to as alpha1. Only one variable can be allowed to have never-consumers.
	
Log-likelihood:

	An estimate of the log-likelihood marginalized over
	the random effects (u) can be output after the MCMC iterations have
	finished. The random effects are integrated out using the Laplace
	approximation. The posterior mode of the random effects is found by
	the BFGS algorithm with derivatives approximated by the central difference
	method. Since the model is fit using an MCMC method, the marginal
	log-likelihood is never used in fitting the model. The log-likelihood
	output by this function corresponds well to the log-likelihood found by
	maximum likelihood for univariate and bivariate models. For
	higher-dimensional models and models allowing never-consumers that cannot
	be fit using maximum likelihood, the estimated likelihood is experimental
	and other measures of model fit should be considered. Estimating the
	log-likelihood increases computation time, especially as the number of
	variables increases.
	
Regression Calibration:

	The model fit by nci_multivar_mcmc can
	be used in regression calibration for measurement error correction. When
	doing regression calibration, random effect (u) matrices must be drawn
	after the MCMC chain conditional on the posterior means of the parameters.
	This can be done by setting num_post to the number of post-MCMC u draws
	to be made.
*/
%macro nci_multivar_mcmc(pre_mcmc_lib=WORK,
												 pre_mcmc_data=,
												 id=,
												 repeat_obs=,
												 weight=,
												 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,
												 mcmc_seed=,
												 num_mcmc_iterations=12000,
												 num_burn=2000,
												 num_thin=25,
												 num_post=0,
												 sigma_u_prior=,
												 sigma_u_constant=N,
												 save_mcmc_u=N,
												 save_all_u=N,
												 do_log_likelihood=N,
												 low_memory_mode=N,
												 outlib=WORK,
												 outname=,
												 DEBUG=0,
												 v143code=,
												 cat_covariates=,
												 cat_reference=,
												 cat_maxn=25);
	
	%local 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;	
		
	** Create dummy variables for categorical covariates;
	%let dat      = &pre_mcmc_lib..&pre_mcmc_data._mcmc_in; /* input data */
	%let catvdata = _catv_info_tmp; /* output data of categorical variables */
	%let allv     = &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 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 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);
	%updateDataWithInteractions(&individual_covariates, &intvdata, var="covariates");
	
	**1. Check input parameters and process variable and covariate lists;
	%check_and_process_input(pre_mcmc_data=&pre_mcmc_lib..&pre_mcmc_data.,
													 id=&id.,
													 weight=&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=&num_post.,
													 save_mcmc_u=&save_mcmc_u.,
													 save_all_u=&save_all_u.,
													 outname=&outlib..&outname.);
													 
	data &outlib..&outname._backtran;
		set &pre_mcmc_lib..&pre_mcmc_data._backtran;
	run;
													 
	**2. Sort input dataset;
	proc sort data=&pre_mcmc_lib..&pre_mcmc_data._mcmc_in out=_mcmc_input; by &id. &repeat_obs.; run;
	
	proc iml;
	
		reset storage=WORK.v143_modules;
		load module=_all_;
		
		reset storage=WORK.mcmc_modules;
		load module=_all_;
		
		**3. Load number of variables and covariates into IML;
		run load_mcmc_variables(episodic_indicators,
														episodic_amounts,
														daily_amounts,
														covariates,
														intercepts,
														num_episodic,
														num_daily,
														num_variables,
														num_covariates,
														num_coefficients,
														has_never_consumers,
														never_consumer_covariates,
														never_consumer_intercept,
														num_never_consumer_covariates,
														num_never_consumer_coefficients,
														"&outlib..&outname._vars",
														"&outlib..&outname._covars",
														"&outlib..&outname._covarsnc");
	
		**4. Initialize parameters for main loop: set seed, variable and covariate matrices, and priors;
		mcmc_parameters = initialize_mcmc("_mcmc_input",
																			"&id.",
																			"&repeat_obs.",
																			"&weight.",
																			episodic_indicators,
																			episodic_amounts,
																			daily_amounts,
																			covariates,
																			intercepts,
																			"&sigma_u_prior.",
																			num_episodic,
																			num_daily,
																			num_variables,
																			num_covariates,
																			num_coefficients,
																			has_never_consumers,
																			never_consumer_covariates,
																			never_consumer_intercept,
																			num_never_consumer_covariates,
																			num_never_consumer_coefficients,
																			&mcmc_seed.);
		
		**4. Main MCMC loop;
		run mcmc_main_loop("&outlib..&outname.",
											 mcmc_parameters$"priors"$"beta_prior"$"mean",
											 mcmc_parameters$"priors"$"beta_prior"$"covariance",
											 mcmc_parameters$"priors"$"r_matrix_prior",
											 mcmc_parameters$"priors"$"theta_matrix_prior",
											 mcmc_parameters$"priors"$"v_matrix_prior",
											 mcmc_parameters$"priors"$"sigma_u_prior",
											 mcmc_parameters$"priors"$"sigma_e_prior",
											 mcmc_parameters$"priors"$"u_matrix_prior",
											 mcmc_parameters$"priors"$"w_matrix_prior",
											 mcmc_parameters$"priors"$"xbeta_prior",
											 mcmc_parameters$"priors"$"xbeta_u_prior",
											 mcmc_parameters$"subject_data"$"recall_availability",
											 mcmc_parameters$"subject_data"$"subject_weighting",
											 mcmc_parameters$"variable_matrix_data"$"episodic_indicator_matrices",
											 mcmc_parameters$"covariate_matrix_data"$"covariate_matrices",
											 mcmc_parameters$"covariate_matrix_data"$"wt_covariate_matrices",
											 mcmc_parameters$"covariate_matrix_data"$"wt_covariate_sq_sums",
											 &num_mcmc_iterations.,
											 &num_burn.,
											 &num_thin.,
											 &num_post.,
											 mcmc_parameters$"subject_data"$"num_subjects",
											 mcmc_parameters$"subject_data"$"num_recalls",
											 num_episodic,
											 num_daily,
											 num_variables,
											 num_coefficients,
											 "&sigma_u_constant.",
											 "&save_mcmc_u.",
											 "&save_all_u.",
											 "&do_log_likelihood.",
											 mcmc_parameters$"variable_matrix_data"$"variable_names",
											 mcmc_parameters$"covariate_matrix_data"$"coefficient_names",
											 has_never_consumers,
											 mcmc_parameters$"priors"$"alpha1_prior"$"mean",
											 mcmc_parameters$"priors"$"alpha1_prior"$"covariance",
											 mcmc_parameters$"priors"$"consumer_probabilities_prior",
											 mcmc_parameters$"covariate_matrix_data"$"never_consumer_covariate_matrix",
											 num_never_consumer_coefficients,
											 mcmc_parameters$"covariate_matrix_data"$"never_consumer_coefficient_names",
											 "&low_memory_mode.");
											 
		**5. Save MCMC model information;
		**subjects;
		subject_id = mcmc_parameters$"subject_data"$"subject_id";
		subject_weighting = mcmc_parameters$"subject_data"$"subject_weighting";
		
		create &outlib..&outname._subjects var {subject_id subject_weighting};
			append;
		close &outlib..&outname._subjects;
		
		**iterations;
		num_traces = &num_mcmc_iterations.;
		
		num_mcmc_iterations = &num_mcmc_iterations.;
		num_burn = &num_burn.;
		num_thin = &num_thin.;
		
		num_post = &num_post.;
		
		if upcase("&save_mcmc_u.") = "Y" then has_saved_mcmc_u = 1;
		else has_saved_mcmc_u = 0;
		
		create &outlib..&outname._iters var {num_traces num_mcmc_iterations num_burn num_thin num_post has_saved_mcmc_u};
			append;
		close &outlib..&outname._iters;
		
		**seed;
		seed = mcmc_parameters$"mcmc_seed";
		
		create &outlib..&outname._seed var {seed};
			append;
		close &outlib..&outname._seed;
	quit;
	
	**Create listing of output datasets;
	data _NULL_;
		set &outlib..&outname._vars;
		
		call symputx("num_variables", num_variables, "L");
		
		call symputx("has_never_consumers", has_never_consumers, "L");
	run;
	
	data &outlib..&outname._obj (keep = name dataset);
		
		length name dataset $32.;
		
		%do var_j = 1 %to &num_variables.;
		
			name = "beta&var_j._mean";
			dataset = "&outname._beta&var_j._mean";
			output;
			
			name = "beta&var_j._var";
			dataset = "&outname._beta&var_j._var";
			output;
			
			name = "beta&var_j.";
			dataset = "&outname._beta&var_j.";
			output;
		%end;
		
		name = "sigma_u_mean";
		dataset = "&outname._sigma_u_mean";
		output;
		
		name = "sigma_u_var";
		dataset = "&outname._sigma_u_var";
		output;
		
		name = "sigma_u";
		dataset = "&outname._sigma_u";
		output;
		
		name = "sigma_e_mean";
		dataset = "&outname._sigma_e_mean";
		output;
		
		name = "sigma_e_var";
		dataset = "&outname._sigma_e_var";
		output;
		
		name = "sigma_e";
		dataset = "&outname._sigma_e";
		output;
		
		%if %upcase(&save_mcmc_u.) = Y %then %do;
			
			name = "u_mcmc";
			dataset = "&outname._u_mcmc";
			output;
		%end;
		
		%if &num_post. > 0 %then %do;
			
			name = "u_post";
			dataset = "&outname._u_post";
			output;
		%end;
		
		name = "subjects";
		dataset = "&outname._subjects";
		output;
		
		name = "iters";
		dataset = "&outname._iters";
		output;
		
		name = "vars";
		dataset = "&outname._vars";
		output;
		
		name = "covars";
		dataset = "&outname._covars";
		output;
		
		name = "backtran";
		dataset = "&outname._backtran";
		output;
		
		name = "seed";
		dataset = "&outname._seed";
		output;
		
		%if &has_never_consumers. = 1 %then %do;
		
			name = "alpha1_mean";
			dataset = "&outname._alpha1_mean";
			output;
			
			name = "alpha1_var";
			dataset = "&outname._alpha1_var";
			output;
			
			name = "alpha1";
			dataset = "&outname._alpha1";
			output;
			
			name = "conprob1_mean";
			dataset = "&outname._conprob1_mean";
			output;
			
			name = "conprob1_var";
			dataset = "&outname._conprob1_var";
			output;
			
			name = "conprob1";
			dataset = "&outname._conprob1";
			output;
			
			name = "covarsnc";
			dataset = "&outname._covarsnc";
			output;
		%end;
		
		%if %upcase(&do_log_likelihood.) = Y %then %do;
		
			name = "ll";
			dataset = "&outname._ll";
			output;
		%end;
	run;
%mend nci_multivar_mcmc;