/*****************************************************************************/
/*** SAS version 9.2 or higher and the following files are required for    ***/
/*** this example:                                                         ***/
/***   std_cov_boxcox24hr_conday_minamt_macro_v1.0.sas                     ***/
/***   multivar_mcmc_macro_v2.1.sas                                        ***/
/***   data.hei2010.nonsmok.smok.ag20p.xpt                                 ***/
/***   bclam.hei2010.nonsmok.smok.ag20p.xpt.                               ***/
/***                                                                       ***/
/*** This example includes 5 programs.  Programs 1a, 1b, 1c can be used to ***/
/*** fit a multivariate measurement error model, calculate mean usual      ***/
/*** Healthy Eating Index–2010 (HEI–2010) scores, and perform t-tests      ***/
/*** comparing mean usual HEI–2010 scores using NHANES 2003-2004 data.     ***/
/*** Programs 1d and 1e can be used to extend the analysis to also         ***/
/*** calculate percentiles and cutpoint probablities for both the HEI-2010 ***/
/*** component scores and the HEI-2010 densities, conditional on           ***/
/*** categories of HEI-2010 total scores.                                  ***/
/***                                                                       ***/
/*** This SAS program illustrates the use of the                           ***/
/*** STD_COV_BOXCOX24HR_CONDAY_MINAMT and MULTIVAR_MCMC macros to fit a    ***/
/*** multivariate measurement error model for the Healthy Eating HEI–2010  ***/
/*** using NHANES 2003-2004 data for smokers, age 20 and older.  This      ***/
/*** program can also be modified to fit a multivariate measurement error  ***/
/*** model for nonsmokers, age 20 and older.                               ***/
/***                                                                       ***/
/*** The HEI-2010 nonsmokers and smokers 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 16 BRR weights (i.e. weights ***/
/*** for replicates 1, 2, ..., 16).  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;


/***************************************************************/
/*** Specify a title and label for this example and stratum. ***/
/***************************************************************/

title1 "Fit Multivariate Measurement Error Model Using MCMC with 24-Hour Recall as Main Instrument";
title2 "Calculate Mean Usual HEI–2010 Scores and Perform T-Tests Comparing Mean Usual HEI–2010 Scores";
title3 "Calculate Conditional Distributions and Cutpoint Probabilities for HEI-2010 Densities and Scores";
%let example_stratum_label = ex1smok;


/*****************************************************************************/
/*** 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.  Note that this library is specified    ***/
/*** using the label provided in the above definition of the               ***/
/*** example_stratum_label macro variable.                                 ***/
/*****************************************************************************/

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

*** Include the required macros ***;
%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";

*** Input data ***;
filename trandat "&home/multivar_surveillance_example1_mcmc_nnc_main24hr/data/data.hei2010.nonsmok.smok.ag20p.xpt";
filename tranlam "&home/multivar_surveillance_example1_mcmc_nnc_main24hr/data/bclam.hei2010.nonsmok.smok.ag20p.xpt";

*** Output data library ***;
libname outlib "&home/multivar_surveillance_example1_mcmc_nnc_main24hr/outlib/&example_stratum_label";



/*****************************************************************************/
/*** Import the data file and file of transformation parameter values for  ***/
/*** this example and create separate data sets for nonsmokers and smokers ***/
/*****************************************************************************/

proc cimport data=data_hei2010_nonsmok_smok_ag20p infile=trandat;
run;

proc cimport data=bclam_hei2010_nonsmok_smok_ag20p infile=tranlam;
run;

data data_hei2010_nonsmok_ag20p;
  set data_hei2010_nonsmok_smok_ag20p(where = (smoker=0));
run;

data data_hei2010_smok_ag20p;
  set data_hei2010_nonsmok_smok_ag20p(where = (smoker=1));
run;

data bclam_hei2010_nonsmok_ag20p(keep = tran_paramindex tran_lambda);
  set bclam_hei2010_nonsmok_smok_ag20p(where = (smoker=0));
run;

data bclam_hei2010_smok_ag20p(keep = tran_paramindex tran_lambda);
  set bclam_hei2010_nonsmok_smok_ag20p(where = (smoker=1));
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                       = data_hei2010_smok_ag20p,
                                  prestand_continuous_covars = ,
                                  rec24hr_epis_vars          = frtjuice wholefrt obeangrn olegumes g_whl oseaplant,
                                  rec24hr_daily_vars         = ononbeangrn g_nwhl d_total ononseaplant tsodi tmfat
                                                               tpfat tsfat emptycal10 nonemptycal10,
                                  boxcox_tran_lambda_data    = bclam_hei2010_smok_ag20p,
                                  minamount_data             = ,
                                  print                      = y,
                                  titles                     = 3
                                  );



/*****************************************************************************/
/*** Use the FIT_MODELS_REPLICATE_LOOP macro to call the MULTIVAR_MCMC     ***/
/*** macro for each of the replicates specified in the above definition of ***/
/*** the replicfirst and repliclast macro variables.                       ***/
/*****************************************************************************/

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



%macro fit_models_replicate_loop;

  %do replicnum = &replicfirst %to &repliclast;

    %let loopseed = %eval(42941 + &replicnum * 1000);
    %let graph_pdf_name = &example_stratum_label..trace.replic&replicnum..pdf;


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

    %multivar_mcmc(data                        = stdcov_stdbc24hr_conday_out,
                   subject                     = seqn,
                   weight_var                  = w0304_&replicnum,
                   repeat                      = recseq,
                   conday_epis_vars            = conday_frtjuice conday_wholefrt conday_obeangrn conday_olegumes
                                                 conday_g_whl conday_oseaplant,
                   gst_rec24hr_epis_vars       = stdbc_frtjuice stdbc_wholefrt stdbc_obeangrn stdbc_olegumes
                                                 stdbc_g_whl stdbc_oseaplant,
                   gst_rec24hr_daily_vars      = stdbc_ononbeangrn stdbc_g_nwhl stdbc_d_total stdbc_ononseaplant
                                                 stdbc_tsodi stdbc_tmfat stdbc_tpfat stdbc_tsfat stdbc_emptycal10
                                                 stdbc_nonemptycal10,
                   covars_epis_prob            = constant1 recallday2 weekend age20_30 age31_50 age51_70
                                                 raceeth_mexam raceeth_nhsblk male male_age20_30 male_age31_50
                                                 male_age51_70 male_raceeth_mexam male_raceeth_nhsblk,
                   covars_epis_amt             = constant1 recallday2 weekend age20_30 age31_50 age51_70
                                                 raceeth_mexam raceeth_nhsblk male male_age20_30 male_age31_50
                                                 male_age51_70 male_raceeth_mexam male_raceeth_nhsblk,
                   covars_daily_amt            = constant1 recallday2 weekend age20_30 age31_50 age51_70
                                                 raceeth_mexam raceeth_nhsblk male male_age20_30 male_age31_50
                                                 male_age51_70 male_raceeth_mexam male_raceeth_nhsblk,
                   set_seed_mcmc               = &loopseed,
                   set_number_mcmc_iterations  = ,
                   set_number_burn_iterations  = ,
                   set_thin                    = ,
                   prior_sigmau_mean_data      = ,
                   sigmau_constant             = ,
                   gen_inverse                 = y,
                   print                       = y,
                   titles                      = 3,
                   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    = constant1 recallday2 weekend age20_30 age31_50 age51_70
                                                 raceeth_mexam raceeth_nhsblk male male_age20_30 male_age31_50
                                                 male_age51_70 male_raceeth_mexam male_raceeth_nhsblk
                                                 tran_paramindex tran_lambda tran_center tran_scale minamount
                   );


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


%mend;



%fit_models_replicate_loop;
