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

/*
This example demonstrates estimating a univariate distribution 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 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.

The NHANES dataset strata are masked pseudo-strata that hide the confidential information used to make the true design strata.
These pseudo-strata are designed to produce the same variance estimates as the true strata.

Individuals subjects are identified by SEQN. 
Each subject has up to two dietary recalls which are identified by DAY.

The nutrient being analyzed is sodium (TSODI).

The covariates being examined are smoking status (SMK_REC), age (RIDAGEYR), sex (RIAGENDR).
Two nuisance covariates are being accounted for, whether the recall was on a weekend (`Weekend`) and whether the recall is on day 2 (`Day2`). 
They are considered nuisance covariates because they are not necessarily of interest to study but must be accounted for to get an accurate distribution.

The WTDRD1 variable is the survey weighting of each observation.

The data is structured in a tall format where each row represents one observation of a subject. 
In some datasets, not all subjects will have the same number of observations. 
Not all subjects need to have the same number of observations as long as there are sufficient subjects that have two or more observations.
A good rule is to have at least 50 subjects with two or more observations where each variable of interest is non-zero.

Looking at the first 10 observations shows the structure of the data. 
*/

**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;

**First 10 observations;
proc print data=input_dataset (obs = 10); 

	var SEQN DAY TSODI SMK_REC RIDAGEYR RIAGENDR Day2 Weekend WTDRD1;
	
	title "Structure of NHANES Data";
run;

/*
The data should now be prepared for analysis.

All subjects used in the analysis must be complete cases (i.e, each subject must have all of the covariates, nutrients, and foods that are being analyzed). 
For this example, subjects who are missing one or more variables being analyzed will be removed from the dataset.
*/

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(TSODI);
	
	**rename sodium variable for readability;
	Sodium = TSODI;
run;

/*
To calculate standard errors, sets of BRR replicate weights must be generated based on the base survey weights. 
BRR is similar to bootstrap, but it is performed using a structured set of weights generated using a Hadamard matrix. 
For standard BRR, half of the sample has double weight and the other half has zero weight for each set of weights. 

The number of weight sets is determined by dimension of the Hadamard matrix, which must be 1, 2, or a multiple of 4 and greater than the number of strata.
Since this example uses a subset of the NHANES data with 6 Strata, a Hadamard matrix of dimension 8 is used, and 8 BRR replicates are required.
If the full dataset were used, there are 46 strata so 48 BRR replicates would be needed.
Since ordinary bootstrap usually requires 200 or more replicates, BRR is often far more efficient for datasets where it can be applied.

Post-stratification adjustment will also be performed so that the sum of the weights in each post-stratification cell is the same as for the base weights. 

To use BRR, each strata in the dataset must have exactly two primary sampling units (PSUs).
Most NHANES data, including the dataset used in this example, meets this criteria.

Some NHANES cycles contain unusual strata that do not have exactly two PSUs.
These strata can be corrected the %fix_nhanes_strata() utility.
The function call is demonstrated below even though all strata in the dataset used in this example have 2 PSUs.
*/

%fix_nhanes_strata(input_data=input_dataset,
									 outname=input_dataset);
									 
/*
BRR replicate weight sets can be generated using the %brr_weights() macro. 

The %brr_weights() macro takes in the strata (SDMVSTRA), PSU (SDMVPSU), base weight (WTDRD1), and post-stratification cell (PSCELL) variables. 

The output of the function will contain the base weight variable (RepWt_0) and a RepWt_ variable for each BRR weight set.
The output weights are integerized by rounding them to the nearest integer and distributing the remainders to minimize rounding error.

In this tutorial, the Fay method of BRR will be used. 
This variation of BRR adjusts the weights using a specified value called the Fay factor (f), which must be between 0 and 1.
The weights in half of the sample are multiplied by 1 + f and the weights in the other half are multiplied by 1 - f for every replicate weight set.
Using a Fay factor below 1 ensures that every observation will be used in every weight set.
A Fay factor of 1 is the same as standard BRR.

For this example, a Fay factor of 0.7 will be used. 
This multiplies the weights in half of the sample by 1.7 (1 + 0.7) and multiplies the weights in the other half by 0.3 (1 - 0.7) for every replicate weight set.
It is important to record the Fay factor used to generate the weights since it will become important when calculating the BRR variance and standard error.
*/

%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);
						 
/*
The next step is to Winsorize extreme values of the sodium variable. 
The %boxcox_survey() macro has a helpful utility to suggest cutoffs for Winsorization.

By default, extreme values are defined as being beyond the range three times the interquartile range (IQR) below the 25th or above the 75th percentile on the transformed scale. 
This can be changed using the iqr_multiple parameter which sets how many times the IQR away from 25th and 75th percentiles a value can be before being considered an outlier.
Lower IQR multiples lead to stricter Winsorization of the data since more values are identified as outliers.
For this example analysis, a cutoff of twice the IQR is used to illustrate the structure of the Winsorization report. 
For an actual analysis, the ideal Winsorization cutoff should be the least strict value that still removes problematic outliers.
This will have to be found through experimentation, though the default value of three times the IQR will usually suffice.

Note that the suggested Winsorized values are found using values in the first recall but are applied to the entire dataset. 
Additionally, only positive values can be used since the Box-Cox transformation is invalid for zero and negative values. 
The %boxcox_survey() macro will only use positive values to find the best lambda. 
Since sodium is a daily consumed nutrient, values are expected to be above zero for all subjects.

For Winsorization, the id and repeat_obs parameters are needed to specify the subject and observation of each row of the dataset.
This allows the algorithm to identify which rows are in need of Winsorization and produce a report.
*/

**run Box-Cox survey on sodium variable to get suggested Winsorization cutoffs using the weighting for the point estimate;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Sodium,
							 weight=RepWt_0,
							 do_winsorization=Y,
							 iqr_multiple=2,
							 id=SEQN,
							 repeat_obs=DAY);
							 
/*
The Winsorization report has identified extreme values of Sodium as being below 387.79 mg or above 13598.56 mg. 
Using these values, it is now possible to Winsorize the Sodium variable.
*/

data input_dataset;
	set input_dataset;
	
	if Sodium < 387.7924 then Sodium = 387.7924;
	if Sodium > 13598.5605 then Sodium = 13598.5605;
run;

/*
The NCI method assumes that food and nutrient amounts have a normal distribution. 
The Sodium variable can be normalized using a Box-Cox transformation. 
A suitable transformation parameter can be found using the same %boxcox_survey() macro used to find the Winsorization cutoffs. 
Unlike with Winsorization, the function will be run using the covariates that will be used in the model.

The covariates are specified as a space-delimited list in the covariates argument.
Categorical covariates are specified in cat_covariates and their corresponding reference values in cat_reference.
Smoking status (SMK_REC) will be categorical using never-smokers (3) as the reference.

Setting do_plots to Y generates Q-Q plots of the original scale and Box-Cox transformed variable against the standard normal distribution.
The solid line represents the hypothetical quantiles of the variable if it were normally distributed.
The points will be close to the line if the actual distribution of the variable is close to normal.
*/

ods listing close;

ods pdf file="daily_nutrient_boxcox_plots.pdf";

**run Box-Cox survey with covariates to find best lambda parameter;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Sodium,
							 weight=RepWt_0,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3,
							 do_plots=Y);
							 
ods pdf close;
							 
ods listing;
							 
/*
Next, the minimum amount for sodium is calculated using the %calculate_minimum_amount() macro. 
The minimum amount is calculated as half the smallest non-zero intake in the dataset.
This value becomes the minimum allowed usual intake when simulating usual intakes later in the workflow.
As with the Box-Cox lambda, the minimum amount is found using the first recall.
*/

%calculate_minimum_amount(input_data=input_dataset,
													row_subset=%quote(Day2 = 0),
													daily_variables=Sodium);
													
/*
Now, the %nci_multivar_preprocessor() macro can be used to generate the input dataset for the MCMC model. 
The preprocessor function will apply the Box-Cox transformation from %boxcox_survey() to the sodium variable.
The transformed variable is then standardized to a mean of zero and a variance of 2 as required by the MCMC algorithm. 

All continuous covariates (such as age) should be standardized to a mean of zero and a variance of 1 for best results in the MCMC algorithm. 
The covariates to be standardized should be given in continuous_covariates as a space-delimited list.
The standardized covariates created by %nci_multivar_preprocessor() have the prefix std_.
*/

%nci_multivar_preprocessor(input_data=input_dataset,
													 daily_variables=Sodium,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=bc_Sodium,
													 minimum_amount_data=minimum_amount_data,
													 outname=model);
													 
/*
Before fitting the MCMC model, it is worthwhile to explore the parts of the MCMC input.

The mcmc_in element is the input dataset with the transformed and standardized variables and covariates added. 

The backtran element contains the data used to transform the sodium variable. 
This data is necessary to interpret the MCMC results and will be utilized later when backtransforming simulated usual intakes to the original scale.

It is also important to be familiar with the backtran data columns.

	tran_lambda: The Box-Cox lambda used to transform the variable.
	minamount: The minimum allowed usual intake, defined as half of the smallest non-zero intake in the observed data.
	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 indicator of whether the variable is a biomarker assumed to be unbiased on the transformed scale.
						 If 0, a bias correction factor will be added and a 9-point approximation will be used for backtransformation.
						 If 1, an exact backtransformation will be used with no correction.
*/

proc print data=model_backtran; 

	title "Backtransformation Data";
run;

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

Covariates to be used for every variable in the model are specified in default_covariates as a space-delimited list.
Note that the standardized version of the age covariate generated by the preprocessor (std_RIDAGEYR) is used in the MCMC model.
Like %boxcox_survey(), categorical covariates are specified in cat_covariates and their reference values in cat_reference.

The number of iterations to run the MCMC must also be selected.
The iteration parameters are:

  num_mcmc_iterations: The total number of iterations to run the MCMC including burn-in
	num_burn: The number of burn-in iterations that will be discarded when calculating posterior means
	num_thin: The spacing between iterations that will be used to calculate posterior means to ensure independence between iterations
	
It is important to choose iteration parameters so that there are enough thinned samples to obtain stable estimates of the posterior means.
A good starting point is to make sure that the number of thinned iterations, calculated as (num_mcmc_iterations - num_burn)/num_thin, is at least 400.

The num_post argument specifies how many conditional random effect (U) matrix draws should be made after the MCMC is finished.
These U matrices are conditional on the recalls and MCMC parameters for the subjects in the MCMC.

Conditional U matrices can be used for generating usual intakes if the same population is used to fit the MCMC model and generate usual intakes.
Using conditional U matrices does not substantially change the results for distributions, but it is required for regression calibration.
It can be helpful to output conditional U matrices even for distributions to make it easier to add regression calibration later on without re-running the MCMC.

Drawing more conditional U matrices will give a more stable distribution at the cost of using more memory to store the U matrices.
A good starting point is between 200-500 conditional U matrix draws.

The random seed for the MCMC is given in mcmc_seed.
The seed should be changed each time the MCMC is run within the analysis to avoid biasing the results.
Re-running the MCMC models using the same sequence of seeds will produce the same results each time.

For BRR to work properly, nearly all BRR replicates should be used in variance calculation.
It is good practice to verify that all of the BRR replicates have converged using techniques such as trace plots and the Gelman-Rubin test.
Replicates that did not converge should not be included in the variance calculation.
*/

%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.,
											 daily_variables=Sodium,
											 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,
											 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 MCMC model estimates the following parameters:

  Beta: Coefficients of the covariates and intercept in the model
  Sigma-u: Between-person variance-covariance matrix
  Sigma-e: Within-person variance-covariance matrix

Each iteration of the MCMC algorithm represents a sample from the distribution of the parameters.
The posterior means are the means of the parameter samples at the thinned iterations after burn-in.
They are used as the point estimates of the MCMC parameters.

The posterior means of each parameter can be viewed in the model output.
*/

proc print data=model_brr0_beta1_mean; 

	title "Beta";
run;
proc print data=model_brr0_sigma_u_mean; 

	title "Sigma-u";
run;
proc print data=model_brr0_sigma_e_mean; 

	title "Sigma-e";
run;

/*
The next step is simulating usual intakes using the MCMC model parameters.

A population-based dataset must be created that contains all of the subjects and additional variables that will be used for generating usual intakes and downstream analysis. 
For this example, the population base will be the first instance of each subject in the MCMC input dataset.
*/

**get first instance of each subject;
proc sort data=model_mcmc_in; by SEQN; run;

data distrib_pop;
	set model_mcmc_in;
	by SEQN;
	if first.SEQN; 
run;

/*
Once the initial population base is created, nuisance variables must be accounted for.
To factor out the effect of the recall being conducted on day 2, the Day2 variable will be set to zero for all subjects.
*/

data distrib_pop;
	set distrib_pop;
	
	Day2 = 0;
run;

/*
To account for weekend vs weekday consumption, the simulated usual intakes for weekends and weekdays will be weighted and averaged for each subject.
To accomplish this, a repeat of each subject will be created in the population base. 
The first instance of each subject corresponds to weekday consumption (Weekend = 0) and the second instance corresponds to weekend consumption (Weekend = 1).
Since weekends are defined as three days in this dataset (Friday, Saturday, and Sunday), weekend consumption will be given a weight of 3 while weekday consumption will be given a weight of 4.
*/

**create a repeat of each subject and set weekend indicators and weights;
data distrib_pop;
	set distrib_pop;
	
	Weekend = 0;
	Weekend_Weight = 4;
	output;
	
	Weekend = 1;
	Weekend_Weight = 3;
	output;
run;

/*
The population base and MCMC parameters are used to create a distribution of simulated usual intakes using %nci_multivar_distrib().

The subjects in the MCMC and population base are the same, so the conditional U matrices from the MCMC can be used for generating usual intakes.
This done by setting use_mcmc_u_matrices to Y.
Since there are 500 conditional U matrices, 500 usual intakes will be generated for each subject in the population base.

The id and weight parameters specify the subject identifier and weight from the population base to include in the output.
For population bases with multiple nuisance variable levels (such as Weekend), nuisance_weight specifies the weight for each level.

The random seed for usual intake simulation is given by distrib_seed.
Re-simulating the usual intakes with the same seeds and the same MCMC models will produce the same simulated values.
As with the MCMC, the seed should be changed for each call to simulate usual intakes in the analysis to avoid biasing the results.

Summary statistics can then be calculated from the simulated usual intakes using the %nci_multivar_summary() macro.
This is done in the same loop as the usual intake simulation so that the usual intake distribution dataset can be overwritten each time to save memory.

The proportions are specified as SAS datasets for the lower and upper intake thresholds.
The variable column represents the variable(s) to calculate proportions for while the threshold column represents the lower or upper intake limits.
	
The lower_thresholds parameter is used for calculating the proportions of intakes above a lower limit.
Likewise, upper_thresholds is used for calculating the proportions of intakes below an upper limit.

This process is repeated for the base weight and each set of replicate weights.
*/

**SAS datasets for lower and upper sodium intake thresholds;
data lower_thresholds;
			
	variable = "usl_Sodium";
	threshold = 2200;
	output;
run;
		
data upper_thresholds;

	variable = "usl_Sodium";
	threshold = 3600;
	output;
run;

%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);
													
		**compute means, quantiles, and proportions for simulated sodium intakes;
		%nci_multivar_summary(input_data=model_distrib_out,
													variables=usl_Sodium,
													weight=RepWt_&brr_rep.,
													do_means=Y,
													do_quantiles=Y,
													quantiles=5 25 50 75 95,
													do_proportions=Y,
													lower_thresholds=lower_thresholds,
													upper_thresholds=upper_thresholds,
													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.);

/*
With a point estimate and BRR replicates for every summary statistic, standard errors and confidence intervals can now be calculated. 
With the data set up as one column per replicate, this can be done easily and efficiently using a SAS DATA step.

If using Fay's method, the variance of the BRR replicates is divided by the square of the Fay factor.
It is important to verify the Fay factor used to generate the BRR weights so that the correct variance can be calculated.
If ordinary bootstrapping is used instead of BRR, the same standard error formula can be used with the Fay factor set to 1.

The number of degrees of freedom for confidence intervals is the total number of PSUs minus the number of strata.
Since BRR uses exactly two PSUs per strata, the number of degrees of freedom is simply the number of strata.
*/

**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 Sodium Usual Intake";
run;

