/*
Title: Using the NCI Method in Low-Memory Environments
Author: Amer Moosa
Date: 2/3/2026
*/

/*
This example demonstrates an NCI method workflow using the low-memory features of the NCI method macros.

Workflow:

  1. Data Cleaning
  2. Pre-Processing
  3. MCMC Model
  4. Simulate Usual Intakes
  5. Calculate Ratios
  6. Regression Calibration
*/

libname indata "../data";

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

/*
The NCI method can use a lot of computer memory, especially for MCMC model parameters and simulated usual intake datasets.
Fortunately, the NCI method macros have built-in ways to reduce the amount of data held in memory.
These features often come at the cost of computation speed, but they can be useful for large models or when running the software on systems with fewer resources.

The effect of the ratio of sodium (TSODI) to potassium (TPOTA) on systolic blood pressure (BPSY_AVG) will be measured from 2005-2010 NHANES data.
To calculate the sodium/potassium ratio, a bivariate NCI method model of sodium and potassium will be fit.
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 survey weighting of each subject.

A simple way to conserve memory is subsetting the data to only the variables needed for the analysis.
While SAS datasets are stored on the disk, they are loaded into memory for some PROCs including those used by NCI method macros.

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

**subset data;
data input_dataset (keep = SEQN DAY
													 TSODI TPOTA
													 SMK_REC RIDAGEYR RIAGENDR Weekend Day2
													 BPSY_AVG
													 WTDRD1);
	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, variables, or outcomes;
	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)	 	 and
		 not missing(TPOTA)		 and
		 not missing(BPSY_AVG);
	
	**rename sodium and potassium variables for readability;
	Sodium = TSODI;
	Potassium = TPOTA;
run;

/*
The sodium and potassium variables can now be transformed and standardized for use in the MCMC algorithm.
*/

**Winsorize extreme values of sodium and potassium intake;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Sodium,
							 weight=WTDRD1,
							 do_winsorization=Y,
							 id=SEQN,
							 repeat_obs=DAY);
							 
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Potassium,
							 weight=WTDRD1,
							 do_winsorization=Y,
							 id=SEQN,
							 repeat_obs=DAY);
							 
data input_dataset;
	set input_dataset;
	
	if Sodium < 113.1129 then Sodium = 113.1129;
	
	if Potassium < 42.45263 then Potassium = 42.45263;
run;

**run Box-Cox survey and create Box-Cox lambda data using the first recall;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Sodium,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3,
							 weight=WTDRD1);
							 
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Potassium,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3,
							 weight=WTDRD1);
							 
data boxcox_lambda_data;
	set bc_Sodium
			bc_Potassium;
run;

**Calculate minimum amount of sodium and potassium intake in the first recall;
%calculate_minimum_amount(input_data=input_dataset,
													row_subset=%quote(Day2 = 0),
													daily_variables=Sodium Potassium);
													
**Run MCMC pre-preprocessor;
%nci_multivar_preprocessor(input_data=input_dataset,
													 daily_variables=Sodium Potassium,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=boxcox_lambda_data,
													 minimum_amount_data=minimum_amount_data,
													 outname=model);
													 
/*
The MCMC measurement error model can now be fit using %nci_multivar_mcmc(), drawing 500 conditional U matrices.

The MCMC parameters from each iteration and U matrices are saved in memory by default.
Covariance matrix size increases with the square of the number of variables, and U matrix size increases with the number of subjects.
This means that memory requirements can become very high for multivariate models and large sample sizes.

The low_memory_mode option in %nci_multivar_mcmc() can be used to reduce the memory requirements of the algorithm.
Setting low_memory_mode to Y writes parameters and U matrices to SAS datasets at each iteration instead of holding them in memory.
The memory usage of the algorithm is significantly reduced at the cost of slower run time due to the extra I/O overhead.
*/

%nci_multivar_mcmc(pre_mcmc_data=model,
									 id=SEQN,
									 repeat_obs=DAY,
									 weight=WTDRD1,
									 daily_variables=Sodium Potassium,
									 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=2,
									 num_post=500,
									 mcmc_seed=9999,
									 outname=model,
									 low_memory_mode=Y);
									 
/*
The population base for simulating usual intakes is generated from the input data.
*/

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;

/*
The %nci_multivar_distrib() macro can now be used to generate 500 usual intakes for each subject using the conditional U matrices from the MCMC.
The outcome variable and covariates for the regression are included in the output dataset using additional_output.

Just like the MCMC model, the size of the simulated usual intake dataset increases with the number of subjects and variables.

The %nci_multivar_distrib() macro also contains a low_memory_mode option to reduce memory use.
When low_memory_mode is Y, each simulated usual intake replicate is written to the output dataset instead of keeping it in memory.
As with the MCMC algorithm, memory usage is reduced significantly in exchange for a slower run time.
*/

**simulate usual intakes for each subject using MCMC U matrices;
%nci_multivar_distrib(multivar_mcmc_model=model,
											distrib_population=distrib_pop,
											id=SEQN,
											weight=WTDRD1,
											nuisance_weight=Weekend_Weight,
											use_mcmc_u_matrices=Y,
											additional_output=SMK_REC RIDAGEYR RIAGENDR BPSY_AVG,
											distrib_seed=99999,
											outname=model_distrib_out,
											low_memory_mode=Y);
											
/*
The usual intakes of sodium and potassium are used to calculate the sodium/potassium ratio.
*/

**Calculate sodium/potassium ratio;
data model_distrib_out;
	set model_distrib_out;
	
	Sodium_Potassium_Ratio = usl_Sodium/usl_Potassium;
run;

/*
The distribution of the sodium/potassium ratio is summarized by using %nci_multivar_summary().

Calculating quantiles for a variable requires loading the full column into memory.
As a result, memory requirements are higher if quantile calculation is enabled.
If quantiles are not needed, do_quantiles can be set to N.
*/

%nci_multivar_summary(input_data=model_distrib_out,
											variables=Sodium_Potassium_Ratio,
											weight=WTDRD1,
											do_means=Y,
											do_quantiles=Y,
											quantiles=5 25 50 75 95,
											outname=summary_stats);
											
proc print data=summary_stats;

	title "Distribution of Sodium/Potassium Ratio";
run;

/*
For regression calibration, the sodium/potassium ratio is averaged for each subject to calculate the proxy intakes.
*/

**Average sodium/potassium ratio for each subject;
proc sort data=model_distrib_out; by SEQN; run;

proc means data=model_distrib_out noprint;
	by SEQN;
	
	var Sodium_Potassium_Ratio;
	
	id SMK_REC RIDAGEYR RIAGENDR BPSY_AVG WTDRD1;
	
	output out=regression_data mean=Sodium_Potassium_Ratio;
run;

/*
The linear model is fit on the dataset of proxy intakes.
*/

ods select none;

**fit linear model;
proc surveyreg data=regression_data;

	class SMK_REC(ref="3");

	model BPSY_AVG = Sodium_Potassium_Ratio SMK_REC RIDAGEYR RIAGENDR / solution;
	
	weight WTDRD1;
	
	ods output ParameterEstimates=bp_model;
run;

ods select all;

**summary dataset of model parameters;
%summary_coef_surveyreg(parameter_estimates=bp_model,
												response=BPSY_AVG,
												outname=bp_parameters);
												
proc print data=bp_parameters; 

	title "Systolic Blood Pressure vs. Sodium/Potassium Ratio";
run;