/*
Title: NCI Method Food and Nutrient
Author: NCI/Information Management Services
Date: 6/9/2025
*/

/*
This example demonstrates estimating a multivariate distribution and ratios with standard errors from balanced repeated replication (BRR).

Workflow:

  1. Data Cleaning
  2. Replicate Weights
  3. Pre-Processing
  
  For each replicate:
    4. MCMC Model
    5. Simulate Usual Intakes
    6. Summary Statistics
    
  7. Variances and Confidence Intervals
*/

libname indata "../data";

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

/*
The NCI method can be used to obtain distributions of usual intake ratios.

Ratios should be calculated from the usual intake distributions of the components (ratio of usuals).
This requires fitting a multivariate model in order to account for relationships among the components.

It is NOT recommended to take the ratio of raw intakes and apply the NCI method to the ratio (usual of ratio).

Whole grain density, the ratio of  whole grain (G_WHOLE) to energy (TKCAL), will be calculated from the 2005-2010 NHANES data.
A bivariate model will fit with whole grains as an episodically consumed food and energy as a daily consumed nutrient.
A subset of six strata (SDMVSTRA) will be used to reduce computation time and allow this example to run in real time.

The covariates being examined are smoking status (SMK_REC), age (RIDAGEYR), and sex (RIAGENDR). 
Two nuisance covariates will be adjusted for, whether the recall was on a weekend (Weekend) and and whether the recall is on day 2 (Day2).

The WTDRD1 variable is the base weighting for each observation.

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_WHOLE) and
		 not missing(TKCAL);
	
	**rename whole grain and energy variables for readability;
	Whole_Grain = G_WHOLE;
	Energy = TKCAL;
run;

/*
BRR weights will be added to the dataset.
*/

%let fay_factor = 0.7;

%brr_weights(input_data=input_dataset,
						 id=SEQN,
						 strata=SDMVSTRA,
						 psu=SDMVPSU,
						 cell=PSCELL,
						 weight=WTDRD1,
						 fay_factor=&fay_factor.,
						 outname=input_dataset);

/*					 
Next, extreme observations will be identified for Winsorization. 
The %boxcox_survey() macro must be run for both variables. 
The non-zero values of the first recall will be used. 

Since whole grain is an episodic variable, Winsorization for low outliers is done slightly differently. 
Instead of changing values that are too small to a minimum threshold, they are set to zero and treated as a non-consumption observation. 
This behavior is toggled by setting the is_episodic parameter to Y.
*/

%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Whole_Grain,
							 is_episodic=Y,
							 weight=RepWt_0,
							 do_winsorization=Y,
							 iqr_multiple=2,
							 id=SEQN,
							 repeat_obs=DAY);

%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Energy,
							 weight=RepWt_0,
							 do_winsorization=Y,
							 iqr_multiple=2,
							 id=SEQN,
							 repeat_obs=DAY);
							 
data input_dataset;
	set input_dataset;
	
	**Winsorize whole grain;
	if Whole_Grain > 10.71163 then Whole_Grain = 10.71163;
	
	**Winsorize energy;
	if Energy < 269.0701 then Energy = 269.0701;
	if Energy > 8026.0436 then Energy = 8026.0436;
run;

/*
Next, the best Box-Cox lambda parameter for each variable can be found using the Winsorized data in the presence of covariates. 
A Box-Cox survey must be run for each variable used in the model, and the output datasets are concatenated.
*/

%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Whole_Grain,
							 weight=RepWt_0,
							 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=Energy,
							 weight=RepWt_0,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3);
							 
data boxcox_lambda_data;
	set bc_Whole_Grain
			bc_Energy;
run;

/*
Next, the minimum amounts for whole grain and energy can be calculated using the first recall. 

The %calculate_minimum_amount() function can handle multiple variables at the same time.
Episodically consumed foods are specified with episodic_variables while daily consumed nutrients are specified with daily_variables.
If there are multiple daily and/or episodic variables, they should be specified as space-delimited lists in their respective parameters.
*/

%calculate_minimum_amount(input_data=input_dataset,
													row_subset=%quote(Day2 = 0),
													episodic_variables=Whole_Grain,
													daily_variables=Energy);
													
/*
The %nci_multivar_preprocessor() macro can now generate the MCMC input data. 
Whole grains should be input in the episodic_variables parameter to tell the macro to generate a consumption indicator variable in addition to an amount variable.

The episodic_variables and daily_variables parameters are space-delimited lists of episodically consumed foods and daily consumed nutrients, respectively.

Continuous covariates to be standardized are specified with continuous_covariates just like the daily nutrient vignette (daily_nutrient_analysis.sas).
*/

%nci_multivar_preprocessor(input_data=input_dataset,
													 episodic_variables=Whole_Grain,
													 daily_variables=Energy,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=boxcox_lambda_data,
													 minimum_amount_data=minimum_amount_data,
													 outname=model);

/*
The preprocessor creates the indicator variable ind_Whole_Grain` for and amount variable amt_Whole_Grain for Whole_Grain.
The indicator variable shows whether or not the subject reported consuming whole grains on a particular day. 
The amount is always non-zero when the indicator is 1 and missing when the indicator is 0.

The first 10 observations demonstrate the relationship between indicators and amounts.
*/

proc print data=model_mcmc_in (obs = 10);
	
	var SEQN DAY ind_Whole_Grain amt_Whole_Grain;
	
	title "Indicator and Amount Variables";
run;

/*
MCMC models are fit using both whole grains and energy using the base weight and each set of replicate weights.

Episodic and daily variables are specified in episodic_variables and daily_variables as in the previous steps.

Note that the number of iterations and burn-in is higher than for a single daily consumed nutrient.
This is because models with episodic variables and multivariate models take more iterations to converge.

Since the same population will be used to fit the MCMC model and generate usual intakes, 500 conditional U matrix draws will be made.
*/

%let num_brr = 8;

%macro mcmc_brr(num_brr=);

	%do brr_rep = 0 %to &num_brr.;
	
		%put Starting Iteration &brr_rep.;
		
		%nci_multivar_mcmc(pre_mcmc_data=model,
											 id=SEQN,
											 repeat_obs=DAY,
											 weight=RepWt_&brr_rep.,
											 episodic_variables=Whole_Grain,
											 daily_variables=Energy,
											 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend,
											 cat_covariates=SMK_REC,
											 cat_reference=3,
											 num_mcmc_iterations=4000,
											 num_burn=2000,
											 num_thin=5,
											 num_post=500,
											 mcmc_seed=%eval(9999 + &brr_rep.),
											 outname=model_brr&brr_rep.);
	%end;
%mend mcmc_brr;

%mcmc_brr(num_brr=&num_brr.);

/*
The population base for simulating usual intakes is generated from the input data in the same way as the daily nutrient vignette (daily_nutrient_analysis.sas).
*/

proc sort data=model_mcmc_in; by SEQN; run;

data distrib_pop;
	set model_mcmc_in;
	by SEQN;
	
	**get first instance of each subject;
	if first.SEQN then do;
	
		**Set Day 2 to zero to factor out the effect of Day 2 recalls;
		Day2 = 0;
	
		**create repeats of each subject for weekday and weekend consumption;
		Weekend = 0;
		Weekend_Weight = 4;
		output;
	
		Weekend = 1;
		Weekend_Weight = 3;
		output;
	end;
run;

/*
Usual intakes are simulated with %nci_multivar_distrib() from the population base for each replicate using the 500 conditional U matrices from the MCMC models.

Whole grain density is calculated as the ratio of the usual intakes of whole grains and energy.

The %nci_multivar_summary() macro is used to calculate means and quantiles. 
Multiple variables can be summarized at the same time by specifying them as a space-delimited list in the variables parameter.
*/

%macro summary_brr(num_brr=);

	%do brr_rep = 0 %to &num_brr.;
	
		%put Starting Iteration &brr_rep.;
		
		**simulate usual intakes using conditional U matrices from MCMC;
		%nci_multivar_distrib(multivar_mcmc_model=model_brr&brr_rep.,
													distrib_population=distrib_pop,
													id=SEQN,
													weight=RepWt_&brr_rep.,
													nuisance_weight=Weekend_Weight,
													use_mcmc_u_matrices=Y,
													distrib_seed=%eval(99999 + &brr_rep.),
													outname=model_distrib_out);
													
		**calculate whole grain density;
		data model_distrib_out;
			set model_distrib_out;
			
			density_Whole_Grain = 1000 * usl_Whole_Grain/usl_Energy;
		run;
													
		**compute means, quantiles, and proportions for simulated whole grain and energy intakes;
		%nci_multivar_summary(input_data=model_distrib_out,
													variables=density_Whole_Grain
																		usl_Whole_Grain 
																		usl_Energy,
													weight=RepWt_&brr_rep.,
													do_means=Y,
													do_quantiles=Y,
													quantiles=5 25 50 75 95,
													outname=summary_brr&brr_rep.);
	%end;
	
	**extract point estimate and BRR replicates;
	data summary_brr_data;
		set summary_brr0;
		%do brr_rep = 1 %to &num_brr.;
			set summary_brr&brr_rep. (keep = value rename=(value = brr&brr_rep.));
		%end;
	run;
%mend summary_brr;

%summary_brr(num_brr=&num_brr.);

/*
Variances and confidence intervals are calculated using the BRR replicates.
*/

**calculate degrees of freedom;
proc sort data=input_dataset; by SDMVSTRA; run;

data _NULL_;
	set input_dataset end=last_obs;
	by SDMVSTRA;
	
	retain num_strata 0;;
	
	if first.SDMVSTRA then num_strata = num_strata + 1;
	
	if last_obs = 1 then call symputx("df", num_strata);
run;

**create summary report;
data summary_report (keep = population variable statistic value std_error confidence_lower confidence_upper);
	set summary_brr_data;
	
	array reps{&num_brr.} brr1-brr&num_brr.;
	
	**calculate BRR standard error;
	sum_sq_diff = 0;
	do i = 1 to &num_brr.;
	
		sum_sq_diff = sum_sq_diff + (reps{i} - value)**2;
	end;
	
	std_error = sqrt(sum_sq_diff/(&num_brr.*&fay_factor.**2));
	
	**95% confidence intervals;
	confidence_lower = value + tinv(0.025, &df.)*std_error;
	confidence_upper = value + tinv(0.975, &df.)*std_error;
run;

proc print data=summary_report; 

	title "Distribution of Whole Grain Density";
run;