/*****************************************************************************/
/*****************************************************************************/
/* The following files are required to run this example:                     */
/*   boxcox_survey.macro.v1.2.sas                                            */
/*   std_cov_boxcox24hr_conday_minamt_macro_v2.0.sas                         */
/*   multivar_mcmc_macro_v2.1.sas                                            */
/*   multivar_distrib_macro_v2.1.sas                                         */
/*   simulated_example_data.sas7bdat                                         */
/*****************************************************************************/
/*                                                                           */
/*****************************************************************************/
/* This example analysis uses regression calibration and fits a logistic     */
/* regression model to assess the relationships between seven dietary        */
/* components and a health outcome.  For this analysis, the example data     */
/* were simulated to include data from 100,000 participants, and the 24-hour */
/* dietary recall is the main dietary instrument.  This analysis allows for  */
/* participants that never consume alcohol.  The dietary components are      */
/* alcohol, red meat, whole grains, sodium, potassium, protein, and energy,  */
/* and the health outcome is a binary event variable.  The simulated data    */
/* include repeated intake measurements from 24-hour dietary recalls.        */
/*                                                                           */
/* This analysis uses bootstrap variance estimation, so the multivariate     */
/* measurement error model and logistic model are fit using the original     */
/* data set and 200 replicate data sets.  The following replicfirst and      */
/* repliclast macro variables allow specification of a range of replicates.  */
/* This example program considers the original data set (i.e. replicate 0).  */
/* When the replicfirst macro variable is set to 0, as it is in this         */
/* example, the 200 replicate data sets are generated using the SAS          */
/* Surveyselect procedure, and these 200 replicate data sets are combined    */
/* with the original data set and are saved for later use.                   */
/*****************************************************************************/

%let replicfirst = 0;
%let repliclast  = 0;



title1 "Fit Multivariate Measurement Error Model Using MCMC with 24-Hour Recall as Main Instrument, Allowing for Never Consumers";
title2 "Predict Intake and Perform Regression Calibration";
title3 "Assess Relationships between Seven Dietary Components and a Health Outcome";



***********************************************************************;
*** The path must be modified in the following lines to specify the ***;
*** location of the SAS macros, the SAS data files, and the SAS     ***;
*** library that will be used for output                            ***;
***********************************************************************;

%let home = /prj/dcp/statprog/meas.err/develop.public.resources.stat.meth;

*** Include the required macros ***;
%include "&home/include.files.macros/boxcox_survey.macro.v1.2.sas";
%include "&home/include.files.macros/std_cov_boxcox24hr_conday_minamt_macro_v2.0.sas";
%include "&home/include.files.macros/multivar_mcmc_macro_v2.1.sas";
%include "&home/include.files.macros/multivar_distrib_macro_v2.1.sas";

*** Input data library ***;
libname inlib "&home/data";

*** Output data library ***;
libname outlib "&home/multivar_epidemiology_example5_mcmc_nc_main24hr/outlib";



********************************************************************;
*** Comparison values, t0 and t1, for each risk model covariate. ***;
*** The values, t0 and t1, are the 25th and 75th percentiles for ***;
*** each covariate except alcohol.                               ***;
********************************************************************;

%let t0_alcohol = 0;
%let t1_alcohol = 14;

%let t0_redmeat = 46;
%let t1_redmeat = 114;

%let t0_whlgrain = 0.49;
%let t1_whlgrain = 1.44;

%let t0_sodium = 3431;
%let t1_sodium = 4779;

%let t0_potassium = 2913;
%let t1_potassium = 4135;

%let t0_protein = 79;
%let t1_protein = 106;

%let t0_energy = 1915;
%let t1_energy = 2666;

%let t0_std_entry_age = -0.8637;
%let t1_std_entry_age =  0.8406;

%let t0_std_log_bmi = -0.6309;
%let t1_std_log_bmi =  0.6162;



************************************************************************************;
*** The macro fit_models_replicate_loop is used to call the MULTIVAR_MCMC and    ***;
*** MULTIVAR_DISTRIB macros using the original study data (i.e. replicate data   ***;
*** set 0) and each of the replicate data sets for bootstrap variance estimation ***;
*** (i.e. replicate data sets 1, 2, ..., 200)                                    ***;
************************************************************************************;

%macro fit_models_replicate_loop;


  **************************************************************************;
  *** For replicate 0 create bootstrap data sets and save them to outlib ***;
  **************************************************************************;

  %if &replicfirst = 0 %then %do;

    *********************************************************************************;
    *** Import the simulated data file for this example and set it as replicate 0 ***;
    *********************************************************************************;

    data simdata;
      set inlib.simulated_example_data;
      replicate = 0;
    run;



    **************************************************************************************************;
    *** Use proc surveyselect to generate 200 replicate bootstrap data sets for the simulated data ***;
    **************************************************************************************************;

    proc surveyselect data=simdata(drop = replicate) out=bootdata seed=1563783 method=urs sampsize=100000 outhits reps=200;
    run;



    **********************************************************;
    *** Combine the simulated data with the bootstrap data ***;
    **********************************************************;

    data simdata_and_bootdata;
      set simdata bootdata;
    run;



    ************************************************************************;
    *** Create the unique record ID variable for each replicate data set ***;
    ************************************************************************;

    proc sort data=simdata_and_bootdata;
      by replicate;
    run;

    data simdata_and_bootdata;
      set simdata_and_bootdata;
      by replicate;

      retain replicaterowid;

      if first.replicate then replicaterowid = 1;
      else replicaterowid = replicaterowid + 1;
    run;



    **************************************************;
    *** Save the simulated data and bootstrap data ***;
    **************************************************;

    data outlib.simdata_and_bootdata(keep = replicate replicaterowid r1_alcohol_g r2_alcohol_g r1_redmeat_g r2_redmeat_g
                                            r1_whlgrain_ozeq r2_whlgrain_ozeq r1_sodium_mg r2_sodium_mg r1_potassium_mg
                                            r2_potassium_mg r1_protein_g r2_protein_g r1_energy_kcal r2_energy_kcal
                                            ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                                            std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg
                                            std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi
                                            constant1 ybin_event);
      set simdata_and_bootdata;
    run;

  %end;

  %else %do;

    *******************************************************************;
    *** For other replicates, read in the saved replicate data sets ***;
    *******************************************************************;

    data simdata_and_bootdata;
      set outlib.simdata_and_bootdata;
    run;

  %end;



  %do replicnum = &replicfirst %to &repliclast;

    title4 "Replicate &replicnum";
    data select_replicate;
      set simdata_and_bootdata (where = (replicate = &replicnum));
    run;


    ************************************************************************************;
    *** Create separate data sets for the MULTIVAR_MCMC and MULTIVAR_DISTRIB macros. ***;
    ************************************************************************************;

    data datamrec(keep = replicaterowid repeat r_alcohol_g r_redmeat_g r_whlgrain_ozeq r_potassium_mg r_protein_g r_sodium_mg
                         r_energy_kcal ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                         std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g
                         std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi constant1)

         subj1rec(keep = replicaterowid r1_alcohol_g r2_alcohol_g r1_redmeat_g r2_redmeat_g r1_whlgrain_ozeq r2_whlgrain_ozeq
                         r1_potassium_mg r2_potassium_mg r1_protein_g r2_protein_g r1_sodium_mg r2_sodium_mg r1_energy_kcal
                         r2_energy_kcal ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                         std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g
                         std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi ybin_event constant1);

      set select_replicate;

      ********************************************;
      *** Output record for 1st 24-hour recall ***;
      ********************************************;

      if nmiss(of r1_alcohol_g r1_redmeat_g r1_whlgrain_ozeq r1_sodium_mg r1_potassium_mg r1_protein_g r1_energy_kcal) = 0 then do;

        repeat = 1;

        r_alcohol_g     = r1_alcohol_g;
        r_redmeat_g     = r1_redmeat_g;
        r_whlgrain_ozeq = r1_whlgrain_ozeq;
        r_sodium_mg     = r1_sodium_mg;
        r_potassium_mg  = r1_potassium_mg;
        r_protein_g     = r1_protein_g;
        r_energy_kcal   = r1_energy_kcal;

        output datamrec;

      end;

      ********************************************;
      *** Output record for 2nd 24-hour recall ***;
      ********************************************;

      if nmiss(of r2_alcohol_g r2_redmeat_g r2_whlgrain_ozeq r2_sodium_mg r2_potassium_mg r2_protein_g r2_energy_kcal) = 0 then do;

        repeat = 2;

        r_alcohol_g     = r2_alcohol_g;
        r_redmeat_g     = r2_redmeat_g;
        r_whlgrain_ozeq = r2_whlgrain_ozeq;
        r_sodium_mg     = r2_sodium_mg;
        r_potassium_mg  = r2_potassium_mg;
        r_protein_g     = r2_protein_g;
        r_energy_kcal   = r2_energy_kcal;

        output datamrec;

      end;


      ****************************************;
      *** Output 1 record for each subject ***;
      ****************************************;

      output subj1rec;

    run;



    *******************************************************************************************;
    *** Create a data set containing only positive values for each episodic recall variable ***;
    *******************************************************************************************;

    data subj1rec_pos_alcohol subj1rec_pos_redmeat subj1rec_pos_whlgrain;
      set subj1rec;

      if r1_alcohol_g     > 0 then output subj1rec_pos_alcohol;
      if r1_redmeat_g     > 0 then output subj1rec_pos_redmeat;
      if r1_whlgrain_ozeq > 0 then output subj1rec_pos_whlgrain;
    run;



    *******************************************************************************;
    *** Call the BOXCOX_SURVEY macro to find the optimal Box-Cox transformation ***;
    *** paramater, lambda, for each recall variable                             ***;
    *******************************************************************************;

    %boxcox_survey(data    = subj1rec_pos_alcohol,
                   subject = replicaterowid,
                   var     = r1_alcohol_g,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec_pos_redmeat,
                   subject = replicaterowid,
                   var     = r1_redmeat_g,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec_pos_whlgrain,
                   subject = replicaterowid,
                   var     = r1_whlgrain_ozeq,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec,
                   subject = replicaterowid,
                   var     = r1_sodium_mg,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec,
                   subject = replicaterowid,
                   var     = r1_potassium_mg,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec,
                   subject = replicaterowid,
                   var     = r1_protein_g,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = subj1rec,
                   subject = replicaterowid,
                   var     = r1_energy_kcal,
                   covars  = ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                             std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   weight  = ,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );



    ****************************************************************************************;
    *** Create dataset with lambda values for the STD_COV_BOXCOX24HR_CONDAY_MINAMT macro ***;
    ****************************************************************************************;

    data boxcox_lambda_data(keep = tran_lambda tran_paramindex);

      set subj1rec_pos_alcohol(obs=1 in=alcohol)
          subj1rec_pos_redmeat(obs=1 in=redmeat)
          subj1rec_pos_whlgrain(obs=1 in=whlgrain)
          subj1rec(obs=1 in=subj1rec);

      if alcohol then do;
        tran_lambda = lambda_r1_alcohol_g;
        tran_paramindex = 1;
        output;
      end;

      else if redmeat then do;
        tran_lambda = lambda_r1_redmeat_g;
        tran_paramindex = 2;
        output;
      end;

      else if whlgrain then do;
        tran_lambda = lambda_r1_whlgrain_ozeq;
        tran_paramindex = 3;
        output;
      end;

      else if subj1rec then do;
        tran_lambda = lambda_r1_sodium_mg;
        tran_paramindex = 4;
        output;

        tran_lambda = lambda_r1_potassium_mg;
        tran_paramindex = 5;
        output;

        tran_lambda = lambda_r1_protein_g;
        tran_paramindex = 6;
        output;

        tran_lambda = lambda_r1_energy_kcal;
        tran_paramindex = 7;
        output;
      end;

    run;



    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* Description of the STD_COV_BOXCOX24HR_CONDAY_MINAMT MACRO                 */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /*                                                                           */
    /* The STD_COV_BOXCOX24HR_CONDAY_MINAMT macro is used to standardize         */
    /* continuous covariates, transform and standardize reported nonzero amounts */
    /* (e.g. from 24-hour recalls), and create variables needed for the          */
    /* MULTIVAR_MCMC macro.                                                      */
    /*                                                                           */
    /* Continuous covariates are standardized to have a mean of 0 and a variance */
    /* of 1, and these new standardized variables are named by adding the prefix */
    /* "std_" to the original variable name.                                     */
    /*                                                                           */
    /* For dietary components consumed episodically, consumption-day indicator   */
    /* variables are created, and these new variables are named by adding the    */
    /* prefix "conday_" to the original variable name.  Also, the following      */
    /* transformation approach, used by Zhang et al. (2011a, 2011b), is          */
    /* performed.  First, a suitable Box-Cox transformation is applied to the    */
    /* nonzero values.  Second, the Box-Cox transformed values are centered and  */
    /* scaled to obtain an approximate normal distribution with a mean of 0 and  */
    /* a variance of 2.  These new variables are named by adding the prefix      */
    /* "stdbc_" to the original variable name.                                   */
    /*                                                                           */
    /* For dietary components that are consumed nearly every day, as specified   */
    /* by the user, zero values are replaced with minimum amount values provided */
    /* by the user or calculated as half of the observed nonzero minimum amount. */
    /* Then the following transformation approach, used by Zhang et al. (2011a,  */
    /* 2011b), is performed for the dietary components that are non-episodically */
    /* consumed.  First, a suitable Box-Cox transformation is applied.  Second,  */
    /* the Box-Cox transformed values are centered and scaled to obtain an       */
    /* approximate normal distribution with a mean of 0 and a variance of 2.     */
    /* These new variables are named by adding the prefix "stdbc_" to the        */
    /* original variable name.                                                   */
    /*                                                                           */
    /* References:                                                               */
    /*                                                                           */
    /*   Zhang S, Krebs-Smith SM, Midthune D, Perez A, Buckman DW, Kipnis V,     */
    /*   Freedman LS, Dodd KW, Carroll RJ. Fitting a bivariate measurement error */
    /*   model for episodically consumed dietary components. Int J Biostat       */
    /*   2011;7(1):Article 1.                                                    */
    /*                                                                           */
    /*   Zhang S, Midthune D, Guenther PM, Krebs-Smith SM, Kipnis V, Dodd KW,    */
    /*   Buckman DW, Tooze JA, Freedman L, Carroll RJ. A new multivariate        */
    /*   measurement error model with zero-inflated dietary data, and its        */
    /*   application to dietary assessment. Ann Appl Stat 2011 Jun;5(2B):        */
    /*   1456-87.                                                                */
    /*                                                                           */
    /*                                                                           */
    /* The syntax for calling the STD_COV_BOXCOX24HR_CONDAY_MINAMT macro is:     */
    /*                                                                           */
    /* %std_cov_boxcox24hr_conday_minamt(data                            =,      */
    /*                                   prestand_continuous_covars      =,      */
    /*                                   rec24hr_epis_vars               =,      */
    /*                                   rec24hr_daily_vars              =,      */
    /*                                   boxcox_tran_lambda_data         =,      */
    /*                                   minamount_data                  =,      */
    /*                                   print                           =,      */
    /*                                   titles                          =       */
    /*                                   );                                      */
    /*                                                                           */
    /*  where                                                                    */
    /*                                                                           */
    /*  "data"                        Specifies an input data set that includes  */
    /*                                one or more observations for each subject. */
    /*                                                                           */
    /*  "prestand_continuous_covars"  Specifies a list of continuous covariates  */
    /*                                that will be standardized to have a mean   */
    /*                                of 0 and a variance of 1.                  */
    /*                                                                           */
    /*  "rec24hr_epis_vars"           Specifies a list of 24-hour recall         */
    /*                                variables for dietary components consumed  */
    /*                                episodically.                              */
    /*                                                                           */
    /*  "rec24hr_daily_vars"          Specifies a list of 24-hour recall         */
    /*                                variables for dietary components consumed  */
    /*                                every day or nearly every day.             */
    /*                                                                           */
    /*  "boxcox_tran_lambda_data"     Specifies an input data set that includes  */
    /*                                the following two variables:               */
    /*                                "tran_paramindex"                          */
    /*                                    an index value of 1, 2, ..., Q where   */
    /*                                    Q is the total number of dietary       */
    /*                                    components specified.  The values 1,   */
    /*                                    2, ..., Q should be assigned according */
    /*                                    to the order of the dietary components */
    /*                                    as specified for the                   */
    /*                                    "rec24hr_epis_vars" and                */
    /*                                    "rec24hr_daily_vars" macro parameters, */
    /*                                    and the "rec24hr_epis_vars" dietary    */
    /*                                    components should precede the          */
    /*                                    "rec24hr_daily_vars" dietary           */
    /*                                    components.                            */
    /*                                "tran_lambda"                              */
    /*                                    a Box-Cox transformation parameter     */
    /*                                    value for the corresponding dietary    */
    /*                                    component.                             */
    /*                                The records in this data set should be     */
    /*                                sorted according to the order described    */
    /*                                for the "tran_paramindex" variable.        */
    /*                                                                           */
    /*  "minamount_data"              Specifies an optional input data set that  */
    /*                                includes the following two variables:      */
    /*                                "tran_paramindex"                          */
    /*                                    an index value of 1, 2, ..., Q where   */
    /*                                    Q is the total number of dietary       */
    /*                                    components specified.  The values 1,   */
    /*                                    2, ..., Q should be assigned according */
    /*                                    to the order described for the         */
    /*                                    "boxcox_tran_lambda_data" macro        */
    /*                                    parameter.                             */
    /*                                "minamount"                                */
    /*                                    a minimum amount value for the         */
    /*                                    corresponding dietary component.       */
    /*                                The records in this data set should be     */
    /*                                sorted according to the order described    */
    /*                                for the "tran_paramindex" variable.        */
    /*                                                                           */
    /*  "print"                       If "print=y" or "print=Y" then macro       */
    /*                                results are printed.  The default value is */
    /*                                "y".                                       */
    /*                                                                           */
    /*  "titles"                      Specifies the number of title lines to be  */
    /*                                reserved for the user's titles.  Two       */
    /*                                additional title lines are used by the     */
    /*                                macro.  The default value is "0".          */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /* Macro Output:                                                             */
    /*                                                                           */
    /* The new variables and the original variables are saved in a SAS data set  */
    /* named "stdcov_stdbc24hr_conday_out" which can be used as the input data   */
    /* set for the MULTIVAR_MCMC macro.                                          */
    /*                                                                           */
    /* The following global macro variables are declared and used to output      */
    /* variable lists that can be used to specify the lists of input variables   */
    /* needed for the MULTIVAR_MCMC macro:                                       */
    /*     conday_var_list                                                       */
    /*     stdbc_epis_var_list                                                   */
    /*     stdbc_daily_var_list                                                  */
    /*     std_continuous_covar_list.                                            */
    /*                                                                           */
    /* The macro also saves the following variables in a SAS data set named      */
    /* "backtran_out" which can be used in subsequent analysis steps that        */
    /* require back-transformation:                                              */
    /*     tran_paramindex tran_lambda tran_center tran_scale minamount.         */
    /*                                                                           */
    /*****************************************************************************/
    /*****************************************************************************/
    /*****************************************************************************/

    %std_cov_boxcox24hr_conday_minamt(data                       = datamrec,
                                      prestand_continuous_covars = ,
                                      rec24hr_epis_vars          = r_alcohol_g r_redmeat_g r_whlgrain_ozeq,
                                      rec24hr_daily_vars         = r_sodium_mg r_potassium_mg r_protein_g r_energy_kcal,
                                      boxcox_tran_lambda_data    = boxcox_lambda_data,
                                      minamount_data             = ,
                                      print                      = y,
                                      titles                     = 4
                                      );



    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* Description of the MULTIVAR_MCMC MACRO                                    */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /*                                                                           */
    /* The MULTIVAR_MCMC macro is used to fit a multivariate measurement error   */
    /* model that can include episodically consumed dietary components and       */
    /* non-episodically consumed dietary components.  The measurement error      */
    /* model is fit using the Markov Chain Monte Carlo (MCMC) algorithm from     */
    /* Zhang et al. (2011a, 2011b).  This methodology and macro also allow for   */
    /* survey-weighted MCMC computations.  Output from this macro can be saved   */
    /* and used in subsequent analysis steps such as the estimation of           */
    /* multivariate usual intake distributions and the use of regression         */
    /* calibration in health outcome models to assess diet and health            */
    /* relationships.                                                            */
    /*                                                                           */
    /* The input data for this macro can be prepared using the                   */
    /* STD_COV_BOXCOX24HR_CONDAY_MINAMT macro which is used to standardize       */
    /* continuous covariates, transform and standardize reported nonzero amounts */
    /* (e.g. from 24-hour recalls), and create variables needed for the          */
    /* MULTIVAR_MCMC macro.                                                      */
    /*                                                                           */
    /* The calculations performed by the STD_COV_BOXCOX24HR_CONDAY_MINAMT macro  */
    /* are summarized as follows:                                                */
    /*    1.  Continuous covariates are standardized to have a mean of 0 and a   */
    /*        variance of 1.                                                     */
    /*    2.  For dietary components consumed episodically, consumption-day      */
    /*        indicator variables are created, and the following transformation  */
    /*        approach, used by Zhang et al. (2011a, 2011b), is performed.       */
    /*        First, a suitable Box-Cox transformation is applied to the nonzero */
    /*        values.  Second, the Box-Cox transformed values are centered and   */
    /*        scaled to obtain an approximate normal distribution with a mean of */
    /*        0 and a variance of 2.                                             */
    /*    3.  For dietary components that are consumed nearly every day, as      */
    /*        specified by the user, zero values are replaced with minimum       */
    /*        amount values provided by the user or calculated as half of the    */
    /*        observed nonzero minimum amount.  Then the following               */
    /*        transformation approach, used by Zhang et al. (2011a, 2011b), is   */
    /*        performed for the dietary components that are non-episodically     */
    /*        consumed.  First, a suitable Box-Cox transformation is applied.    */
    /*        Second, the Box-Cox transformed values are centered and scaled to  */
    /*        obtain an approximate normal distribution with a mean of 0 and a   */
    /*        variance of 2.                                                     */
    /*                                                                           */
    /* The MULTIVAR_MCMC macro requires SAS IML version 9.2 or higher.           */
    /*                                                                           */
    /* References:                                                               */
    /*                                                                           */
    /*   Zhang S, Krebs-Smith SM, Midthune D, Perez A, Buckman DW, Kipnis V,     */
    /*   Freedman LS, Dodd KW, Carroll RJ. Fitting a bivariate measurement error */
    /*   model for episodically consumed dietary components. Int J Biostat       */
    /*   2011;7(1):Article 1.                                                    */
    /*                                                                           */
    /*   Zhang S, Midthune D, Guenther PM, Krebs-Smith SM, Kipnis V, Dodd KW,    */
    /*   Buckman DW, Tooze JA, Freedman L, Carroll RJ. A new multivariate        */
    /*   measurement error model with zero-inflated dietary data, and its        */
    /*   application to dietary assessment. Ann Appl Stat 2011 Jun;5(2B):        */
    /*   1456-87.                                                                */
    /*                                                                           */
    /*                                                                           */
    /* The syntax for calling the MULTIVAR_MCMC macro is:                        */
    /*                                                                           */
    /* %multivar_mcmc(data                        =,                             */
    /*                subject                     =,                             */
    /*                weight_var                  =,                             */
    /*                repeat                      =,                             */
    /*                conday_epis_vars            =,                             */
    /*                gst_rec24hr_epis_vars       =,                             */
    /*                gst_rec24hr_daily_vars      =,                             */
    /*                covars_epis_prob            =,                             */
    /*                covars_epis_amt             =,                             */
    /*                covars_daily_amt            =,                             */
    /*                nev_consumers_epis1         =,                             */
    /*                covars_prob_consumer_epis1  =,                             */
    /*                set_seed_mcmc               =,                             */
    /*                set_number_mcmc_iterations  =,                             */
    /*                set_number_burn_iterations  =,                             */
    /*                set_thin                    =,                             */
    /*                prior_sigmau_mean_data      =,                             */
    /*                sigmau_constant             =,                             */
    /*                gen_inverse                 =,                             */
    /*                print                       =,                             */
    /*                titles                      =,                             */
    /*                std_print_store             =,                             */
    /*                notes_print                 =,                             */
    /*                out_lib                     =,                             */
    /*                out_store_label             =,                             */
    /*                out_save_label_max5char     =,                             */
    /*                set_number_saved_out_data   =,                             */
    /*                save_mcmc_u_out_data        =,                             */
    /*                set_number_post_mcmc_u_out  =,                             */
    /*                traceplots_method1_gpath    =,                             */
    /*                traceplots_method2_file_pdf =,                             */
    /*                optional_iml_store_data     =,                             */
    /*                optional_iml_store_names    =                              */
    /*                );                                                         */
    /*                                                                           */
    /*  where                                                                    */
    /*                                                                           */
    /*  "data"                        Specifies an input data set that includes  */
    /*                                one or more observations for each subject. */
    /*                                An adequate number of subjects should have */
    /*                                at least two days of reported consumption  */
    /*                                for episodically consumed dietary          */
    /*                                components, so covariance matrix           */
    /*                                parameters can be estimated.               */
    /*                                                                           */
    /*  "subject"                     Specifies a variable that uniquely         */
    /*                                identifies each subject.                   */
    /*                                                                           */
    /*  "weight_var"                  Specifies a weight variable that is used   */
    /*                                to create a new variable "weight_nw_sumw"  */
    /*                                defined as:                                */
    /*                                    the number of subjects, n, multiplied  */
    /*                                    by the input weight value and divided  */
    /*                                    by the sum of the input weight values. */
    /*                                If "weight_var" is not specified, an       */
    /*                                unweighted analysis is performed where the */
    /*                                value of "weight_nw_sumw" is 1 for each    */
    /*                                subject.                                   */
    /*                                                                           */
    /*  "repeat"                      Specifies a variable that indexes repeated */
    /*                                observations for each subject.             */
    /*                                                                           */
    /*  "conday_epis_vars"            Specifies a list of variables that         */
    /*                                includes a consumption-day indicator       */
    /*                                variable for each episodically consumed    */
    /*                                dietary component.                         */
    /*                                                                           */
    /*  "gst_rec24hr_epis_vars"       Specifies a list of variables that         */
    /*                                includes an appropriately transformed      */
    /*                                24-hour recall variable for each           */
    /*                                episodically consumed dietary component.   */
    /*                                With the priors used in the MCMC, the      */
    /*                                following transformation approach has been */
    /*                                used by Zhang et al. (2011a, 2011b).       */
    /*                                First, a suitable Box-Cox transformation   */
    /*                                is applied to the nonzero values.          */
    /*                                Second, the Box-Cox transformed values are */
    /*                                centered and scaled to obtain an           */
    /*                                approximate normal distribution with a     */
    /*                                mean of 0 and a variance of 2.             */
    /*                                                                           */
    /*  "gst_rec24hr_daily_vars"      Specifies a list of variables that         */
    /*                                includes an appropriately transformed      */
    /*                                24-hour recall variable for each daily     */
    /*                                consumed dietary component.                */
    /*                                With the priors used in the MCMC, the      */
    /*                                following transformation approach has been */
    /*                                used by Zhang et al. (2011a, 2011b).       */
    /*                                First, a suitable Box-Cox transformation   */
    /*                                is applied.                                */
    /*                                Second, the Box-Cox transformed values are */
    /*                                centered and scaled to obtain an           */
    /*                                approximate normal distribution with a     */
    /*                                mean of 0 and a variance of 2.             */
    /*                                                                           */
    /*  "covars_epis_prob"            Specifies a list of covariates for each of */
    /*                                the consumption-day indicator variables    */
    /*                                for the episodically consumed dietary      */
    /*                                components.  If one list of covariates is  */
    /*                                given, these covariates will be used for   */
    /*                                all of the episodically consumed dietary   */
    /*                                components.  More generally, separate      */
    /*                                lists of covariates can be specified for   */
    /*                                each dietary component if each list is     */
    /*                                preceded with the following keyword:       */
    /*                                STARTCOVARLIST.  To include an intercept,  */
    /*                                the list of covariates should include a    */
    /*                                constant such as the variable "constant1"  */
    /*                                that is created by the macro and equals 1  */
    /*                                for each observation of the input data     */
    /*                                set.                                       */
    /*                                                                           */
    /*  "covars_epis_amt"             Specifies a list of covariates for each of */
    /*                                the transformed 24-hour recall variables   */
    /*                                for the episodically consumed dietary      */
    /*                                components.  If one list of covariates is  */
    /*                                given, these covariates will be used for   */
    /*                                all of the episodically consumed dietary   */
    /*                                components.  More generally, separate      */
    /*                                lists of covariates can be specified for   */
    /*                                each dietary component if each list is     */
    /*                                preceded with the following keyword:       */
    /*                                STARTCOVARLIST.  To include an intercept,  */
    /*                                the list of covariates should include a    */
    /*                                constant such as the variable "constant1"  */
    /*                                that is created by the macro and equals 1  */
    /*                                for each observation of the input data     */
    /*                                set.                                       */
    /*                                                                           */
    /*  "covars_daily_amt"            Specifies a list of covariates for each of */
    /*                                the transformed 24-hour recall variables   */
    /*                                for the daily consumed dietary components. */
    /*                                If one list of covariates is given, these  */
    /*                                covariates will be used for all of the     */
    /*                                daily consumed dietary components.  More   */
    /*                                generally, separate lists of covariates    */
    /*                                can be specified for each dietary          */
    /*                                component if each list is preceded with    */
    /*                                the following keyword:                     */
    /*                                STARTCOVARLIST.  To include an intercept,  */
    /*                                the list of covariates should include a    */
    /*                                constant such as the variable "constant1"  */
    /*                                that is created by the macro and equals 1  */
    /*                                for each observation of the input data     */
    /*                                set.                                       */
    /*                                                                           */
    /*  "nev_consumers_epis1"         If "nev_consumers_epis1=y" or              */
    /*                                "nev_consumers_epis1=Y" then the           */
    /*                                calculations will allow for never-         */
    /*                                consumers for the first dietary component  */
    /*                                listed in the "conday_epis_vars" and the   */
    /*                                "gst_rec24hr_epis_vars" variable lists.    */
    /*                                The default value is "n".                  */
    /*                                                                           */
    /*  "covars_prob_consumer_epis1"  Specifies a list of covariates used to     */
    /*                                model the probability of being a consumer  */
    /*                                of the first dietary component listed in   */
    /*                                the "conday_epis_vars" and the             */
    /*                                "gst_rec24hr_epis_vars" variable lists.    */
    /*                                To include an intercept, the list of       */
    /*                                covariates should include a constant such  */
    /*                                as the variable "constant1" that is        */
    /*                                created by the macro and equals 1 for each */
    /*                                observation of the input data set.  This   */
    /*                                list of covariates is only used if         */
    /*                                "nev_consumers_epis1=y" or                 */
    /*                                "nev_consumers_epis1=Y".                   */
    /*                                                                           */
    /*  "set_seed_mcmc"               Specifies a seed for random number         */
    /*                                generation in the MCMC algorithm.  If      */
    /*                                "set_seed_mcmc" is not specified, the SAS  */
    /*                                IML procedure will generate an initial     */
    /*                                seed value from the system clock.          */
    /*                                                                           */
    /*  "set_number_mcmc_iterations"  Specifies the total number of iterations   */
    /*                                for the MCMC algorithm including the       */
    /*                                burn-in.  The default value is 12000.      */
    /*                                                                           */
    /*  "set_number_burn_iterations"  Specifies the burn-in value for the MCMC   */
    /*                                algorithm.  The default value is 2000.     */
    /*                                                                           */
    /*  "set_thin"                    Specifies the thinning value for the MCMC  */
    /*                                algorithm.  The default value is 25.       */
    /*                                                                           */
    /*  "prior_sigmau_mean_data"      Specifies an optional SAS data set used to */
    /*                                define the mean of the prior distribution  */
    /*                                for the sigmau covariance matrix.  The     */
    /*                                values from this SAS data set are read     */
    /*                                into a SAS IML procedure matrix,           */
    /*                                prior_sigmau_mean, so the optional input   */
    /*                                SAS data set can only include the values   */
    /*                                needed to define this matrix.  Also, these */
    /*                                values should be appropriately arranged    */
    /*                                with each data set variable corresponding  */
    /*                                to a column of this matrix.  Advanced      */
    /*                                users can use this macro parameter in      */
    /*                                conjunction with the "sigmau_constant"     */
    /*                                macro parameter.  The description for the  */
    /*                                "out_store_label" macro parameter and      */
    /*                                references provide notational details.     */
    /*                                                                           */
    /*  "sigmau_constant"             If "sigmau_constant=y" or                  */
    /*                                "sigmau_constant=Y" then the sigmau        */
    /*                                covariance matrix is not updated during    */
    /*                                the MCMC iterative calculations.  Advanced */
    /*                                users can use this macro parameter in      */
    /*                                conjunction with the                       */
    /*                                "prior_sigmau_mean_data" macro parameter.  */
    /*                                The description for the "out_store_label"  */
    /*                                macro parameter and the references provide */
    /*                                notational details.  The default value is  */
    /*                                "n".                                       */
    /*                                                                           */
    /*  "gen_inverse"                 If "gen_inverse=n" or "gen_inverse=N" then */
    /*                                the SAS IML "inv" function is used to      */
    /*                                compute the inverse of matrices as needed  */
    /*                                throughout the algorithm.  This            */
    /*                                specification may save computational time, */
    /*                                but the user may encounter numerical       */
    /*                                problems and an error message regarding a  */
    /*                                matrix that should be non-singular.  The   */
    /*                                default value is "y".  By default, the SAS */
    /*                                IML "ginv" function is used to compute the */
    /*                                Moore-Penrose generalized inverse of       */
    /*                                matrices as needed throughout the          */
    /*                                algorithm.                                 */
    /*                                                                           */
    /*  "print"                       If "print=n" or "print=N" then macro       */
    /*                                results are not printed.  The default      */
    /*                                value is "y".                              */
    /*                                                                           */
    /*  "titles"                      Specifies the number of title lines to be  */
    /*                                reserved for the user's titles.  One       */
    /*                                additional title line is used by the       */
    /*                                macro.  The default value is "0".          */
    /*                                                                           */
    /*  "std_print_store"             If "std_print_store=y" or                  */
    /*                                "std_print_store=Y" then sample standard   */
    /*                                deviations are printed and stored in the   */
    /*                                IML storage catalog.  The default value is */
    /*                                "n" because some analyses require more     */
    /*                                sophisticated methods of variance          */
    /*                                estimation.  For example, replication      */
    /*                                methods such as balanced repeated          */
    /*                                replication (BRR) are used for analysis of */
    /*                                data from a survey with complex sampling.  */
    /*                                                                           */
    /*  "notes_print"                 If "notes_print=n" or "notes_print=N" then */
    /*                                notes are not printed to the SAS log.  The */
    /*                                default value is "y".                      */
    /*                                                                           */
    /*  "out_lib"                     Specifies a SAS library that the macro     */
    /*                                uses when saving the IML storage catalog   */
    /*                                file and optional SAS data sets.           */
    /*                                The description for the "out_store_label"  */
    /*                                macro parameter includes additional        */
    /*                                details regarding the IML storage catalog. */
    /*                                The following SAS output data sets         */
    /*                                    "multivar_mcmc_samples_out"            */
    /*                                    "multivar_mcmc_samples_u_out"          */
    /*                                    "multivar_post_mcmc_u_out"             */
    /*                                are requested using the                    */
    /*                                "set_number_saved_out_data" and            */
    /*                                "save_mcmc_u_out_data" and                 */
    /*                                "set_number_post_mcmc_u_out" macro         */
    /*                                parameters and are saved in the library    */
    /*                                specified by the "out_lib" macro           */
    /*                                parameter.                                 */
    /*                                                                           */
    /*  "out_store_label"             Specifies a label for the SAS IML storage  */
    /*                                catalog.  The catalog name is defined by   */
    /*                                adding the prefix "iml" to the user        */
    /*                                supplied label.  The catalog is stored in  */
    /*                                the library specified by the "out_lib"     */
    /*                                macro parameter.  After deleting the       */
    /*                                current contents of the catalog, the       */
    /*                                following matrices (including vectors and  */
    /*                                scalars) are stored.                       */
    /*                                "weight_nw_sumw"                           */
    /*                                    The number of subjects, n, multiplied  */
    /*                                    by the input weight value and divided  */
    /*                                    by the sum of the input weight values. */
    /*                                "num_epis_diet_comp"                       */
    /*                                    The number of episodically consumed    */
    /*                                    dietary components.                    */
    /*                                "num_daily_diet_comp"                      */
    /*                                    The number of daily consumed dietary   */
    /*                                    components.                            */
    /*                                "num_rows_covmat"                          */
    /*                                    The number of rows (or columns) of the */
    /*                                    LxL covariance matrices, sigmau_mean   */
    /*                                    and sigmae_mean, where                 */
    /*                                    L = (2 * num_epis_diet_comp)           */
    /*                                         + num_daily_diet_comp.            */
    /*                                "beta1_mean", ..., "betaL_mean"            */
    /*                                "sigmau_mean"                              */
    /*                                "sigmae_mean"                              */
    /*                                    Sample means calculated using          */
    /*                                    generated values from MCMC iterations  */
    /*                                    selected according to the              */
    /*                                    "set_number_burn_iterations" and       */
    /*                                    "set_thin" macro parameters.  The rows */
    /*                                    of the LxL covariance matrices,        */
    /*                                    sigmau_mean and sigmae_mean, are       */
    /*                                    arranged, so the episodically consumed */
    /*                                    dietary components precede the daily   */
    /*                                    consumed dietary components.  For      */
    /*                                    example, 3 episodically consumed       */
    /*                                    dietary components and 2 daily         */
    /*                                    consumed dietary components yield      */
    /*                                    (2*3 + 2 = 8 = L) rows sorted as:      */
    /*                                    1. consumption probability (epis 1),   */
    /*                                    2. amount consumed         (epis 1),   */
    /*                                    3. consumption probability (epis 2),   */
    /*                                    4. amount consumed         (epis 2),   */
    /*                                    5. consumption probability (epis 3),   */
    /*                                    6. amount consumed         (epis 3),   */
    /*                                    7. amount consumed         (daily 1),  */
    /*                                    8. amount consumed         (daily 2),  */
    /*                                "beta1_std", ..., "betaL_std"              */
    /*                                "sigmau_std"                               */
    /*                                "sigmae_std"                               */
    /*                                    Sample standard deviations calculated  */
    /*                                    using generated values from MCMC       */
    /*                                    iterations selected according to the   */
    /*                                    "set_number_burn_iterations" and       */
    /*                                    "set_thin" macro parameters.  These    */
    /*                                    sample standard deviations are stored  */
    /*                                    when "std_print_store=y" or            */
    /*                                    "std_print_store=Y".                   */
    /*                                                                           */
    /*  "out_save_label_max5char"     Specifies an optional label to append to   */
    /*                                the end of the following SAS output data   */
    /*                                set names:                                 */
    /*                                    "multivar_mcmc_samples_out"            */
    /*                                    "multivar_mcmc_samples_u_out"          */
    /*                                    "multivar_post_mcmc_u_out".            */
    /*                                These SAS output data sets are requested   */
    /*                                using the                                  */
    /*                                "set_number_saved_out_data" and            */
    /*                                "save_mcmc_u_out_data" and                 */
    /*                                "set_number_post_mcmc_u_out" macro         */
    /*                                parameters and are saved in the library    */
    /*                                specified by the "out_lib" macro           */
    /*                                parameter.  Only the first 5 characters of */
    /*                                the label are used.                        */
    /*                                                                           */
    /*  "set_number_saved_out_data"   Specifies the number of iterations to save */
    /*                                in a SAS output data set                   */
    /*                                "multivar_mcmc_samples_out" that is        */
    /*                                saved in the library specified by the      */
    /*                                "out_lib" macro parameter.  The data set   */
    /*                                includes MCMC samples for the Beta,        */
    /*                                Sigma_u, and Sigma_e parameters with       */
    /*                                values selected from among the iterations  */
    /*                                specified according to the                 */
    /*                                "set_number_burn_iterations" and           */
    /*                                "set_thin" macro parameters.  The output   */
    /*                                data set name can be extended to include   */
    /*                                an optional label specified using the      */
    /*                                "out_save_label_max5char" macro parameter. */
    /*                                The default value is "0".                  */
    /*                                                                           */
    /*  "save_mcmc_u_out_data"        If "save_mcmc_u_out_data=y" or             */
    /*                                "save_mcmc_u_out_data=Y" then the MCMC     */
    /*                                samples for the u matrix are saved in a    */
    /*                                SAS output data set                        */
    /*                                "multivar_mcmc_samples_u_out" in the       */
    /*                                library specified by the "out_lib" macro   */
    /*                                parameter.  The variables include:         */
    /*                                    - "Iteration" which identifies the     */
    /*                                      MCMC sample,                         */
    /*                                    - the variable specified by the        */
    /*                                      "subject" macro parameter,           */
    /*                                    - u_col1, u_col2, ..., u_colL where L  */
    /*                                      is the number of columns of the u    */
    /*                                      matrix.  L is defined in the         */
    /*                                      description of the "out_store_label" */
    /*                                      macro parameter.                     */
    /*                                The data set includes MCMC samples         */
    /*                                selected from among the iterations         */
    /*                                specified according to the                 */
    /*                                "set_number_burn_iterations" and           */
    /*                                "set_thin" macro parameters.  The output   */
    /*                                data set name can be extended to include   */
    /*                                an optional label specified using the      */
    /*                                "out_save_label_max5char" macro parameter. */
    /*                                The default value is "n".                  */
    /*                                                                           */
    /*  "set_number_post_mcmc_u_out"  Specifies the number of post MCMC          */
    /*                                iterations and the number of post MCMC     */
    /*                                samples for the u matrix that are saved in */
    /*                                a SAS output data set                      */
    /*                                "multivar_post_mcmc_u_out" in the          */
    /*                                library specified by the "out_lib" macro   */
    /*                                parameter.  The variables include:         */
    /*                                    - "Iteration" which identifies the     */
    /*                                      post MCMC sample,                    */
    /*                                    - the variable specified by the        */
    /*                                      "subject" macro parameter,           */
    /*                                    - u_col1, u_col2, ..., u_colL where L  */
    /*                                      is the number of columns of the u    */
    /*                                      matrix.  L is defined in the         */
    /*                                      description of the "out_store_label" */
    /*                                      macro parameter.                     */
    /*                                The output data set name can be extended   */
    /*                                to include an optional label specified     */
    /*                                using the "out_save_label_max5char" macro  */
    /*                                parameter.  The default value is "0".      */
    /*                                                                           */
    /*  "traceplots_method1_gpath"    Specifies a valid SAS fileref to indicate  */
    /*                                a folder used for storing the MCMC trace   */
    /*                                plot files.  The fileref is used by the    */
    /*                                "gpath" option of the "ods listing"        */
    /*                                statement in the SAS Output Delivery       */
    /*                                System (ODS).  Each trace plot file        */
    /*                                includes up to 4 trace plots               */
    /*                                (i.e. 4 panels).  Example file names       */
    /*                                include:                                   */
    /*                                    "Beta_Trace_Plot_Panels4_Image.png"    */
    /*                                    "Sigmau_Trace_Plot_Panels4_Image.png"  */
    /*                                    "Sigmae_Trace_Plot_Panels4_Image.png"  */
    /*                                and SAS ODS will add an image index number */
    /*                                to ensure that each file produced has a    */
    /*                                unique name.  MCMC trace plots are         */
    /*                                produced for all of the Beta, Sigma_u, and */
    /*                                Sigma_e parameters.  The plots include     */
    /*                                MCMC samples selected according to the     */
    /*                                "set_number_burn_iterations" and           */
    /*                                "set_thin" macro parameters.  If no value  */
    /*                                is specified, these trace plot files are   */
    /*                                not produced.  Trace plots can be produced */
    /*                                using another method as described for the  */
    /*                                "traceplots_method2_file_pdf" macro        */
    /*                                parameter.                                 */
    /*                                                                           */
    /*  "traceplots_method2_file_pdf" Specifies a file name with a "pdf"         */
    /*                                extension, such as                         */
    /*                                "example1.traceplot.pdf".  This pdf file   */
    /*                                will include the MCMC trace plots for all  */
    /*                                of the Beta, Sigma_u, and Sigma_e          */
    /*                                parameters.  The plots include MCMC        */
    /*                                samples selected according to the          */
    /*                                "set_number_burn_iterations" and           */
    /*                                "set_thin" macro parameters.  If no value  */
    /*                                is specified, this trace plot file is not  */
    /*                                produced.  Trace plots can be produced     */
    /*                                using another method as described for the  */
    /*                                "traceplots_method1_gpath" macro           */
    /*                                parameter.  The pdf file produced by trace */
    /*                                plot method 2 tends to have a much larger  */
    /*                                file size when compared to the set of      */
    /*                                files produced by trace plot method 1.     */
    /*                                                                           */
    /*  "optional_iml_store_data"     Specifies an optional SAS data set that    */
    /*                                includes useful variables to be stored in  */
    /*                                the IML storage catalog.  The names of     */
    /*                                these variables must be specified in the   */
    /*                                "optional_iml_store_names" macro           */
    /*                                parameter.                                 */
    /*                                                                           */
    /*  "optional_iml_store_names"    Specifies a list of optional matrices      */
    /*                                (including vectors and scalars) to be      */
    /*                                stored in the IML storage catalog.  This   */
    /*                                macro parameter can be used in conjunction */
    /*                                with the "optional_iml_store_data" macro   */
    /*                                parameter, and it can be used by advanced  */
    /*                                users that want to save values that are    */
    /*                                defined in the SAS IML procedure within    */
    /*                                the macro.                                 */
    /*                                                                           */
    /*****************************************************************************/
    /*****************************************************************************/
    /*****************************************************************************/

    *******************************************************;
    *** Only create trace plots for every 20 replicates ***;
    *******************************************************;

    %let create_traceplot_every = 20;
    %if %sysfunc(mod(&replicnum, &create_traceplot_every)) = 0 %then %let traceplot = epidemiology_mcmc_nc_main24hr_traceplot.rep&replicnum..pdf;
    %else %let traceplot = %str();


    ******************************************************************************************************;
    *** The num_subsets_monte_carlo_distrib and subset_number_post_mcmc_u_out macro variables are      ***;
    *** used to illustrate an approach that requires less memory.  The num_subsets_monte_carlo_distrib ***;
    *** macro variable is used to set the number of subsets.  The subset_number_post_mcmc_u_out macro  ***;
    *** variable is used to set the number of iterations per subset.                                   ***;
    ******************************************************************************************************;

    %let num_subsets_monte_carlo_distrib = 4;
    %let subset_number_post_mcmc_u_out = 250;
    %let number_post_mcmc_u_out = %eval(&num_subsets_monte_carlo_distrib * &subset_number_post_mcmc_u_out);


    %let loopseed_mcmc = %eval(85600495 + &replicnum * 10);

    ****************************************************;
    *** Fit the multivariate measurement error model ***;
    ****************************************************;

    %multivar_mcmc(data                        = stdcov_stdbc24hr_conday_out,
                   subject                     = replicaterowid,
                   weight_var                  = ,
                   repeat                      = repeat,
                   conday_epis_vars            = conday_r_alcohol_g conday_r_redmeat_g conday_r_whlgrain_ozeq,
                   gst_rec24hr_epis_vars       = stdbc_r_alcohol_g stdbc_r_redmeat_g stdbc_r_whlgrain_ozeq,
                   gst_rec24hr_daily_vars      = stdbc_r_sodium_mg stdbc_r_potassium_mg stdbc_r_protein_g stdbc_r_energy_kcal,
                   covars_epis_prob            = constant1 ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                                                 std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg
                                                 std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   covars_epis_amt             = constant1 ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                                                 std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg
                                                 std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   covars_daily_amt            = constant1 ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                                                 std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg
                                                 std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                   nev_consumers_epis1         = y,
                   covars_prob_consumer_epis1  = constant1 ffq_eqzero_alcohol_g,
                   set_seed_mcmc               = &loopseed_mcmc,
                   set_number_mcmc_iterations  = 2500,
                   set_number_burn_iterations  = 500,
                   set_thin                    = 10,
                   prior_sigmau_mean_data      = ,
                   sigmau_constant             = ,
                   gen_inverse                 = ,
                   print                       = y,
                   titles                      = 4,
                   std_print_store             = y,
                   notes_print                 = y,
                   out_lib                     = outlib,
                   out_store_label             = epi_mcmc_nc_main24hr&replicnum,
                   out_save_label_max5char     = _r&replicnum,
                   set_number_saved_out_data   = ,
                   save_mcmc_u_out_data        = ,
                   set_number_post_mcmc_u_out  = &number_post_mcmc_u_out,
                   traceplots_method1_gpath    = ,
                   traceplots_method2_file_pdf = &traceplot,
                   optional_iml_store_data     = backtran_out,
                   optional_iml_store_names    = tran_paramindex tran_lambda tran_center tran_scale minamount
                   );



    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* Description of the MULTIVAR_DISTRIB macro                                 */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /*                                                                           */
    /* The MULTIVAR_DISTRIB macro uses parameter estimates from the              */
    /* MULTIVAR_MCMC macro to generate a multivariate Monte Carlo distribution   */
    /* of the usual intakes for the dietary components specified in the          */
    /* multivariate measurement error model fit using the MULTIVAR_MCMC macro.   */
    /* The MULTIVAR_DISTRIB macro can also use values that were stored using     */
    /* the macro parameter "optional_iml_store_names" of the MULTIVAR_MCMC       */
    /* macro.  The MULTIVAR_DISTRIB macro allows specification of 1 or 2 lists   */
    /* of covariates.  If 2 lists of covariates are specified, the usual intake  */
    /* is calculated as the weighted average of the usual intake calculated      */
    /* using covariate list 1 and the usual intake calculated using covariate    */
    /* list 2.  This weighted average is calculated using the values specified   */
    /* for the "set_value_for_weight_cov_list1" and                              */
    /* "set_value_for_weight_cov_list2" macro parameters.                        */
    /*                                                                           */
    /* References:                                                               */
    /*                                                                           */
    /*   Zhang S, Krebs-Smith SM, Midthune D, Perez A, Buckman DW, Kipnis V,     */
    /*   Freedman LS, Dodd KW, Carroll RJ. Fitting a bivariate measurement error */
    /*   model for episodically consumed dietary components. Int J Biostat       */
    /*   2011;7(1):Article 1.                                                    */
    /*                                                                           */
    /*   Zhang S, Midthune D, Guenther PM, Krebs-Smith SM, Kipnis V, Dodd KW,    */
    /*   Buckman DW, Tooze JA, Freedman L, Carroll RJ. A new multivariate        */
    /*   measurement error model with zero-inflated dietary data, and its        */
    /*   application to dietary assessment. Ann Appl Stat 2011 Jun;5(2B):        */
    /*   1456-87.                                                                */
    /*                                                                           */
    /*                                                                           */
    /* The syntax for calling the MULTIVAR_DISTRIB macro is:                     */
    /*                                                                           */
    /* %multivar_distrib(multivar_mcmc_out_lib           =,                      */
    /*                   multivar_mcmc_out_store_label   =,                      */
    /*                   t_weightavg_covariates_list1    =,                      */
    /*                   t_weightavg_covariates_list2    =,                      */
    /*                   set_value_for_weight_cov_list1  =,                      */
    /*                   set_value_for_weight_cov_list2  =,                      */
    /*                   nev_consumers_epis1             =,                      */
    /*                   covars_prob_consumer_epis1      =,                      */
    /*                   optional_input_data             =,                      */
    /*                   optional_input_data_var_list    =,                      */
    /*                   optional_input_mcmc_u_out_data  =,                      */
    /*                   additional_output_var_list      =,                      */
    /*                   additional_output_subject_var   =,                      */
    /*                   output_mcmc_weight_var          =,                      */
    /*                   set_seed_distrib                =,                      */
    /*                   set_number_monte_carlo_rand_obs =,                      */
    /*                   print                           =                       */
    /*                   );                                                      */
    /*                                                                           */
    /*  where                                                                    */
    /*                                                                           */
    /*  "multivar_mcmc_out_lib"            Specifies the same SAS library that   */
    /*                                     was specified for the "out_lib" macro */
    /*                                     parameter of the MULTIVAR_MCMC macro. */
    /*                                     The MULTIVAR_DISTRIB macro will use   */
    /*                                     the parameter estimates and can use   */
    /*                                     covariates that were stored by the    */
    /*                                     MULTIVAR_MCMC macro.  Covariates can  */
    /*                                     be stored using the                   */
    /*                                     "optional_iml_store_names" macro      */
    /*                                     parameter of the MULTIVAR_MCMC macro. */
    /*                                                                           */
    /*  "multivar_mcmc_out_store_label"    Specifies the same label that was     */
    /*                                     specified for the "out_store_label"   */
    /*                                     macro parameter of the MULTIVAR_MCMC  */
    /*                                     macro.  The MULTIVAR_DISTRIB macro    */
    /*                                     will use the parameter estimates and  */
    /*                                     can use covariates that were stored   */
    /*                                     by the MULTIVAR_MCMC macro.           */
    /*                                     Covariates can be stored using the    */
    /*                                     "optional_iml_store_names" macro      */
    /*                                     parameter of the MULTIVAR_MCMC macro. */
    /*                                                                           */
    /*  "t_weightavg_covariates_list1"     Specifies a list of covariates that   */
    /*                                     correspond to the beta parameter      */
    /*                                     estimates from the multivariate       */
    /*                                     measurement error model fit using the */
    /*                                     MULTIVAR_MCMC macro.  The             */
    /*                                     calculations in the MULTIVAR_DISTRIB  */
    /*                                     macro are based on the assumption     */
    /*                                     that the multivariate measurement     */
    /*                                     error model was fit using the same    */
    /*                                     set of covariates for each part of    */
    /*                                     the model.  The MULTIVAR_DISTRIB      */
    /*                                     macro creates "constant0" which can   */
    /*                                     be used as a variable that equals 0   */
    /*                                     for each observation, and the macro   */
    /*                                     creates "constant1" which can be used */
    /*                                     as a variable that equals 1 for each  */
    /*                                     observation.  "Constant0" and         */
    /*                                     "constant1" can be used when          */
    /*                                     specifying covariates of interest.    */
    /*                                     If the optional macro parameter       */
    /*                                     "t_weightavg_covariates_list2" is     */
    /*                                     used to specify a second list of      */
    /*                                     covariates, then usual intake is      */
    /*                                     calculated as the weighted average of */
    /*                                     the usual intake calculated using     */
    /*                                     covariate list 1 and the usual intake */
    /*                                     calculated using covariate list 2.    */
    /*                                     This weighted average is calculated   */
    /*                                     using the values specified for the    */
    /*                                     "set_value_for_weight_cov_list1" and  */
    /*                                     "set_value_for_weight_cov_list2"      */
    /*                                     macro parameters.  If the second list */
    /*                                     of covariates is not specified, then  */
    /*                                     the Monte Carlo distribution of usual */
    /*                                     intakes is only based on the          */
    /*                                     covariates specified for the          */
    /*                                     "t_weightavg_covariates_list1" macro  */
    /*                                     parameter.                            */
    /*                                                                           */
    /*  "t_weightavg_covariates_list2"     Specifies an optional second list of  */
    /*                                     covariates.  If this list of          */
    /*                                     covariates is specified, the usual    */
    /*                                     intake is calculated as the weighted  */
    /*                                     average of the usual intake           */
    /*                                     calculated using covariate list 1 and */
    /*                                     the usual intake calculated using     */
    /*                                     covariate list 2.  This weighted      */
    /*                                     average is calculated using the       */
    /*                                     values specified for the              */
    /*                                     "set_value_for_weight_cov_list1" and  */
    /*                                     "set_value_for_weight_cov_list2"      */
    /*                                     macro parameters.                     */
    /*                                                                           */
    /*  "set_value_for_weight_cov_list1"   Specifies a value that is used to     */
    /*                                     calculate a weighted average when 2   */
    /*                                     covariate lists are provided as       */
    /*                                     described for the                     */
    /*                                     "t_weightavg_covariates_list1" and    */
    /*                                     "t_weightavg_covariates_list2" macro  */
    /*                                     parameters.  The value specified is   */
    /*                                     used as the numerator of the weight.  */
    /*                                     The denominator of the weight is      */
    /*                                     defined as the sum of the values from */
    /*                                     the "set_value_for_weight_cov_list1"  */
    /*                                     and "set_value_for_weight_cov_list2"  */
    /*                                     macro parameters.  The default value  */
    /*                                     is 1 which corresponds to the         */
    /*                                     situation where the                   */
    /*                                     "t_weightavg_covariates_list1" macro  */
    /*                                     parameter is used to provide          */
    /*                                     covariate list 1 and the              */
    /*                                     "t_weightavg_covariates_list2" macro  */
    /*                                     parameter is not used.                */
    /*                                                                           */
    /*  "set_value_for_weight_cov_list2"   Specifies a value that is used to     */
    /*                                     calculate a weighted average when 2   */
    /*                                     covariate lists are provided as       */
    /*                                     described for the                     */
    /*                                     "t_weightavg_covariates_list1" and    */
    /*                                     "t_weightavg_covariates_list2" macro  */
    /*                                     parameters.  The value specified is   */
    /*                                     used as the numerator of the weight.  */
    /*                                     The denominator of the weight is      */
    /*                                     defined as the sum of the values from */
    /*                                     the "set_value_for_weight_cov_list1"  */
    /*                                     and "set_value_for_weight_cov_list2"  */
    /*                                     macro parameters.  The default value  */
    /*                                     is 0 which corresponds to the         */
    /*                                     situation where the                   */
    /*                                     "t_weightavg_covariates_list2" macro  */
    /*                                     parameter is not used.                */
    /*                                                                           */
    /*  "nev_consumers_epis1"              Specifies the same setting that was   */
    /*                                     specified for the                     */
    /*                                     "nev_consumers_epis1" macro parameter */
    /*                                     of the MULTIVAR_MCMC macro.  The      */
    /*                                     default value is "n".                 */
    /*                                                                           */
    /*  "covars_prob_consumer_epis1"       Specifies the same list of covariates */
    /*                                     that was specified for the            */
    /*                                     "covars_prob_consumer_epis1"          */
    /*                                     macro parameter of the MULTIVAR_MCMC  */
    /*                                     macro.                                */
    /*                                                                           */
    /*  "optional_input_data"              Specifies an optional SAS data set    */
    /*                                     that provides the additional input    */
    /*                                     variables listed in the               */
    /*                                     "optional_input_data_var_list" macro  */
    /*                                     parameter.  If this data set is       */
    /*                                     specified and the                     */
    /*                                     "additional_output_subject_var" macro */
    /*                                     parameter specifies a subject         */
    /*                                     identification variable, then the     */
    /*                                     subject identification variable is    */
    /*                                     included in the output Monte Carlo    */
    /*                                     data set.  When an optional input     */
    /*                                     data set is not specified, the Monte  */
    /*                                     Carlo "mc_t_distrib_out" data set     */
    /*                                     includes the "weight_nw_sumw"         */
    /*                                     variable created and stored by the    */
    /*                                     MULTIVAR_MCMC macro (see the          */
    /*                                     MULTIVAR_MCMC "weight_var" macro      */
    /*                                     parameter for details).  When an      */
    /*                                     optional input data set is specified, */
    /*                                     an advanced user can use the          */
    /*                                     "output_mcmc_weight_var" macro        */
    /*                                     parameter to include the              */
    /*                                     "weight_nw_sumw" variable in the      */
    /*                                     "mc_t_distrib_out" data set if        */
    /*                                     appropriate.  When an optional input  */
    /*                                     data set and both the                 */
    /*                                     "additional_output_var_list" and      */
    /*                                     "additional_output_subject_var" macro */
    /*                                     parameters are specified, an advanced */
    /*                                     user can include these additional     */
    /*                                     variables in the "mc_t_distrib_out"   */
    /*                                     data set.                             */
    /*                                                                           */
    /*  "optional_input_data_var_list"     Specifies optional input variables.   */
    /*                                     For details, see the description for  */
    /*                                     the "optional_input_data" macro       */
    /*                                     parameter.  This variable list should */
    /*                                     not include names of variables stored */
    /*                                     by the MULTIVAR_MCMC macro since the  */
    /*                                     stored values are brought into the    */
    /*                                     macro as described for the            */
    /*                                     "multivar_mcmc_out_lib" and           */
    /*                                     "multivar_mcmc_out_store_label" macro */
    /*                                     parameters.                           */
    /*                                                                           */
    /*  "optional_input_mcmc_u_out_data"   Specifies an optional SAS data set    */
    /*                                     that provides u matrices saved from   */
    /*                                     iterations of the MULTIVAR_MCMC       */
    /*                                     macro.  Considering the MULTIVAR_MCMC */
    /*                                     "set_number_post_mcmc_u_out" macro    */
    /*                                     parameter, the variables required for */
    /*                                     this input data set include:          */
    /*                                     - "Iteration" which indexes the first */
    /*                                       u matrix through last u matrix in   */
    /*                                       the data set,                       */
    /*                                     - the variable specified by the       */
    /*                                       "subject" macro parameter from the  */
    /*                                       MULTIVAR_MCMC macro,                */
    /*                                     - u_col1, u_col2, ..., u_colL where L */
    /*                                       is the number of columns of the u   */
    /*                                       matrix.  L is defined in the        */
    /*                                       description of the                  */
    /*                                       "out_store_label" macro parameter   */
    /*                                       of the MULTIVAR_MCMC macro.         */
    /*                                     The MULTIVAR_DISTRIB macro will use   */
    /*                                     the first (i.e. iteration 1) u matrix */
    /*                                     under the assumption that the rows of */
    /*                                     this matrix are still sorted in the   */
    /*                                     same order as when the MULTIVAR_MCMC  */
    /*                                     macro generated this u matrix, and    */
    /*                                     the MULTIVAR_DISTRIB macro will       */
    /*                                     generate and append additional rows   */
    /*                                     to this matrix as needed to obtain    */
    /*                                     the dimensions as specified and       */
    /*                                     implied by other user input to this   */
    /*                                     macro.  The MULTIVAR_DISTRIB macro    */
    /*                                     will use this same approach for each  */
    /*                                     subsequent u matrix and iteration.    */
    /*                                                                           */
    /*  "additional_output_var_list"       Specifies additional output           */
    /*                                     variables.  For details, see the      */
    /*                                     description for the                   */
    /*                                     "optional_input_data" macro           */
    /*                                     parameter.                            */
    /*                                                                           */
    /*  "additional_output_subject_var"    Specifies a subject identification    */
    /*                                     variable to be included in the output */
    /*                                     Monte Carlo data set.  For details,   */
    /*                                     see the description for the           */
    /*                                     "optional_input_data" macro           */
    /*                                     parameter.                            */
    /*                                                                           */
    /*  "output_mcmc_weight_var"           When an optional input data set is    */
    /*                                     specified, an advanced user can       */
    /*                                     specify "output_mcmc_weight_var=y" or */
    /*                                     "output_mcmc_weight_var=Y" to include */
    /*                                     the "weight_nw_sumw" variable in the  */
    /*                                     "mc_t_distrib_out" data set if        */
    /*                                     appropriate.  For details, see the    */
    /*                                     description for the                   */
    /*                                     "optional_input_data" macro           */
    /*                                     parameter.  The default value is "n". */
    /*                                                                           */
    /*  "set_seed_distrib"                 Specifies a seed for random number    */
    /*                                     generation for the Monte Carlo method */
    /*                                     used in the macro.  If                */
    /*                                     "set_seed_distrib" is not specified,  */
    /*                                     the SAS IML procedure will generate   */
    /*                                     an initial seed value from the system */
    /*                                     clock.                                */
    /*                                                                           */
    /*  "set_number_monte_carlo_rand_obs"  Specifies the number of observations  */
    /*                                     of the random effects vector to       */
    /*                                     generate for an individual.  If an    */
    /*                                     "optional_input_mcmc_u_out_data" SAS  */
    /*                                     data set is specified, then the       */
    /*                                     "set_number_monte_carlo_rand_obs"     */
    /*                                     macro parameter is not used, and the  */
    /*                                     number of observations of the random  */
    /*                                     effects vector is determined by the   */
    /*                                     maximum value of the "Iteration"      */
    /*                                     index variable.  Additional details   */
    /*                                     are provided in the description of    */
    /*                                     the "optional_input_mcmc_u_out_data"  */
    /*                                     macro parameter.  If an               */
    /*                                     "optional_input_mcmc_u_out_data" SAS  */
    /*                                     data set is not specified, the        */
    /*                                     number of observations in the Monte   */
    /*                                     Carlo data set is the product of the  */
    /*                                     "number of individuals" and the       */
    /*                                     specified value for the               */
    /*                                     "set_number_monte_carlo_rand_obs"     */
    /*                                     macro parameter.  The default value   */
    /*                                     is 500.                               */
    /*                                                                           */
    /*  "print"                            If "print=n" or "print=N" then macro  */
    /*                                     information is not printed.           */
    /*                                     Otherwise the macro prints            */
    /*                                     information regarding the IML storage */
    /*                                     catalog used and prints some stored   */
    /*                                     and specified values.  The default    */
    /*                                     value is "y".                         */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /* Macro Output Data Set:  mc_t_distrib_out                                  */
    /*                                                                           */
    /* The MULTIVAR_DISTRIB macro produces a SAS data set "mc_t_distrib_out"     */
    /* that includes usual intake variables named:                               */
    /*          mc_t1   mc_t2   ...   mc_tp                                      */
    /* where p = num_epis_diet_comp + num_daily_diet_comp, and the dietary       */
    /* component order (1, 2, ..., p) is equivalent to the order used when the   */
    /* multivariate measurement error model was fit using the MULTIVAR_MCMC      */
    /* macro.  When the "t_weightavg_covariates_list2" macro parameter is not    */
    /* used, the output data set "mc_t_distrib_out" includes the following       */
    /* variables used to calculate the usual intake variables:                   */
    /*          mc_prob1   mc_prob2   ...   mc_probq                             */
    /*          mc_amount1   mc_amount2   ...   mc_amountq                       */
    /* where q = num_epis_diet_comp, and these variables correspond to the       */
    /* probability to consume and the intake amount on consumption days.         */
    /*                                                                           */
    /* When the "optional_input_data" macro parameter is not used to specify an  */
    /* optional input data set, the "mc_t_distrib_out" data set also             */
    /* includes the variable:                                                    */
    /*          weight_nw_sumw                                                   */
    /* created and stored by the MULTIVAR_MCMC macro (see the MULTIVAR_MCMC      */
    /* "weight_var" macro parameter for details).                                */
    /*                                                                           */
    /* Additional variables can be stored in the "mc_t_distrib_out" data set     */
    /* (see the macro parameter descriptions for details).                       */
    /*                                                                           */
    /* If an "optional_input_mcmc_u_out_data" SAS data set is not specified, the */
    /* number of observations in the output data set is the product of the       */
    /* "number of individuals" and the specified value for the                   */
    /* "set_number_monte_carlo_rand_obs" macro parameter (see the macro          */
    /* parameter descriptions for details).                                      */
    /*                                                                           */
    /*****************************************************************************/
    /*****************************************************************************/
    /*****************************************************************************/

    ****************************************************;
    *** Call the MULTIVAR_DISTRIB macro.             ***;
    *** Note that the multivariate measurement error ***;
    *** model was fit using the MULTIVAR_MCMC macro. ***;
    ****************************************************;

    %do subset_monte_carlo_distrib = 1 %to &num_subsets_monte_carlo_distrib;

      %let loopseed_distrib = %eval(89009890 + (&replicnum * 10000) + (&subset_monte_carlo_distrib * 111));

      data outlib.sub_multivar_post_mcmc_u_r&replicnum;
        set outlib.multivar_post_mcmc_u_out_r&replicnum;
        if (&subset_number_post_mcmc_u_out * (&subset_monte_carlo_distrib - 1)) < iteration <= (&subset_number_post_mcmc_u_out * &subset_monte_carlo_distrib);
        iteration = iteration - (&subset_number_post_mcmc_u_out * (&subset_monte_carlo_distrib - 1));
      run;

      %multivar_distrib(multivar_mcmc_out_lib           = outlib,
                        multivar_mcmc_out_store_label   = epi_mcmc_nc_main24hr&replicnum,
                        t_weightavg_covariates_list1    = constant1 ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g std_boxcox_ffq_redmeat_g
                                                          std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg std_boxcox_ffq_potassium_mg
                                                          std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal std_entry_age std_log_bmi,
                        t_weightavg_covariates_list2    = ,
                        set_value_for_weight_cov_list1  = ,
                        set_value_for_weight_cov_list2  = ,
                        nev_consumers_epis1             = y,
                        covars_prob_consumer_epis1      = constant1 ffq_eqzero_alcohol_g,
                        optional_input_data             = subj1rec,
                        optional_input_data_var_list    = replicaterowid constant1 ffq_eqzero_alcohol_g std_boxcox_ffq_0mean_alcohol_g
                                                          std_boxcox_ffq_redmeat_g std_boxcox_ffq_whlgrain_ozeq std_boxcox_ffq_sodium_mg
                                                          std_boxcox_ffq_potassium_mg std_boxcox_ffq_protein_g std_boxcox_ffq_energy_kcal
                                                          std_entry_age std_log_bmi,
                        optional_input_mcmc_u_out_data  = sub_multivar_post_mcmc_u_r&replicnum,
                        additional_output_var_list      = ,
                        additional_output_subject_var   = replicaterowid,
                        output_mcmc_weight_var          = ,
                        set_seed_distrib                = &loopseed_distrib,
                        set_number_monte_carlo_rand_obs = ,
                        print                           = y
                        );

      data allsubsets_mc_t_distrib_out;
        %if &subset_monte_carlo_distrib = 1 %then %do;
          set mc_t_distrib_out;
        %end;
        %else %do;
          set allsubsets_mc_t_distrib_out mc_t_distrib_out;
        %end;
      run;

    %end;


    *********************************;
    *** Delete unneeded data sets ***;
    *********************************;

    proc datasets lib=outlib nolist;
      delete sub_multivar_post_mcmc_u_r&replicnum multivar_post_mcmc_u_out_r&replicnum;
    quit;
    run;


    ***************************************************;
    *** Rename the MC distribution output variables ***;
    ***************************************************;

    data usintake_mc_t(drop = mc_t1 - mc_t7);
      set allsubsets_mc_t_distrib_out;

      mc_alcohol   = mc_t1;
      mc_redmeat   = mc_t2;
      mc_whlgrain  = mc_t3;
      mc_sodium    = mc_t4;
      mc_potassium = mc_t5;
      mc_protein   = mc_t6;
      mc_energy    = mc_t7;
    run;



    ******************************************************************;
    *** Calculate mean of predicted intake values for each subject ***;
    ******************************************************************;

    proc sort data=usintake_mc_t;
      by replicaterowid;
    run;

    proc means data=usintake_mc_t noprint;
      by replicaterowid;
      var mc_alcohol mc_redmeat mc_whlgrain mc_sodium mc_potassium mc_protein mc_energy;
      output out=predint_out mean=predint_alcohol predint_redmeat predint_whlgrain predint_sodium predint_potassium
                                  predint_protein predint_energy;
    run;



    ******************************************************************************;
    *** Merge mean of predicted intake values back onto the replicate data set ***;
    ******************************************************************************;

    data subj1recres;
      merge subj1rec predint_out(keep = replicaterowid predint_alcohol predint_redmeat predint_whlgrain predint_sodium
                                        predint_potassium predint_protein predint_energy);
      by replicaterowid;
    run;



    *****************************************************************************************;
    *** Standardize the covariates so that one unit equals the distance between t0 and t1 ***;
    *****************************************************************************************;

    data subj1recres;
      set subj1recres;

      unit_alcohol = &t1_alcohol - &t0_alcohol;
      predint_alcohol_unit_iqr = predint_alcohol / unit_alcohol;

      unit_redmeat = &t1_redmeat - &t0_redmeat;
      predint_redmeat_unit_iqr = predint_redmeat / unit_redmeat;

      unit_whlgrain = &t1_whlgrain - &t0_whlgrain;
      predint_whlgrain_unit_iqr = predint_whlgrain / unit_whlgrain;

      unit_sodium = &t1_sodium - &t0_sodium;
      predint_sodium_unit_iqr = predint_sodium / unit_sodium;

      unit_potassium = &t1_potassium - &t0_potassium;
      predint_potassium_unit_iqr = predint_potassium / unit_potassium;

      unit_protein = &t1_protein - &t0_protein;
      predint_protein_unit_iqr = predint_protein / unit_protein;

      unit_energy = &t1_energy - &t0_energy;
      predint_energy_unit_iqr = predint_energy / unit_energy;

      unit_std_entry_age = &t1_std_entry_age - &t0_std_entry_age;
      std_entry_age_unit_iqr = std_entry_age / unit_std_entry_age;

      unit_std_log_bmi = &t1_std_log_bmi - &t0_std_log_bmi;
      std_log_bmi_unit_iqr = std_log_bmi / unit_std_log_bmi;
    run;



    ********************************************;
    *** Fit a logistic regression risk model ***;
    ********************************************;

    proc logistic data=subj1recres outest=regout;
      model ybin_event(event='1') = predint_alcohol_unit_iqr predint_redmeat_unit_iqr predint_whlgrain_unit_iqr
                                    predint_sodium_unit_iqr predint_potassium_unit_iqr predint_protein_unit_iqr
                                    predint_energy_unit_iqr std_entry_age_unit_iqr std_log_bmi_unit_iqr;
      title5 "Logistic Regression Risk Model: Adjusted for Measurement Error";
    run;
    title5;



    *************************************************************;
    *** Save the regression results for the current replicate ***;
    *************************************************************;

    data outlib.diethealthout&replicnum(keep = replicate intercept rc_beta_alcohol_unit_iqr rc_or_alcohol_unit_iqr
                                               rc_beta_redmeat_unit_iqr rc_or_redmeat_unit_iqr rc_beta_whlgrain_unit_iqr
                                               rc_or_whlgrain_unit_iqr rc_beta_sodium_unit_iqr rc_or_sodium_unit_iqr
                                               rc_beta_potassium_unit_iqr rc_or_potassium_unit_iqr rc_beta_protein_unit_iqr
                                               rc_or_protein_unit_iqr rc_beta_energy_unit_iqr rc_or_energy_unit_iqr
                                               rc_beta_std_entry_age_unit_iqr rc_or_std_entry_age_unit_iqr rc_beta_std_log_bmi_unit_iqr
                                               rc_or_std_log_bmi_unit_iqr t0_alcohol t1_alcohol t0_redmeat t1_redmeat t0_whlgrain
                                               t1_whlgrain t0_sodium t1_sodium t0_potassium t1_potassium t0_protein t1_protein t0_energy
                                               t1_energy t0_std_entry_age t1_std_entry_age t0_std_log_bmi t1_std_log_bmi);

      set regout(rename = (predint_alcohol_unit_iqr   = rc_beta_alcohol_unit_iqr
                           predint_redmeat_unit_iqr   = rc_beta_redmeat_unit_iqr
                           predint_whlgrain_unit_iqr  = rc_beta_whlgrain_unit_iqr
                           predint_sodium_unit_iqr    = rc_beta_sodium_unit_iqr
                           predint_potassium_unit_iqr = rc_beta_potassium_unit_iqr
                           predint_protein_unit_iqr   = rc_beta_protein_unit_iqr
                           predint_energy_unit_iqr    = rc_beta_energy_unit_iqr
                           std_entry_age_unit_iqr     = rc_beta_std_entry_age_unit_iqr
                           std_log_bmi_unit_iqr       = rc_beta_std_log_bmi_unit_iqr
                           ));

      replicate = &replicnum;



      *********************************;
      *** Calculate the odds ratios ***;
      *********************************;

      rc_or_alcohol_unit_iqr       = exp(rc_beta_alcohol_unit_iqr);
      rc_or_redmeat_unit_iqr       = exp(rc_beta_redmeat_unit_iqr);
      rc_or_whlgrain_unit_iqr      = exp(rc_beta_whlgrain_unit_iqr);
      rc_or_sodium_unit_iqr        = exp(rc_beta_sodium_unit_iqr);
      rc_or_potassium_unit_iqr     = exp(rc_beta_potassium_unit_iqr);
      rc_or_protein_unit_iqr       = exp(rc_beta_protein_unit_iqr);
      rc_or_energy_unit_iqr        = exp(rc_beta_energy_unit_iqr);
      rc_or_std_entry_age_unit_iqr = exp(rc_beta_std_entry_age_unit_iqr);
      rc_or_std_log_bmi_unit_iqr   = exp(rc_beta_std_log_bmi_unit_iqr);



      ****************************************************************************;
      *** Save the comparison values, t0 and t1, for each risk model covariate ***;
      ****************************************************************************;

      t0_alcohol = &t0_alcohol;
      t1_alcohol = &t1_alcohol;

      t0_redmeat = &t0_redmeat;
      t1_redmeat = &t1_redmeat;

      t0_whlgrain = &t0_whlgrain;
      t1_whlgrain = &t1_whlgrain;

      t0_sodium = &t0_sodium;
      t1_sodium = &t1_sodium;

      t0_potassium = &t0_potassium;
      t1_potassium = &t1_potassium;

      t0_protein = &t0_protein;
      t1_protein = &t1_protein;

      t0_energy = &t0_energy;
      t1_energy = &t1_energy;

      t0_std_entry_age = &t0_std_entry_age;
      t1_std_entry_age = &t1_std_entry_age;

      t0_std_log_bmi = &t0_std_log_bmi;
      t1_std_log_bmi = &t1_std_log_bmi;

    run;

  %end;   *** End of the replicate loop ***;

%mend fit_models_replicate_loop;



************************************************;
*** Call the fit_models_replicate_loop macro ***;
************************************************;

%fit_models_replicate_loop;
