options mergenoby=nowarn;

/***************************************************************************/
/* The following files are required to run this example:                   */
/*   nlmixed_bivariate_macro_v1.2.sas                                      */
/*   t_0104nc_may7.sas7bdat                                                */
/***************************************************************************/
/*                                                                         */
/***************************************************************************/
/* This example analysis fits univariate and bivariate measurement error   */
/* models for two dietary components that are consumed almost daily.  It   */
/* also estimates the distribution of the ratio of usual intakes for two   */
/* dietary components.  For this analysis, the example data is from the    */
/* 2001-2004 NHANES public use data and includes children ages 1-8, males  */
/* aged 9+ years, and females aged 9+ years.  The 24-hour dietary recall   */
/* is the main instrument and the two dietary components are saturated fat */
/* and energy.  The example data include repeated intake measurements from */
/* 24-hour dietary recalls.                                                */
/*                                                                         */
/* This analysis uses balanced repeated replication (BRR) variance         */
/* estimation, so the bivariate measurement error model and distribution   */
/* of usual intake use the original weight and 32 BRR weights.  The        */
/* following replicfirst and repliclast macro variables allow              */
/* specification of a range of replicates.  This program considers the     */
/* original data set (i.e. replicate 0).                                   */
/***************************************************************************/

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



title1 "Fit Bivariate Measurement Error Model Using MLE with 24-Hour Recall as Main Instrument";
title2 "Estimate Distribution of the Ratio of Usual Intakes for Two Dietary Components";



***********************************************************************;
*** 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/nlmixed_bivariate_macro_v1.2.sas";

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



***************;
*** Formats ***;
***************;

proc format;
  value sexfmt
    1 = "Male"
    2 = "Female";
  value stratfmt
    1 = "Children 1-8"
    2 = "Males 9+"
    3 = "Females 9+";
run;



************************************************************************************;
*** A macro is used to loop the original weight variable and the 32 BRR weights  ***;
*** through calls to the NLMIXED_UNIVARIATE macro.  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;



data nhanes (keep = stratum id sex female age race repeat repeat2 dayofweek weekend agegrp rndw0 rndw1-rndw32 recall_dc1
                    recall_dc2 age_1to3 age_4to8 age_9to13 age_14to18 age_19to30 age_31to50 age_51to70 age_71plus);

  set outlib.t_0104nc_may7;


  ***********************************************************************;
  *** 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;


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

  stratum   = stra;
  id        = seqn;
  sex       = riagendr;
  age       = ridageyr;
  race      = ridreth1;
  repeat    = drddaycd;
  dayofweek = drdday;
  weekend   = wkend;


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

  if stratum in (1:3);


  ****************************************;
  *** Specify recall dietary component ***;
  ****************************************;

  recall_dc1 = DRXTSFAT;
  recall_dc2 = DRXTKCAL;


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

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

  if sex = 2 then female = 1;
  else if sex = 1 then female = 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;


  format sex sexfmt. stratum stratfmt.;

run;

proc sort data=nhanes;
  by stratum id repeat;
run;



**********************************************************************************;
*** The macro fit_models_replicate_loop is used to call the NLMIXED_BIVARIATE  ***;
*** macro and fit a bivariate measurement error model using the original study ***;
*** data (i.e. replicate data set 0) and each of the replicate datasets (i.e.  ***;
*** replicate datasets 1, 2, ..., 32) for balanced repeated replication (BRR)  ***;
*** variance estimation                                                        ***;
**********************************************************************************;

%macro fit_models_replicate_loop;

  %do replicnum = &replicfirst %to &repliclast;

    title4 "Replicate &replicnum";

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

    %let replicvar=&weight_var&replicnum ;



    ***********************************************************************************************;
    *** Get initial parameter estimates from univariate models for the first dietary component. ***;
    *** Rename the univariate parameters for bivariate model.                                   ***;
    ***********************************************************************************************;

    data init_parms_dc1;
      set outlib.parms_u_sfat&replicnum;

      rename a_intercept        = A1_intercept
             a_repeat2          = A1_repeat2
             a_weekend          = A1_weekend
             a_female           = A1_female
             a_age_4to8         = A1_age_4to8
             a_age_14to18       = A1_age_14to18
             a_age_19to30       = A1_age_19to30
             a_age_31to50       = A1_age_31to50
             a_age_51to70       = A1_age_51to70
             a_age_71plus       = A1_age_71plus
             a_logsde           = A1_LogSDe
             a_lambda           = A1_Lambda
             logsdu2            = LogSDu2
             min_a              = min_a1;
    run;



    ************************************************************************************************;
    *** Get initial parameter estimates from univariate models for the second dietary component. ***;
    *** Rename the univariate parameters for bivariate model.                                    ***;
    ************************************************************************************************;

    data init_parms_dc2;
      set outlib.parms_u_energy&replicnum;

      rename a_intercept        = A2_intercept
             a_repeat2          = A2_repeat2
             a_weekend          = A2_weekend
             a_female           = A2_female
             a_age_4to8         = A2_age_4to8
             a_age_14to18       = A2_age_14to18
             a_age_19to30       = A2_age_19to30
             a_age_31to50       = A2_age_31to50
             a_age_51to70       = A2_age_51to70
             a_age_71plus       = A2_age_71plus
             a_logsde           = A2_LogSDe
             a_lambda           = A2_Lambda
             logsdu2            = LogSDu3
             min_a              = min_a2;
      run;



    ******************************************************************************;
    *** Combine the initial parameter estimates for the two dietary components ***;
    ******************************************************************************;

    data init_parms;
      merge init_parms_dc1 init_parms_dc2;
        by stratum;

      keep stratum
           A1_intercept A1_repeat2 A1_weekend A1_female A1_age_4to8 A1_age_14to18 A1_age_19to30 A1_age_31to50 A1_age_51to70 A1_age_71plus
           A2_intercept A2_repeat2 A2_weekend A2_female A2_age_4to8 A2_age_14to18 A2_age_19to30 A2_age_31to50 A2_age_51to70 A2_age_71plus
           A1_LogSDe A2_LogSDe LogSDu2 LogSDu3;

      format stratum stratfmt.;
    run;



    ***************************************************************************************************;
    *** Add minimum amount on consumption day to data set.                                          ***;
    *** For one-part models and for energy, set 24hr zero values to half the minimum nonzero value. ***;
    ***************************************************************************************************;

    data min_a;
      merge init_parms_dc1 init_parms_dc2;
        by stratum;

      keep stratum min_a1 min_a2;
    run;

    data nhanes_min_a;
      merge nhanes min_a;
        by stratum;

      modeltype = "ONEPART";
      if modeltype = "ONEPART" & recall_dc1 = 0 then recall_dc1 = min_a1 / 2;
      if recall_dc2 = 0 then recall_dc2 = min_a2 / 2;
    run;



    *********************************************************************************************;
    *** Create macro variables for the Box-Cox transformations for the two dietary components ***;
    *********************************************************************************************;

    data lambdas;
      merge init_parms_dc1 init_parms_dc2;
        by stratum;

      if stratum = 1 then call symput("lambda11",strip(put(a1_lambda,4.2)));
      if stratum = 2 then call symput("lambda21",strip(put(a1_lambda,4.2)));
      if stratum = 3 then call symput("lambda31",strip(put(a1_lambda,4.2)));

      if stratum = 1 then call symput("lambda12",strip(put(a2_lambda,4.2)));
      if stratum = 2 then call symput("lambda22",strip(put(a2_lambda,4.2)));
      if stratum = 3 then call symput("lambda32",strip(put(a2_lambda,4.2)));

      keep stratum a1_lambda a2_lambda;
    run;



    ***********************************************************************************;
    *** Create separate data sets each stratum (children 1-8, males 9+, females 9+) ***;
    ***********************************************************************************;

    data nhanes1 nhanes2 nhanes3;
      set nhanes_min_a;
      if stratum = 1 then output nhanes1;
      if stratum = 2 then output nhanes2;
      if stratum = 3 then output nhanes3;
    run;


    data init_parms1 (drop = stratum A1_age_14to18 A1_age_19to30 A1_age_31to50 A1_age_51to70 A1_age_71plus
                                     A2_age_14to18 A2_age_19to30 A2_age_31to50 A2_age_51to70 A2_age_71plus)
         init_parms2 (drop = stratum A1_female A1_age_4to8 A2_female A2_age_4to8)
         init_parms3 (drop = stratum A1_female A1_age_4to8 A2_female A2_age_4to8);

      set init_parms;
      if stratum = 1 then output init_parms1;
      if stratum = 2 then output init_parms2;
      if stratum = 3 then output init_parms3;
    run;



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

    proc datasets lib=work nolist;
      delete nhanes_min_a init_parms_dc1 init_parms_dc2 init_parms;
    run;



    /*******************************************************************************/
    /*                                                                             */
    /* SAS macro NLMixed_Bivariate fits a bivariate model for two foods/nutrients. */
    /* The first food/nutrient can be episodically consumed or consumed every      */
    /* day, while the second food/nutrient is assumed to be consumed every day.    */
    /*                                                                             */
    /* Model for episodically consumed foods/nutrients (two-part model):           */
    /* For episodically consumed foods/nutrients, the macro fits a two-part        */
    /* nonlinear mixed model, where the first part is the probability to           */
    /* consume and the second part is the amount consumed on a consumption day.    */
    /* The model allows for covariates in each part, includes a random effect      */
    /* for each part, and allows the random effects to be correlated.              */
    /*                                                                             */
    /* Model for foods/nutrients consumed every day (one-part model):              */
    /* For foods/nutrients consumed every day, the macro fits a one-part           */
    /* nonlinear mixed model of the amount consumed (the probability to consume    */
    /* is assumed to be 1). The model allows for covariates and includes a         */
    /* random effect.                                                              */
    /*                                                                             */
    /* For a food/nutrient that is consumed nearly every day by nearly everyone,   */
    /* so that the number of zero values is small, it may be preferable to use     */
    /* the one-part (consumed every day) model, since the two-part model may       */
    /* have trouble modeling the probability to consume in such a situation.       */
    /*                                                                             */
    /* Note, however, that the one-part model requires all responses to be         */
    /* greater than zero (zero values are treated as missing values).              */
    /* Before fitting the one-part model to a food/nutrient that has some zero     */
    /* values, replace the zero values with a small positive value, such as        */
    /* half the smallest observed nonzero value.                                   */
    /*                                                                             */
    /* Note: Initial parameter estimates must be supplied by the user.             */
    /* They can be estimated using SAS macro NLMixed_Univariate.                   */
    /*                                                                             */
    /* The macro calls the NLMixed procedure to fit the model.                     */
    /*                                                                             */
    /*******************************************************************************/
    /*                                                                             */
    /* Macro Parameters:                                                           */
    /*                                                                             */
    /*    Required Parameters:                                                     */
    /*       data          = name of SAS data set containing the data to be        */
    /*                       analyzed. The data set has multiple observations      */
    /*                       for each subject, one for each reptition of the       */
    /*                       24-hour recall (or other dietary instrument).         */
    /*       subject       = name of the variable that uniquely identifies each    */
    /*                       subject (i.e., ID variable).                          */
    /*       repeat        = name of the variable that indexes repeated            */
    /*                       observations for each subject.                        */
    /*       response1     = name of first food/nutrient variable to be modeled    */
    /*                       (24-hour recall variable for first food/nutrient).    */
    /*       response2     = name of second food/nutrient variable to be modeled   */
    /*                       (24-hour recall variable for second food/nutrient).   */
    /*       modeltype    = model for first food/nutrient:                         */
    /*                       to fit the two-part (epsisodic) model, specify        */
    /*                          modeltype = TWOPART                                */
    /*                       to fit the one-part (every day) model, specify        */
    /*                          modeltype = ONEPART                                */
    /*       init_parms    = name of SAS data set that contains initial            */
    /*                       parameter estimates. See the description of output    */
    /*                       data set parms_b (below) for further information.     */
    /*                                                                             */
    /*    Optional Parameters:                                                     */
    /*       covars_prob1  = list of variables that are covariates in the          */
    /*                       probability part of the two-part model for the        */
    /*                       first food/nutrient.                                  */
    /*                       if modeltype = ONEPART, then covars_prob is ignored.  */
    /*       covars_amt1   = list of variables that are covariates in the          */
    /*                       one-part model or the amount part of the              */
    /*                       two-part model for the first food/nutrient.           */
    /*       covars_amt2   = list of variables that are covariates in the          */
    /*                       one-part model for the second food/nutrient           */
    /*       link          = link function for the probability part of the two-    */
    /*                       part model for the first food/nutrient.               */
    /*                       to fit a logistic model, specify                      */
    /*                          link = logit                                       */
    /*                       to fit a probit model, specify                        */
    /*                          link = probit                                      */
    /*                       by default, link = probit.                            */
    /*                       if modeltype = ONEPART, then link is ignored.         */
    /*       lambda1       = Box-Cox transformation parameter for the first        */
    /*                       food/nutrient. If lambda1 is not specified, then      */
    /*                       it is estimated as part of the model.                 */
    /*       lambda2       = Box-Cox transformation parameter for the second       */
    /*                       food/nutrient. If lambda2 is not specified, then      */
    /*                       it is estimated as part of the model.                 */
    /*       var_u1        = variance of the random effect in the probability      */
    /*                       part of the model for the first food/nutrient.        */
    /*                       If var_u1 is not specified, then it is estimated      */
    /*                       as part of the model.                                 */
    /*                       if modeltype = ONEPART, then var_u1 is ignored.       */
    /*       var_u2        = variance of the random effect in the amount           */
    /*                       part of the model for the first food/nutrient.        */
    /*                       If var_u2 is not specified, then it is estimated      */
    /*                       as part of the model.                                 */
    /*       var_u3        = variance of the random effect in the amount           */
    /*                       part of the model for the second food/nutrient.       */
    /*                       If var_u3 is not specified, then it is estimated      */
    /*                       as part of the model.                                 */
    /*       indep_u1      = Y if random effect u1 is independent of u2 and u3.    */
    /*                     = N otherwise. by default, indep_u1 = N.                */
    /*                       if modeltype = ONEPART, then indep_u1 is ignored.     */
    /*       indep_u2      = Y if random effect u2 is independent of u1 and u3.    */
    /*                     = N otherwise. by default, indep_u2 = N.                */
    /*       indep_u3      = Y if random effect u3 is independent of u1 and u2.    */
    /*                     = N otherwise. by default, indep_u3 = N.                */
    /*       threshold     = Y to fit a latent variable threshold model.           */
    /*                     = N otherwise. by default, threshold = Y.               */
    /*                       if threshold = Y, then the probit model is fit.       */
    /*                       if modeltype = ONEPART, then threshold is ignored.    */
    /*       replicate_var = name of the sampling weight variable if the data      */
    /*                       is from a complex survey with weights.                */
    /*                       by default, the macro performs an unweighted          */
    /*                       analysis (assumes a simple random sample).            */
    /*       nloptions     = options for the NLMixed procedure that are            */
    /*                       appended to the PROC NLMIXED statement, e.g.,         */
    /*                        nloptions = technique=newrap maxiter=200,            */
    /*       print         = Y to print the output from the model.                 */
    /*                     = N to supress printing the output from the model.      */
    /*                     = V (verbose) to print extra output.                    */
    /*                       by default, print = Y.                                */
    /*       ntitle        = number of titles defined by the user.                 */
    /*                       by default, ntitle = 2.                               */
    /*                                                                             */
    /*******************************************************************************/
    /*                                                                             */
    /* Output Data Sets:                                                           */
    /*                                                                             */
    /*   parms_b = data set containing parameter estimates for the model.          */
    /*             parms_b contains the following variables:                       */
    /*                                                                             */
    /*                 A1_Intercept = intercept in the amount part of the model    */
    /*                                for the first food/nutrient.                 */
    /*                 A1_varname   = regression slope for covariate "varname"     */
    /*                                in the amount part of the model for the      */
    /*                                first food/nutrient.                         */
    /*                 A2_Intercept = intercept for the second food/nutrient.      */
    /*                 A2_varname   = regression slope for covariate "varname"     */
    /*                                for the second food/nutrient.                */
    /*                 A1_LogSDe    = Log(Sqrt(Var_e2))                            */
    /*                 A2_LogSDe    = Log(Sqrt(Var_e3))                            */
    /*                 z_e2e3       = Fisher transformation of Corr_e2e3:          */
    /*                                  z = ln[(1+corr)/(1-corr)] / 2              */
    /*                 Var_e2       = variance of within-person error e2 (amount   */
    /*                                part of model for first food/nutrient).      */
    /*                 Var_e3       = variance of within-person error e3 (second   */
    /*                                food/nutrient).                              */
    /*                 Var_u2       = variance of random effect u2 (amount part    */
    /*                                of model for first food/nutrient).           */
    /*                 Var_u3       = variance of random effect u3 (second         */
    /*                                food/nutrient).                              */
    /*                                of the model for the second food).           */
    /*                 Cov_e2e3     = covariance of random errors e2 and e3.       */
    /*                 Corr_e2e3    = correlation of random errors e2 and e3.      */
    /*                 Cov_u2u3     = covariance of random effects u2 and u3.      */
    /*                 Corr_u2u3    = correlation of random effects u2 and u3.     */
    /*                                                                             */
    /*             if fitting the two-part model for the first food/nutrient,      */
    /*             then parms_b also contains the following variables:             */
    /*                                                                             */
    /*                 P1_Intercept = intercept in the prob. part of the model     */
    /*                                for the first food/nutrient.                 */
    /*                 P1_varname   = regression slope for covariate "varname"     */
    /*                                in the prob. part of the model for the       */
    /*                                first food/nutrient.                         */
    /*                 Var_u1       = variance of random effect u1 (prob. part     */
    /*                                of model for first food/nutrient).           */
    /*                 Cov_u1u2     = covariance of random effects u1 and u2.      */
    /*                 Cov_u1u3     = covariance of random effects u1 and u3.      */
    /*                 Corr_u1u2    = correlation of random effects u1 and u2.     */
    /*                 Corr_u1u3    = correlation of random effects u1 and u3.     */
    /*                                                                             */
    /*             note: initial parameter estimates must be supplied by the       */
    /*                   user using the init_parms option.                         */
    /*                   the user-supplied data set will have the same variables   */
    /*                   as data set parms_b, except it should not include the     */
    /*                   following variables:                                      */
    /*                         z_e2e3                                              */
    /*                         Cov_e2e3                                            */
    /*                         Corr_e2e3                                           */
    /*                         Cov_u1u2 Cov_u1u3 Cov_u2u3                          */
    /*                         Corr_u1u3 Corr_u2u3                                 */
    /*                   All the necessary initial parameter estimates can be      */
    /*                   estimated using the SAS macro NLMixed_Univariate.         */
    /*                                                                             */
    /*   pred_x_b = data set containing predicted values for the model.            */
    /*              pred_x_b contains all the variables in the input data set,     */
    /*              plus the following variables:                                  */
    /*                                                                             */
    /*                 pred_x_a1 = predicted mean amount on consumption day for    */
    /*                             the first food/nutrient.                        */
    /*                 pred_x_a2 = predicted mean amount for the second            */
    /*                             food/nutrient.                                  */
    /*                                                                             */
    /*             if fitting the two-part model for the first food/nutrient,      */
    /*             then pred_x_b also contains the following variable:             */
    /*                                                                             */
    /*                 pred_x_p1 = predicted probability of consumption for        */
    /*                             the first food/nutrient.                        */
    /*                                                                             */
    /*******************************************************************************/

    *******************************************************************************************;
    *** Call macro NLMIXED_BIVARIATE to fit nonlinear mixed model in each stratum.  In this ***;
    *** example, we specify fixed values for parameters lambda1 (Box-Cox transformation for ***;
    *** dietary component 1) and lambda2 (Box-Cox transformation for dietary component 2).  ***;
    *** Alternatively, one can let the macro estimate lambda1 and lambda2 by not sepcifying ***;
    *** values for macro parameters "lambda1" and "lambda2".                                ***;
    *******************************************************************************************;

    title3 "Stratum 1 = Children 1-8";

    %nlmixed_bivariate(data          = nhanes1,
                       init_parms    = init_parms1,
                       subject       = id,
                       repeat        = repeat,
                       response1     = recall_dc1,
                       response2     = recall_dc2,
                       modeltype     = ONEPART,
                       covars_prob1  = ,
                       covars_amt1   = repeat2 weekend female age_4to8,
                       covars_amt2   = repeat2 weekend female age_4to8,
                       link          = ,
                       lambda1       = &lambda11,
                       lambda2       = &lambda12,
                       var_u1        = ,
                       var_u2        = ,
                       var_u3        = ,
                       corr_u1u2     = ,
                       corr_u1u3     = ,
                       corr_u2u3     = ,
                       replicate_var = &replicvar,
                       nloptions     = technique=trureg,
                       print         = Y,
                       ntitle        = 3
                       );

    data parms_b1;
      stratum = 1;
      set parms_b;
    run;

    data pred_b1;
      stratum = 1;
      set pred_x_b;
    run;



    title3 "Stratum 2 = Males 9+";

    %nlmixed_bivariate(data          = nhanes2,
                       init_parms    = init_parms2,
                       subject       = id,
                       repeat        = repeat,
                       response1     = recall_dc1,
                       response2     = recall_dc2,
                       modeltype     = ONEPART,
                       covars_prob1  = ,
                       covars_amt1   = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                       covars_amt2   = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                       link          = ,
                       lambda1       = &lambda21,
                       lambda2       = &lambda22,
                       var_u1        = ,
                       var_u2        = ,
                       var_u3        = ,
                       corr_u1u2     = ,
                       corr_u1u3     = ,
                       corr_u2u3     = ,
                       replicate_var = &replicvar,
                       nloptions     = technique=trureg,
                       print         = Y,
                       ntitle        = 3
                       );

    data parms_b2;
      stratum = 2;
      set parms_b;
    run;

    data pred_b2;
      stratum = 2;
      set pred_x_b;
    run;



    title3 "Stratum 3 = Females 9+";

    %nlmixed_bivariate(data           = nhanes3,
                       init_parms     = init_parms3,
                       subject        = id,
                       repeat         = repeat,
                       response1      = recall_dc1,
                       response2      = recall_dc2,
                       modeltype      = ONEPART,
                       covars_prob1   = ,
                       covars_amt1    = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                       covars_amt2    = repeat2 weekend age_14to18 age_19to30 age_31to50 age_51to70 age_71plus,
                       link           = ,
                       lambda1        = &lambda31,
                       lambda2        = &lambda32,
                       var_u1        = ,
                       var_u2        = ,
                       var_u3        = ,
                       corr_u1u2     = ,
                       corr_u1u3     = ,
                       corr_u2u3     = ,
                       replicate_var  = &replicvar,
                       nloptions      = technique=trureg,
                       print          = Y,
                       ntitle         = 3
                       );

    data parms_b3;
      stratum = 3;
      set parms_b;
    run;

    data pred_b3;
      stratum = 3;
      set pred_x_b;
    run;



    *********************************************************************;
    *** Combine data sets of parameter estimates and predicted values ***;
    *********************************************************************;

    data parms_b;
      set parms_b1 parms_b2 parms_b3;
        by stratum;
    run;

    data pred_b;
      set pred_b1 pred_b2 pred_b3;
        by stratum;
    run;



    **********************************************;
    *** Keep minimum amount on consumption day ***;
    **********************************************;

    data parms_b;
      merge parms_b min_a;
        by stratum;
    run;



    ********************************;
    *** Save parameter estimates ***;
    ********************************;

    data outlib.parms_b_sfat_energy&replicnum;
      set parms_b;
      run;

    data outlib.pred_b_sfat_energy&replicnum;
      set pred_b;
    run;

  %end;

%mend;



%fit_models_replicate_loop;
