/*
Title: Trace Plots for NCI Multivar MCMC
Author: NCI/Information Management Services
Date: 6/9/2025
*/

/*
This example demonstrates the use of trace plots and statistical testing to assess MCMC models.

Workflow:

  1. Data Cleaning
  2. Pre-Processing
  3. MCMC Model
  4. Trace Plots
  5. Gelman-Rubin Test
*/

libname indata "../data";

%include "../macros/ncimultivar.sas";

/*
When performing an analysis with the NCI method, it is important to assess the measurement error model to ensure that it properly represents the distribution of the true parameters. 
The NCI method uses a Markov Chain Monte Carlo (MCMC) method to fit the measurement error model, so it worthwhile to become familiar with common methods used to assess MCMC models. 
These methods can be used to diagnose issues in the model and can help guide refinements to it.

A good model as well as models with some common problems will be shown to build familiarity with diagnosing errors.

The dataset used for this analysis is derived from the 2005-2010 NHANES data.
A subset of six strata (SDMVSTRA) will be used to reduce computation time and allow this example to run in real time.

Intakes of total grain (G_TOTAL) and refined grain (G_REFINED) will be modeled to demonstrate trace plots for different scenarios.

The covariates of interest are smoking status (SMK_REC), age (RIDAGEYR), and sex (RIAGENDR).
The nuisance covariates will be whether the recall was on a weekend (Weekend) or whether the recall was on day 2 (Day2).

The base weight variable (WTDRD1) will be used for all models.
 
For more information on the specifics of the data processing, please see the food and nutrient analysis vignette (food_and_nutrient_analysis.sas). 

Subjects with missing values are removed, and variables are defined for the analysis. 
*/

**subset data;
data input_dataset;
	set indata.nhcvd;
	if SDMVSTRA in (48, 54, 60, 66, 72, 78);
	
	**Define indicator for Day 2;
	Day2 = (DAY = 2);
run;

data input_dataset;
	set input_dataset;
	
	**remove subjects that are missing any covariates or variables;
	if not missing(SMK_REC) and
		 not missing(RIDAGEYR) and
		 not missing(RIAGENDR) and
		 not missing(Weekend) and
		 not missing(Day2) and
		 not missing(G_TOTAL) and
		 not missing(G_REFINED);
	
	**rename variables for readability;
	Total_Grain = G_TOTAL;
	Refined_Grain = G_REFINED;
run;

/*
The variables will now be transformed and standardized for use in the MCMC models.
*/

**Winsorize extreme values;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Total_Grain,
							 weight=WTDRD1,
							 do_winsorization=Y,
							 id=SEQN,
							 repeat_obs=DAY);
							 
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Refined_Grain,
							 weight=WTDRD1,
							 do_winsorization=Y,
							 id=SEQN,
							 repeat_obs=DAY);
							 
data input_dataset;
	set input_dataset;
	
	if Total_Grain > 51.05915 then Total_Grain = 51.05915;
run;

**Find best Box-Cox lambdas for each variable;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Total_Grain,
							 weight=WTDRD1,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3);
							 
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Refined_Grain,
							 weight=WTDRD1,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3);
							 
data boxcox_lambda_data;
	set bc_Total_Grain
			bc_Refined_Grain;
run;

**calculate minimum amounts;
%calculate_minimum_amount(input_data=input_dataset,
													row_subset=%quote(Day2 = 0),
													daily_variables=Total_Grain Refined_Grain);
													
**run pre-processor to get MCMC input data;
%nci_multivar_preprocessor(input_data=input_dataset,
													 daily_variables=Total_Grain Refined_Grain,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=boxcox_lambda_data,
													 minimum_amount_data=minimum_amount_data,
													 outname=model);
													 
/*
The first example will be of a model that has mixed well. 

The %trace_plots() utility is used to generate trace plots for an MCMC model.

An ideal trace plot looks like random movement around an imaginary horizontal line.
The y-value of the line represents the mean of the distribution and the amplitude of the movement represents the variance. 

Each trace plot has a vertical red line at the number of burn-in iterations. 
The area before the line represents the burn-in period where the parameter values are expected to be unstable. 

The values of the mean and variance do not matter for assessing the plot. 
As long as the mean and variance appear stable over time, the model is said to mix well.

Beta parameters tend to stabilize more easily than the variance components Sigma-u and Sigma-e.
The convergence of all parameters should be investigated when assessing a model.
*/

%nci_multivar_mcmc(pre_mcmc_data=model,
									 id=SEQN,
									 repeat_obs=DAY,
									 weight=WTDRD1,
									 daily_variables=Total_Grain,
									 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend,
									 cat_covariates=SMK_REC,
									 cat_reference=3,
									 num_mcmc_iterations=3000,
									 num_burn=1000,
									 num_thin=2,
									 mcmc_seed=9999,
									 outname=good);

ods listing close;

ods pdf file="trace_plots_good_model.pdf";
								 
	%trace_plots(multivar_mcmc_model=good);
ods pdf close;

ods listing;

/*
The following example demonstrates a model with too few burn-in iterations.

Some of the parameter values early on in the chain will not be representative of the parameter distributions.
The inclusion of these values causes the estimates of the posterior means to become less stable.

If a plot does not appear to have settled around a value after the burn-in period, it may indicate that more burn-in iterations are needed.
Multivariate models and episodic variables generally require more burn-in than univariate models and continous variables.
The variance components Sigma-u and Sigma-e often have longer burn-in periods than the beta parameters.
*/

%nci_multivar_mcmc(pre_mcmc_data=model,
									 id=SEQN,
									 repeat_obs=DAY,
									 weight=WTDRD1,
									 daily_variables=Total_Grain,
									 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend,
									 cat_covariates=SMK_REC,
									 cat_reference=3,
									 num_mcmc_iterations=3000,
									 num_burn=50,
									 num_thin=2,
									 mcmc_seed=9999,
									 outname=low_burn);

ods listing close;

ods pdf file="trace_plots_low_burn_in.pdf";		
						 
	%trace_plots(multivar_mcmc_model=low_burn);
ods pdf close;

ods listing;

/*
Another potential issue is an insufficient number of subjects with more than one recall.

The effect of the low sample size can be seen in all of the trace plots, but it is especially apparent in the variance components (Sigma-e and Sigma-u). 
While there isn't a clear trend up or down, the center of the trace appears choppy which shows that the mean is not able to settle on a value even well after the burn-in period.

Note that the number of subjects with one recall does not have any bearing on model convergence. 
The limiting factor is the number of subjects that have two or more observations. 
As a general rule, at least 50 subjects should have two or more observations were consumption is non-zero for proper mixing. 
More complicated relationships in the data (such as high correlation between different variance components) will require more subjects to fit. 
*/

**subset input dataset so that only the first 10 subjects have a second recall;
data first_recall
		 second_recall;
	set input_dataset;
	
	if Day2 = 0 then output first_recall;
	else if Day2 = 1 then output second_recall;
run;

data second_recall;
	set second_recall;
	if _N_ <= 10;
run;

data input_subset;
	set first_recall
			second_recall;
run;

%nci_multivar_preprocessor(input_data=input_subset,
													 daily_variables=Total_Grain,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=boxcox_lambda_data,
													 minimum_amount_data=minimum_amount_data,
													 outname=subset);
													 
%nci_multivar_mcmc(pre_mcmc_data=subset,
									 id=SEQN,
									 repeat_obs=DAY,
									 weight=WTDRD1,
									 daily_variables=Total_Grain,
									 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day Weekend,
									 cat_covariates=SMK_REC,
									 cat_reference=3,
									 num_mcmc_iterations=3000,
									 num_burn=1000,
									 num_thin=2,
									 mcmc_seed=9999,
									 outname=subset);
									 
ods listing close;

ods pdf file="trace_plots_small_sample.pdf";			
						 
	%trace_plots(multivar_mcmc_model=subset);
ods pdf close;

ods listing;

/*
Variables that are highly correlated or subsets of each other will affect the model fit significantly. 
A common cause of this is not properly disaggregating variables during pre-processing.

The trace plots have a long burn-in and have trouble settling around a value.
This can sometimes resemble the plots from having a low sample size.
If the sample size is sufficient, then this pattern could indicate an issue with disaggregation.
*/

%nci_multivar_mcmc(pre_mcmc_data=model,
									 id=SEQN,
									 repeat_obs=DAY,
									 weight=WTDRD1,
									 daily_variables=Total_Grain Refined_Grain,
									 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend,
									 cat_covariates=SMK_REC,
									 cat_reference=3,
									 num_mcmc_iterations=3000,
									 num_burn=1000,
									 num_thin=2,
									 mcmc_seed=9999,
									 outname=corr);
									 
ods listing close;

ods pdf file="trace_plots_correlated_model.pdf";	
								 
	%trace_plots(multivar_mcmc_model=corr);
ods pdf close;

ods listing;

/*
The Gelman-Rubin test is an alternative method assessing convergence that gives numerical values instead of graphs.
The test works by creating multiple MCMC chains with different random seeds and calculating the within-chain and between-chain variation (Gelman and Rubin, 1992). 
If the model parameters converge, there should be little to no difference between different chains which will cause the between-chain variance to fall to zero. 

The Gelman-Rubin statistic is the square root of the ratio between the total variance and the within-chain variance of a parameter. 
Convergence is indicated by the Gelman-Rubin statistics for all of the parameters being close to 1. 
There is no exact rule for how close to 1 the statistics must be, though an upper bound of 1.1 for convergence has been suggested by Gelman, et al. (2004).
In general, iteration parameters should be chosen that get the Gelman-Rubin statistics as close to 1 as possible while balancing computation time.

The %gelman_rubin() function computes Gelman-Rubin statistics for each MCMC parameter.
The Gelman-Rubin test can be used as evidence that the good model shown earlier in the vignette has converged.
*/

%gelman_rubin(num_chains=5,
							pre_mcmc_data=model,
							id=SEQN,
							repeat_obs=DAY,
							weight=WTDRD1,
							daily_variables=Total_Grain,
							default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend,
							cat_covariates=SMK_REC,
							cat_reference=3,
							num_mcmc_iterations=3000,
							num_burn=1000,
							num_thin=2,
							initial_mcmc_seed=9999,
							outname=model_gelman_rubin);
							
proc print data=model_gelman_rubin; 

	title "Gelman-Rubin Statistics";
run;

/*
The advantage to the Gelman-Rubin test is that it gives a numerical result as opposed to a subjective assessment of a graph. 
On the other hand, trace plots are much faster to create and can be used on a model that has already been fit.
Additionally, the specific patterns in trace plots can give clues about why the model has not converged.
The best practice is to use both methods for a more thorough assessment.
*/