options mergenoby=nowarn;

/*****************************************************************************/
/*****************************************************************************/
/* The following file is required to run this example:                       */
/*   predict_intake_density_macro_v1.2.sas                                   */
/*****************************************************************************/
/*                                                                           */
/*****************************************************************************/
/* This example analysis uses regression calibration and fits a logistic     */
/* regression model to assess the relationships between two dietary          */
/* components and a health outcome.  For this analysis, the example data     */
/* were simulated to include data from 100,000 participants, and the food    */
/* frequency questionnaire (FFQ) is the main dietary instrument.  The        */
/* dietary components are red meat and energy, and the health outcome is a   */
/* binary event variable.  The simulated data include a calibration substudy */
/* of 1,000 participants with repeated intake measurements from 24-hour      */
/* dietary recalls.                                                          */
/*                                                                           */
/* This program predicts the usual intake of the dietary components and fits */
/* a logistic model to assess the relationship between the dietary           */
/* components and health outcome.                                            */
/*                                                                           */
/* This analysis uses bootstrap variance estimation, so the bivariate        */
/* measurement error model and logistic model are fit using the original     */
/* data set and 200 replicate data sets.  The following replicfirst and      */
/* repliclast macro variables allow specification of a range of replicates.  */
/* This example program considers the original data set (i.e. replicate 0).  */
/*****************************************************************************/

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



title1 "Fit Univariate and Bivariate Measurement Error Models Using MLE with FFQ as Main Instrument";
title2 "Predict Intake and Perform Regression Calibration";
title3 "Assess Relationships between Two Dietary Components and a Health Outcome";



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

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

*** Include the required macros ***;
%include "&home/include.files.macros/predict_intake_density_macro_v1.2.sas";

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



*****************************************************;
*** Specify Box-Cox transformations in risk model ***;
*****************************************************;

%let lamt_td = 0;
%let lamt_t2 = 0;



*************************************************;
*** Comparison usual intake values: x0 and x1 ***;
*************************************************;

%let x0_td = 10;
%let x1_td = 60;

%let x0_t2 = 1500;
%let x1_t2 = 3000;



*********************************************************************;
*** Import the simulated data and bootstrap data for this example ***;
*********************************************************************;

data simdata_and_bootdata;
  set outlib.simdata_and_bootdata;
run;



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

%macro fit_models_replicate_loop;

  %do replicnum = &replicfirst %to &repliclast;

    ***************************************************************;
    *** Check if the bivariate parameter estimates exist.       ***;
    *** If not then exit the loop and begin the next iteration. ***;
    ***************************************************************;

    %if %sysfunc(exist(outlib._parms_b&replicnum)) %then %do;

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



      *********************************************************************;
      *** Create an input data set for the PREDICT_INTAKE_DENSITY macro ***;
      *********************************************************************;

      data parsubj1rec;

        if _n_ = 1 then set outlib._parms_b&replicnum;

        set select_replicate;

        lamt_td = &lamt_td;
        lamt_t2 = &lamt_t2;

        *****************************************************************;
        *** Calculate the linear predictor values using the parameter ***;
        *** estimates from NLMIXED_BIVARIATE and the covariates       ***;
        *****************************************************************;

        pred_t_xbeta1_prob = p1_intercept + p1_q1*q1 + p1_q2*q2;
        pred_t_xbeta2_amt = a1_intercept + a1_q1*q1 + a1_q2*q2;
        pred_t_xbeta3_amt = a2_intercept + a2_q1*q1 + a2_q2*q2;

        pred_x_p1_rep1 = pred_t_xbeta1_prob;
        pred_x_a1_rep1 = pred_t_xbeta2_amt;
        pred_x_a2_rep1 = pred_t_xbeta3_amt;

        pred_x_p1_rep2 = pred_t_xbeta1_prob;
        pred_x_a1_rep2 = pred_t_xbeta2_amt;
        pred_x_a2_rep2 = pred_t_xbeta3_amt;

      run;



      /*************************************************************************/
      /*************************************************************************/
      /*                                                                       */
      /*                     PREDICT_INTAKE_DENSITY MACRO                      */
      /*                                                                       */
      /*************************************************************************/
      /*************************************************************************/
      /*                                                                       */
      /*                                                                       */
      /* The PREDICT_INTAKE_DENSITY macro calculates predicted intakes or      */
      /* ratios of intakes (e.g. densities such as the ratio of the intake of  */
      /* interest to total energy intake) for regression calibration using     */
      /* results from one of the following models:                             */
      /* 1. amount-only model for a dietary component consumed daily,          */
      /* 2. two-part model for a dietary component consumed episodically,      */
      /* 3. bivariate model for two dietary components consumed daily,         */
      /* 4. bivariate model for a dietary component consumed episodically and  */
      /*    a dietary component consumed daily.                                */
      /*                                                                       */
      /* The macro performs adaptive Gaussian quadrature to predict usual      */
      /* intake, and the macro allows the user to provide a Box-Cox            */
      /* transformation parameter in order to calculate the predicted values   */
      /* on a transformed scale.  The results from this macro are intended for */
      /* use in a subsequent regression model as discussed by Kipnis et al.    */
      /* (2009, 2015).                                                         */
      /*                                                                       */
      /* References:                                                           */
      /*                                                                       */
      /* Kipnis V, Midthune D, Buckman DW, Dodd KW, Guenther PM, Krebs-Smith   */
      /* SM, Subar AF, Tooze JA, Carroll RJ, Freedman LS. Modeling data with   */
      /* excess zeros and measurement error: application to evaluating         */
      /* relationships between episodically consumed foods and health          */
      /* outcomes. Biometrics. 2009 Dec;65(4):1003-10.                         */
      /*                                                                       */
      /* Kipnis V, Freedman LS, Carroll RJ, Midthune D. A bivariate            */
      /* measurement error model for semicontinuous and continuous variables:  */
      /* Application to nutritional epidemiology. Biometrics. 2016 Mar;72(1):  */
      /* 106-15.                                                               */
      /*                                                                       */
      /*                                                                       */
      /* The syntax for calling the PREDICT_INTAKE_DENSITY macro is:           */
      /*                                                                       */
      /* %predict_intake_density(model14=, subj1recdata=, recid=,              */
      /*                         pred_t_dietcomp13=, pred_t_xbeta1_prob=,      */
      /*                         pred_t_xbeta2_amt=, pred_t_xbeta3_amt=,       */
      /*                         pred_t_min_amt2=, pred_t_min_amt3=,           */
      /*                         pred_t_constant=, boxcox_t_lamt=, lamt=,      */
      /*                         r24vars_dietcomp1=, r24vars_dietcomp2=,       */
      /*                         xbeta1vars=, xbeta2vars=, xbeta3vars=,        */
      /*                         var_u1=, var_u2=, var_u3=, cov_u1u2=,         */
      /*                         cov_u1u3=, cov_u2u3=, var_e2=, var_e3=,       */
      /*                         cov_e2e3=, link=, lambda2=, lambda3=,         */
      /*                         dencalc=, titles=, notesprt=);                */
      /*                                                                       */
      /* where the parameters are described as follows.                        */
      /*                                                                       */
      /*  "model14"            Specifies the type of model that was fit prior  */
      /*                       to calling this macro.  The values 1, 2, 3, or  */
      /*                       4 indicate which of the following models was    */
      /*                       fit:                                            */
      /*                       1 = amount-only model for a dietary component   */
      /*                           consumed daily,                             */
      /*                       2 = two-part model for a dietary component      */
      /*                           consumed episodically,                      */
      /*                       3 = bivariate model for two dietary components  */
      /*                           consumed daily,                             */
      /*                       4 = bivariate model for a dietary component     */
      /*                           consumed episodically and a dietary         */
      /*                           component consumed daily.                   */
      /*                                                                       */
      /*  "subj1recdata"       Specifies a data set with 1 record for each     */
      /*                       individual.                                     */
      /*                                                                       */
      /*  "recid"              Specifies an identification (ID) variable that  */
      /*                       uniquely identifies each individual's record.   */
      /*                                                                       */
      /*  "pred_t_dietcomp13"  Specifies whether the macro predicts intake for */
      /*                       dietary component 1, dietary component 2, or    */
      /*                       the ratio of dietary component 1 to dietary     */
      /*                       component 2.                                    */
      /*                                                                       */
      /*  "pred_t_xbeta1_prob" Specifies a variable that provides the linear   */
      /*                       predictor values used to define the probability */
      /*                       part of the function predicted by the macro.    */
      /*                       The linear predictor values are calculated      */
      /*                       using the covariates and estimates of the fixed */
      /*                       effects parameters from the probability part of */
      /*                       the model.  For related notational details, see */
      /*                       the description for the "xbeta1vars" macro      */
      /*                       parameter.                                      */
      /*                                                                       */
      /*  "pred_t_xbeta2_amt"  Specifies a variable that provides the linear   */
      /*                       predictor values used to define the amount part */
      /*                       of the function predicted by the macro.  The    */
      /*                       linear predictor values are calculated using    */
      /*                       the covariates and estimates of the fixed       */
      /*                       effects parameters from the amount part of the  */
      /*                       model.  For related notational details, see the */
      /*                       description for the "xbeta2vars" macro          */
      /*                       parameter.                                      */
      /*                                                                       */
      /*  "pred_t_xbeta3_amt"  Specifies a variable that provides the linear   */
      /*                       predictor values used to define the amount part */
      /*                       of the function predicted by the macro.  The    */
      /*                       linear predictor values are calculated using    */
      /*                       the covariates and estimates of the fixed       */
      /*                       effects parameters from the amount part of the  */
      /*                       model.  For related notational details, see the */
      /*                       description for the "xbeta3vars" macro          */
      /*                       parameter.                                      */
      /*                                                                       */
      /*  "pred_t_min_amt2"    Specifies a variable that provides the minimum  */
      /*                       intake amount used to define the function       */
      /*                       predicted by the macro.  This value can be      */
      /*                       obtained by considering all of the available    */
      /*                       consumption-day amounts for the first dietary   */
      /*                       component and using the smallest value, and the */
      /*                       specified variable provides this value for each */
      /*                       individual.  This value will be divided in half */
      /*                       and used in the calculations.                   */
      /*                                                                       */
      /*  "pred_t_min_amt3"    Specifies a variable that provides the minimum  */
      /*                       intake amount used to define the function       */
      /*                       predicted by the macro.  This value can be      */
      /*                       obtained by considering all of the available    */
      /*                       consumption-day amounts for the second dietary  */
      /*                       component and using the smallest value, and the */
      /*                       specified variable provides this value for each */
      /*                       individual.  This value will be divided in half */
      /*                       and used in the calculations.                   */
      /*                                                                       */
      /*  "pred_t_constant"    Specifies a constant used in the calculation of */
      /*                       predicted densities.  The default value is 1.   */
      /*                                                                       */
      /*  "boxcox_t_lamt"      If "boxcox_t_lamt=y" or "boxcox_t_lamt=Y" then  */
      /*                       intake will be predicted on a transformed scale */
      /*                       where the Box-Cox transformation is used with   */
      /*                       the Box-Cox parameter value provided by the     */
      /*                       "lamt" macro parameter.  The default value for  */
      /*                       "boxcox_t_lamt" is "n".                         */
      /*                                                                       */
      /*  "lamt"               Specifies a variable that provides the Box-Cox  */
      /*                       parameter value when "boxcox_t_lamt=y" or       */
      /*                       "boxcox_t_lamt=Y".  The macro does not allow    */
      /*                       the Box-Cox parameter to be negative.           */
      /*                                                                       */
      /*  "r24vars_dietcomp1"  Specifies the 24-hour recall variables for      */
      /*                       dietary component 1.  Each variable value is    */
      /*                       either non-negative or a SAS missing value.     */
      /*                       Variables must be space delimited as            */
      /*                       illustrated in the following example:           */
      /*                       "r24vars=r24hr1 r24hr2".                        */
      /*                       Note for Advanced Users:  If all 24-hour recall */
      /*                       values are missing for each subject, then the   */
      /*                       denominator integration should not be           */
      /*                       performed, so the "dencalc" macro parameter     */
      /*                       should be specified as "dencalc=n" and no       */
      /*                       variables are provided for the                  */
      /*                       "r24vars_dietcomp1" and "r24vars_dietcomp2"     */
      /*                       macro parameters.                               */
      /*                                                                       */
      /*  "r24vars_dietcomp2"  Specifies the 24-hour recall variables for      */
      /*                       dietary component 2.  For details, see the      */
      /*                       description for the "r24vars_dietcomp1" macro   */
      /*                       parameter.                                      */
      /*                                                                       */
      /*  "xbeta1vars"         Specifies a space delimited list of variables   */
      /*                       that provide the linear predictor values        */
      /*                       calculated using the covariates for each        */
      /*                       24-hour recall period and estimates of the      */
      /*                       fixed effects parameters from the probability   */
      /*                       part of the model for dietary component 1.      */
      /*                                                                       */
      /*  "xbeta2vars"         Specifies a space delimited list of variables   */
      /*                       that provide the linear predictor values        */
      /*                       calculated using the covariates for each        */
      /*                       24-hour recall period and estimates of the      */
      /*                       fixed effects parameters from the amount part   */
      /*                       of the model for dietary component 1.           */
      /*                                                                       */
      /*  "xbeta3vars"         Specifies a space delimited list of variables   */
      /*                       that provide the linear predictor values        */
      /*                       calculated using the covariates for each        */
      /*                       24-hour recall period and estimates of the      */
      /*                       fixed effects parameters from the amount part   */
      /*                       of the model for dietary component 2.           */
      /*                                                                       */
      /*  "var_u1"             Specifies a variable that provides the variance */
      /*                       estimate for u1, the random effect from the     */
      /*                       probability part of the model for dietary       */
      /*                       component 1.                                    */
      /*                                                                       */
      /*  "var_u2"             Specifies a variable that provides the variance */
      /*                       estimate for u2, the random effect from the     */
      /*                       amount part of the model for dietary component  */
      /*                       1.                                              */
      /*                                                                       */
      /*  "var_u3"             Specifies a variable that provides the variance */
      /*                       estimate for u3, the random effect from the     */
      /*                       amount part of the model for dietary component  */
      /*                       2.                                              */
      /*                                                                       */
      /*  "cov_u1u2"           Specifies a variable that provides the estimate */
      /*                       of the covariance(u1, u2).                      */
      /*                                                                       */
      /*  "cov_u1u3"           Specifies a variable that provides the estimate */
      /*                       of the covariance(u1, u3).                      */
      /*                                                                       */
      /*  "cov_u2u3"           Specifies a variable that provides the estimate */
      /*                       of the covariance(u2, u3).                      */
      /*                                                                       */
      /*  "var_e2"             Specifies a variable that provides the variance */
      /*                       estimate for e2, the within-person error term   */
      /*                       from the amount part of the model for dietary   */
      /*                       component 1.                                    */
      /*                                                                       */
      /*  "var_e3"             Specifies a variable that provides the variance */
      /*                       estimate for e3, the within-person error term   */
      /*                       from the amount part of the model for dietary   */
      /*                       component 2.                                    */
      /*                                                                       */
      /*  "corr_e1e3"          Specifies a variable that provides the estimate */
      /*                       of the correlation(e1, e3).                     */
      /*                                                                       */
      /*  "corr_e2e3"          Specifies a variable that provides the estimate */
      /*                       of the correlation(e2, e3).                     */
      /*                                                                       */
      /*  "link"               Specifies the link function for the probability */
      /*                       part of the model for dietary component 1.      */
      /*                         If "link = probit" a probit link function is  */
      /*                         used.                                         */
      /*                         If "link = logit" a logit link function is    */
      /*                         used.                                         */
      /*                       The default value for "link" is "probit".       */
      /*                                                                       */
      /*  "lambda2"            Specifies a variable that provides the estimate */
      /*                       of the Box-Cox parameter, lambda2, from the     */
      /*                       amount part of the model for dietary component  */
      /*                       1.                                              */
      /*                       The macro does not allow the Box-Cox parameter  */
      /*                       to be zero.                                     */
      /*                                                                       */
      /*  "lambda3"            Specifies a variable that provides the estimate */
      /*                       of the Box-Cox parameter, lambda3, from the     */
      /*                       amount part of the model for dietary component  */
      /*                       2.                                              */
      /*                       The macro does not allow the Box-Cox parameter  */
      /*                       to be zero.                                     */
      /*                                                                       */
      /*  "dencalc"            By default, "dencalc=y" so the denominator      */
      /*                       integration is performed.                       */
      /*                       Note for Advanced Users:  If all 24-hour recall */
      /*                       variables are missing for each subject, then    */
      /*                       the denominator integration should not be       */
      /*                       performed, so the "dencalc" option should be    */
      /*                       specified as "dencalc=n".                       */
      /*                                                                       */
      /*  "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".                                   */
      /*                                                                       */
      /*  "notesprt"           If "notesprt=n" or "notesprt=N" then notes are  */
      /*                       not printed to the SAS log.  The default value  */
      /*                       for "notesprt" is "y".                          */
      /*                                                                       */
      /*************************************************************************/

      ************************************************************************;
      *** Call the PREDICT_INTAKE_DENSITY macro to calculate the predicted ***;
      *** intake values for the ratio of red meat to energy                ***;
      ************************************************************************;

      %predict_intake_density(model14            = 4,
                              subj1recdata       = parsubj1rec,
                              recid              = replicaterowid,
                              pred_t_dietcomp13  = 3,
                              pred_t_xbeta1_prob = pred_t_xbeta1_prob,
                              pred_t_xbeta2_amt  = pred_t_xbeta2_amt,
                              pred_t_xbeta3_amt  = pred_t_xbeta3_amt,
                              pred_t_min_amt2    = min_a1,
                              pred_t_min_amt3    = min_a2,
                              pred_t_constant    = 1000,
                              boxcox_t_lamt      = y,
                              lamt               = lamt_td,
                              r24vars_dietcomp1  = r11 r12,
                              r24vars_dietcomp2  = r21 r22,
                              xbeta1vars         = pred_x_p1_rep1 pred_x_p1_rep2,
                              xbeta2vars         = pred_x_a1_rep1 pred_x_a1_rep2,
                              xbeta3vars         = pred_x_a2_rep1 pred_x_a2_rep2,
                              var_u1             = var_u1,
                              var_u2             = var_u2,
                              var_u3             = var_u3,
                              cov_u1u2           = cov_u1u2,
                              cov_u1u3           = cov_u1u3,
                              cov_u2u3           = cov_u2u3,
                              var_e2             = var_e2,
                              var_e3             = var_e3,
                              corr_e1e3          = corr_e1e3,
                              corr_e2e3          = corr_e2e3,
                              link               = probit,
                              lambda2            = a1_lambda,
                              lambda3            = a2_lambda,
                              dencalc            = y,
                              titles             = 4,
                              notesprt           = y
                              );



      *****************************************************************;
      *** Create a data set with the PREDICT_INTAKE_DENSITY results ***;
      *****************************************************************;

      data subj1recres;
        merge select_replicate _resdata(keep = replicaterowid indusint);
          by replicaterowid;
        rename indusint = predint_td;
      run;



      **************************************************************;
      *** Call the PREDICT_INTAKE_DENSITY macro to calculate the ***;
      *** predicted intake values for dietary component 2        ***;
      **************************************************************;

      %predict_intake_density(model14            = 4,
                              subj1recdata       = parsubj1rec,
                              recid              = replicaterowid,
                              pred_t_dietcomp13  = 2,
                              pred_t_xbeta1_prob = ,
                              pred_t_xbeta2_amt  = ,
                              pred_t_xbeta3_amt  = pred_t_xbeta3_amt,
                              pred_t_min_amt2    = ,
                              pred_t_min_amt3    = min_a2,
                              pred_t_constant    = ,
                              boxcox_t_lamt      = y,
                              lamt               = lamt_t2,
                              r24vars_dietcomp1  = r11 r12,
                              r24vars_dietcomp2  = r21 r22,
                              xbeta1vars         = pred_x_p1_rep1 pred_x_p1_rep2,
                              xbeta2vars         = pred_x_a1_rep1 pred_x_a1_rep2,
                              xbeta3vars         = pred_x_a2_rep1 pred_x_a2_rep2,
                              var_u1             = var_u1,
                              var_u2             = var_u2,
                              var_u3             = var_u3,
                              cov_u1u2           = cov_u1u2,
                              cov_u1u3           = cov_u1u3,
                              cov_u2u3           = cov_u2u3,
                              var_e2             = var_e2,
                              var_e3             = var_e3,
                              corr_e1e3          = corr_e1e3,
                              corr_e2e3          = corr_e2e3,
                              link               = probit,
                              lambda2            = a1_lambda,
                              lambda3            = a2_lambda,
                              dencalc            = y,
                              titles             = 4,
                              notesprt           = y
                              );



      *****************************************************************;
      *** Create a data set with the PREDICT_INTAKE_DENSITY results ***;
      *****************************************************************;

      data subj1recres;
        merge subj1recres _resdata(keep = replicaterowid indusint);
          by replicaterowid;
        rename indusint = predint_t2;
      run;



      *****************************************************************************************************;
      *** Standardize the covariates so that one unit equals the distance between transformed x0 and x1 ***;
      *****************************************************************************************************;

      data subj1recres;
        if _n_ = 1 then set parsubj1rec(obs=1 keep = lamt_td lamt_t2);
        set subj1recres;
        
        if lamt_td  = 0 then do;
          x0_td_bc = log(&x0_td);
          x1_td_bc = log(&x1_td);
        end;
        else do;
          x0_td_bc = (&x0_td**lamt_td - 1) / lamt_td;
          x1_td_bc = (&x1_td**lamt_td - 1) / lamt_td;
        end;

        unit_td = x1_td_bc - x0_td_bc;
        predint_td_unit = predint_td / unit_td;


        if lamt_t2  = 0 then do;
          x0_t2_bc = log(&x0_t2);
          x1_t2_bc = log(&x1_t2);
        end;
        else do;
          x0_t2_bc = (&x0_t2**lamt_t2 - 1) / lamt_t2;
          x1_t2_bc = (&x1_t2**lamt_t2 - 1) / lamt_t2;
        end;

        unit_t2 = x1_t2_bc - x0_t2_bc;
        predint_t2_unit = predint_t2 / unit_t2;
      run;



      *******************************************************************;
      *** Fit a logistic regression risk model for the health outcome ***;
      *******************************************************************;

      proc logistic data=subj1recres outest=regout;
        model y(event='1') = predint_td_unit predint_t2_unit;
        title5 "Logistic Regression Risk Model: Adjusted for Measurement Error";
      run;
      title5;



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

      data outlib.diethealthout&replicnum(keep = replicate intercept rc_beta_td_unit rc_or_td_unit
                                                 rc_beta_t2_unit rc_or_t2_unit x0_td x1_td x0_t2 x1_t2);

        set regout(rename = (predint_td_unit = rc_beta_td_unit
                             predint_t2_unit = rc_beta_t2_unit));

        replicate = &replicnum;



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

        rc_or_td_unit = exp(rc_beta_td_unit);
        rc_or_t2_unit = exp(rc_beta_t2_unit);



        ****************************************************************************;
        *** Save the comparison values, x0 and x1, for each risk model covariate ***;
        ****************************************************************************;

        x0_td = &x0_td;
        x1_td = &x1_td;

        x0_t2 = &x0_t2;
        x1_t2 = &x1_t2;

      run;

    %end;



    %else %do;

      **********************************************************************************************;
      *** Create an empty results output data set if this iteration is skipped.                  ***;
      *** These empty data sets will help when summarizing the results from multiple replicates. ***;
      **********************************************************************************************;

      data outlib.diethealthout&replicnum;
        set _null_;
      run;

    %end;

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

%mend fit_models_replicate_loop;



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

%fit_models_replicate_loop;
