/*
Title: Multiple Imputation with BRR Standard Errors
Author: NCI/Information Management Services
Date: 6/3/2025
*/

/*
This example demonstrates multiple imputation 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. Draw Parameter Sets
    
    For each parameter set:
      6. Simulate Usual Intakes
      7. Regression Model
      
    8. Average Model Coefficients
    
  9. 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 nutrient being analyzed is potassium (TPOTA).

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

The outcome variable being measured will be systolic blood pressure (BPSY_AVG).

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, 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(TPOTA) and
		 not missing(BPSY_AVG);
	
	**rename potassium variable for readability;
	Potassium = TPOTA;
run;

/*
BRR replicate weights are created for 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);
						 
/*
The potassium variable can now be transformed and standardized for use in the MCMC algorithm.
*/

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

**find best Box-Cox lambda parameter;
%boxcox_survey(input_data=input_dataset,
							 row_subset=%quote(Day2 = 0),
							 variable=Potassium,
							 weight=RepWt_0,
							 covariates=SMK_REC RIDAGEYR RIAGENDR Weekend,
							 cat_covariates=SMK_REC,
							 cat_reference=3);
							 
**Calculate minimum consumption amount;
%calculate_minimum_amount(input_data=input_dataset,
													row_subset=%quote(Day2 = 0),
													daily_variables=Potassium);
													
**transform and standardize;
%nci_multivar_preprocessor(input_data=input_dataset,
													 daily_variables=Potassium,
													 continuous_covariates=RIDAGEYR,
													 boxcox_lambda_data=bc_Potassium,
													 minimum_amount_data=minimum_amount_data,
													 outname=model);
													 
/*
MCMC models are fit for the base weight and each set of replicate weights.
Unlike in regression calibration, the outcome variable is included as a covariate in the measurement error model.

Instead of conditional U matrices drawn after the MCMC chain, multiple imputation draws U matrices from the chain itself.
The save_mcmc_u parameter is set to Y so that the U matrices created in the main MCMC chain are saved.
*/

%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=Potassium,
											 default_covariates=SMK_REC std_RIDAGEYR RIAGENDR Day2 Weekend BPSY_AVG,
											 cat_covariates=SMK_REC,
											 cat_reference=3,
											 num_mcmc_iterations=3000,
											 num_burn=1000,
											 num_thin=2,
											 save_mcmc_u=Y,
											 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.
*/

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;

/*
To perform multiple imputation, parameter sets will be drawn from the distribution of the MCMC parameters. 
This can be using the %draw_parameters() macro which will select iterations in the MCMC chain that are spaced far enough apart to ensure that they are independent.
This example will use 10 parameter sets, specified using num_draws in %draw_parameters().

A single set of usual intakes are then simulated with %nci_multivar_distrib() for each parameter set.
Since U matrices in the main MCMC chain were saved from %nci_multivar_mcmc(), a U matrix conditional on the parameters will be drawn with each parameter set.

A regression model is fit on the simulated usual intakes for each parameter set.
The model coefficients are averaged across the parameter sets and saved for each BRR replicate.
*/

%macro summary_brr(num_brr=);

	%do brr_rep = 0 %to &num_brr.;
	
		%put Starting Iteration &brr_rep.;
		
		**draw 10 sets of parameters from MCMC parameter distribution;
		%let num_sets = 10;
		%draw_parameters(multivar_mcmc_model=model_brr&brr_rep., num_draws=&num_sets.);
		
		%do set_num = 1 %to &num_sets.;
		
			**simulate 1 usual intake for each subject using the U matrix from the model parameter set;
			%nci_multivar_distrib(multivar_mcmc_model=DRAW&set_num.,
														distrib_population=distrib_pop,
														id=SEQN,
														weight=RepWt_&brr_rep.,
														nuisance_weight=Weekend_Weight,
														additional_output=SMK_REC RIDAGEYR RIAGENDR BPSY_AVG,
														use_mcmc_u_matrices=Y,
														distrib_seed=%eval(99999 + &num_sets.*&brr_rep. + &set_num.),
														outname=regression_data);
														
			**scale down simulated potassium intake by 1000 to show the effect per 1,000 mg of potassium;
			data regression_data;
				set regression_data;
				
				usl_Potassium = usl_Potassium/1000;
			run;
			
			ods select none;
			
			**fit linear model of systolic blood pressure;
			proc surveyreg data=regression_data;
			
				class SMK_REC(ref="3");
			
				model BPSY_AVG = usl_Potassium SMK_REC RIDAGEYR RIAGENDR / solution;
				
				weight RepWt_&brr_rep.;
				
				ods output ParameterEstimates=bp_model;
			run;
			
			%summary_coef_surveyreg(parameter_estimates=bp_model,
															response=BPSY_AVG,
															outname=coefficients&set_num.);
															
			ods select all;
		%end;
		
		**save average of coefficients across parameter sets;
		data coefficient_values;
			%do set_num = 1 %to &num_sets.;
				set coefficients&set_num. (keep=value rename=(value=draw&set_num.));
			%end;
			
			mean_coefficient = mean(of draw1-draw&num_sets.);
		run;
		
		data summary_brr&brr_rep.;
			set coefficients1 (keep = population variable statistic);
			set coefficient_values (keep=mean_coefficient rename=(mean_coefficient=value));
		run;
	%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 "Systolic Blood Pressure vs. Potassium Usual Intake";
run;