/*****************************************************************************/
/*** SAS version 9.2 or higher and the following files are required:       ***/
/***   boxcox_survey.macro.v1.2.sas                                        ***/
/***   std_cov_boxcox24hr_conday_minamt_macro_v1.0.sas                     ***/
/***   multivar_mcmc_macro_v2.1.sas                                        ***/
/***   multivar_distrib_macro_v2.1.sas                                     ***/
/***   percentiles_survey.macro.v1.1.sas                                   ***/
/***   t_0104nc_may7.sas7bdat                                              ***/
/***                                                                       ***/
/*** This program fits a multivariate measurement error model and          ***/
/*** calculates mean usual intake using NHANES 2003-2004 data.  It also    ***/
/*** calculates percentiles and cutpoint probablities for usual intake,    ***/
/*** conditional on quintiles of sodium intake.  For this analysis, the    ***/
/*** example data include males, ages 9+ from NHANES 2003-2004, 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, meat, whole grains, and sodium.               ***/
/***                                                                       ***/
/*** This analysis uses balanced repeated replication (BRR) variance       ***/
/*** estimation, so the multivariate measurement error model is fit using  ***/
/*** the original weights (i.e. weights for replicate 0) and using the 32  ***/
/*** BRR weights (i.e. weights for replicates 1, 2, ..., 32).  The mean    ***/
/*** usual intakes, usual intake conditional distributions, and cutpoint   ***/
/*** probabilities are also calculated for each replicate.  The following  ***/
/*** replicfirst and repliclast macro variables allow specification of a   ***/
/*** range of replicates.  This example program considers replicate 0.     ***/
/*****************************************************************************/

%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 "Calculate Usual Intake Conditional Distributions and Cutpoint Probabilities";



***********************************************************************;
*** 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_v1.0.sas";
%include "&home/include.files.macros/multivar_mcmc_macro_v2.1.sas";
%include "&home/include.files.macros/multivar_distrib_macro_v2.1.sas";
%include "&home/include.files.macros/percentiles_survey.macro.v1.1.sas";

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

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



************************************************************************************;
*** A macro is used to loop through the original weight variable and the 32 BRR  ***;
*** weights.  To accomplish this seamlessly the name of the original weight      ***;
*** variable and the BRR weight variables must have the same root and must end   ***;
*** in 0-32, e.g. rndw0-rndw32, where rndw0 is the original weight. (In this     ***;
*** example the original weight is rndw0.) The following macro variables are     ***;
*** used to ensure proper naming.                                                ***;
************************************************************************************;

%let weight_var_orig = rndw;
%let weight_var      = rndw;



******************************;
*** Read in the input data ***;
******************************;

data nhanes;
  set inlib.t_0104nc_may7;
run;



***************************************************************************************************;
*** Create separate data sets for the BOXCOX_SURVEY and STD_COV_BOXCOX24HR_CONDAY_MINAMT macros ***;
***************************************************************************************************;

data nhanes(keep = stratum id rndw0 rndw1-rndw32 repeat repeat2 weekend agegrp age_1to3 age_4to8
                   age_9to13 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus a_bev m_meat
                   g_whl drxtsodi)

     nhanes1rec(keep = stratum id rndw0 rndw1-rndw32 repeat repeat2 weekend agegrp age_1to3 age_4to8
                   age_9to13 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus a_bev m_meat
                   g_whl drxtsodi);

  set nhanes;


  *************************;
  *** Rename covariates ***;
  *************************;

  stratum   = stra;
  id        = seqn;
  repeat    = drddaycd;
  weekend   = wkend;


  ***********************************************************************;
  *** The survey weight variable name does not end with zero (0),     ***;
  *** so it is renamed.  The new name has the same root name of the   ***;
  *** BRR weight variables and has 0 appended.                        ***;
  *** In this example the original weight variable name is rndw.      ***;
  *** The BRR weights are named rndw1-rndw32.                         ***;
  ***********************************************************************;

  rename &weight_var_orig = &weight_var.0;


  ****************************************;
  *** Subset to the strata of interest ***;
  ****************************************;

  if stratum = 2;  ***Males, 9+;


  ***************************************;
  *** Create covariates for the model ***;
  ***************************************;

  if repeat = 2 then repeat2 = 1;
  else repeat2 = 0;

  array a (*) age_1to3 age_4to8 age_9to13 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus;
  do i = 1 to dim(a);
    a(i) = 0;
  end;
  a(agegrp) = 1;


  **********************;
  *** Output records ***;
  **********************;

  if repeat = 1 then output nhanes1rec;
  output nhanes;

run;



***************************************************************************;
*** The macro fit_models_replicate_loop is used to call the MCMC macros ***;
*** using the original weight and each of the 32 BRR weights            ***;
***************************************************************************;

%macro fit_models_replicate_loop;


  %do replicnum = &replicfirst %to &repliclast;

    **********************************************;
    *** Use macro variable to represent the    ***;
    *** replicate variables in subsequent code ***;
    **********************************************;

    %let replicvar = &weight_var&replicnum;


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

    data nhanes1rec_pos_alcohol nhanes1rec_pos_meat nhanes1rec_pos_whlgrain;
      set nhanes1rec;

      if a_bev  > 0 then output nhanes1rec_pos_alcohol;
      if m_meat > 0 then output nhanes1rec_pos_meat;
      if g_whl  > 0 then output nhanes1rec_pos_whlgrain;
    run;



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

    %boxcox_survey(data    = nhanes1rec_pos_alcohol,
                   subject = id,
                   var     = a_bev,
                   covars  = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   weight  = &replicvar,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = nhanes1rec_pos_meat,
                   subject = id,
                   var     = m_meat,
                   covars  = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   weight  = &replicvar,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = nhanes1rec_pos_whlgrain,
                   subject = id,
                   var     = g_whl,
                   covars  = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   weight  = &replicvar,
                   print   = N,
                   plot    = N,
                   ntitle  = 4
                   );

    %boxcox_survey(data    = nhanes1rec,
                   subject = id,
                   var     = drxtsodi,
                   covars  = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   weight  = &replicvar,
                   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 nhanes1rec_pos_alcohol(obs=1 in=alcohol)
          nhanes1rec_pos_meat(obs=1 in=meat)
          nhanes1rec_pos_whlgrain(obs=1 in=whlgrain)
          nhanes1rec(obs=1 in=nhanes1rec);

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

      else if meat then do;
        tran_lambda = lambda_m_meat;
        tran_paramindex = 2;
        output;
      end;

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

      else if nhanes1rec then do;
        tran_lambda = lambda_drxtsodi;
        tran_paramindex = 4;
        output;
      end;

    run;



    /*****************************************************************************/
    /*** Use the STD_COV_BOXCOX24HR_CONDAY_MINAMT macro to transform and       ***/
    /*** standardize reported nonzero amounts and create variables needed for  ***/
    /*** the MULTIVAR_MCMC macro.                                              ***/
    /*****************************************************************************/

    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* 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                       = nhanes,
                                      prestand_continuous_covars = ,
                                      rec24hr_epis_vars          = a_bev m_meat g_whl,
                                      rec24hr_daily_vars         = drxtsodi,
                                      boxcox_tran_lambda_data    = boxcox_lambda_data,
                                      minamount_data             = ,
                                      print                      = y,
                                      titles                     = 2
                                      );



    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* 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.                                 */
    /*                                                                           */
    /*****************************************************************************/
    /*****************************************************************************/
    /*****************************************************************************/


    %let graph_pdf_name = surveillance_mcmc_nc_main24hr_traceplot.rep&replicnum..pdf;
    %let loopseed_mcmc = %eval(42941 + &replicnum * 10);


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

    %multivar_mcmc(data                        = stdcov_stdbc24hr_conday_out,
                   subject                     = id,
                   weight_var                  = &replicvar,
                   repeat                      = repeat,
                   conday_epis_vars            = conday_a_bev conday_m_meat conday_g_whl,
                   gst_rec24hr_epis_vars       = stdbc_a_bev stdbc_m_meat stdbc_g_whl,
                   gst_rec24hr_daily_vars      = stdbc_drxtsodi,
                   covars_epis_prob            = constant1 repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   covars_epis_amt             = constant1 repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   covars_daily_amt            = constant1 repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                   nev_consumers_epis1         = y,
                   covars_prob_consumer_epis1  = constant1,
                   set_seed_mcmc               = &loopseed_mcmc,
                   set_number_mcmc_iterations  = ,
                   set_number_burn_iterations  = ,
                   set_thin                    = ,
                   prior_sigmau_mean_data      = ,
                   sigmau_constant             = ,
                   gen_inverse                 = ,
                   print                       = y,
                   titles                      = 2,
                   std_print_store             = ,
                   notes_print                 = y,
                   out_lib                     = outlib,
                   out_store_label             = surv_mcmc_nc_main24hr&replicnum,
                   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 = &graph_pdf_name,
                   optional_iml_store_data     = backtran_out,
                   optional_iml_store_names    = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus 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. ***;
    ****************************************************;

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

    %multivar_distrib(multivar_mcmc_out_lib           = outlib,
                      multivar_mcmc_out_store_label   = surv_mcmc_nc_main24hr&replicnum,
                      t_weightavg_covariates_list1    = constant1 constant0 constant0 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                      t_weightavg_covariates_list2    = constant1 constant0 constant1 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                      set_value_for_weight_cov_list1  = 4,
                      set_value_for_weight_cov_list2  = 3,
                      nev_consumers_epis1             = y,
                      covars_prob_consumer_epis1      = constant1,
                      optional_input_data             = ,
                      optional_input_data_var_list    = ,
                      optional_input_mcmc_u_out_data  = ,
                      additional_output_var_list      = ,
                      additional_output_subject_var   = id,
                      output_mcmc_weight_var          = ,
                      set_seed_distrib                = &loopseed_distrib,
                      set_number_monte_carlo_rand_obs = 500,
                      print                           = y
                      );


    data usintake_mc_t;
      set mc_t_distrib_out;

      ***********************************************************************;
      *** Divide weights by set_number_monte_carlo_rand_obs value used in ***;
      *** the MULTIVAR_DISTRIB macro.                                     ***;
      ***********************************************************************;

      weight_nw_sumw_div500 = weight_nw_sumw / 500;


      ********************************************************;
      *** Assign variable names for usual intakes from the ***;
      *** MULTIVAR_DISTRIB macro.                          ***;
      ********************************************************;

      usintake_alcohol  = mc_t1;
      usintake_meat     = mc_t2;
      usintake_whlgrain = mc_t3;
      usintake_sodium   = mc_t4;

      drop mc_t1-mc_t4;
    run;



    ***************************************************;
    *** Group records by quintiles of sodium intake ***;
    ***************************************************;

    proc univariate data=usintake_mc_t noprint;
      var usintake_sodium;
      output pctlpre=P pctlpts = 20,40,60,80 out=pctls;
    run;

    data usintake_mc_t;
      if _n_ = 1 then set pctls;
      set usintake_mc_t;

      if      usintake_sodium <= P20 then sodium_quintile = 1;
      else if usintake_sodium <= P40 then sodium_quintile = 2;
      else if usintake_sodium <= P60 then sodium_quintile = 3;
      else if usintake_sodium <= P80 then sodium_quintile = 4;
      else                                sodium_quintile = 5;
    run;

    proc sort data=usintake_mc_t;
      by sodium_quintile;
    run;



    *********************************************************************************;
    *** Create a list of the dietary components for the percentiles and cutpoints ***;
    *********************************************************************************;

    %let dc_list = usintake_alcohol usintake_meat usintake_whlgrain;
    %let cutpoint_list = %str(0.1 0.5 1.0,
                              1.0 1.5 2.0,
                              0.2 0.4 0.8);


    *********************************************************;
    *** Call the PERCENTILES_SURVEY macro for each of the ***;
    *** dietary components, by quintiles of sodium intake ***;
    *********************************************************;

    %do scani = 1 %to 3;

      %let dc = %scan(&dc_list, &scani, %str( ));
      %let cutpoints = %scan(&cutpoint_list, &scani, %str(,));

      %percentiles_survey(data      = usintake_mc_t,
                          byvar     = sodium_quintile,
                          var       = &dc,
                          weight    = weight_nw_sumw_div500,
                          cutpoints = &cutpoints,
                          print     = N,
                          ntitle    = 3
                          );


      ********************************************************************************;
      ** Concatenate output for every dietary component on each replicate data set ***;
      ********************************************************************************;

      data _percentiles;
        length dc $20;
        set _percentiles;
        dc = "&dc";

        %do cpi = 1 %to 3;
          %let cutpoint = %scan(&cutpoints, &cpi, %str( ));
          Cutpoint&cpi = &cutpoint;
        %end;
      run;

      %if &scani = 1 %then %do;
        data _percentiles_all_dc&replicnum;
          set _percentiles;
        run;
      %end;
      %else %do;
        proc append base=_percentiles_all_dc&replicnum data=_percentiles force;
        run;
      %end;

      proc datasets nolist lib=work;
        delete _percentiles;
      quit;

    %end;



    **************************************************;
    *** Save the final data set for each replicate ***;
    **************************************************;

    data outlib._percentiles_all_dc&replicnum;
      set _percentiles_all_dc&replicnum;
      replicate=&replicnum;
    run;

    proc datasets nolist lib=work;
      delete _percentiles
             _percentiles_all_dc&replicnum;
    quit;

  %end; *** end loop for replicates ***;

%mend;



%fit_models_replicate_loop;
