/*************************************************************************/
/*************************************************************************/
/*                                                                       */
/*                     PREDICT_INTAKE_DENSITY MACRO                      */
/*                                                                       */
/*************************************************************************/
/*************************************************************************/
/*                  VERSION 1.2               2/1/2018                   */
/*                                                                       */
/*                                                                       */
/* 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=,       */
/*                         corr_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".                          */
/*                                                                       */
/*************************************************************************/


%macro 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=1,
                              boxcox_t_lamt=n,
                              lamt=,
                              r24vars_dietcomp1=,
                              r24vars_dietcomp2=,
                              xbeta1vars=,
                              xbeta2vars=,
                              xbeta3vars=,
                              var_u1=,
                              var_u2=,
                              var_u3=,
                              cov_u1u2=,
                              cov_u1u3=,
                              cov_u2u3=,
                              var_e2=,
                              var_e3=,
                              corr_e1e3=,
                              corr_e2e3=,
                              link=probit,
                              lambda2=,
                              lambda3=,
                              dencalc=y,
                              titles=0,
                              notesprt=y);




  /**********************************************************************************/
  /*** Ensure user macro call does not overwrite default values with a null value ***/
  /**********************************************************************************/

  %if (&pred_t_constant=%str()) %then %let pred_t_constant = 1;
  %if (&boxcox_t_lamt=%str()) %then %let boxcox_t_lamt = n;
  %if (&link=%str()) %then %let link = probit;
  %if (&dencalc=%str()) %then %let dencalc = y;
  %if (&titles=%str()) %then %let titles = 0;
  %if (&notesprt=%str()) %then %let notesprt = y;


  /*********************************************************************/
  /***  Check user input and initialize macro options and variables  ***/
  /*********************************************************************/

  %if (&model14 ^= 1) & (&model14 ^= 2) & (&model14 ^= 3) & (&model14 ^= 4) %then %put Error:  Invalid value for model14.;

  %if (&pred_t_dietcomp13 ^= 1) & (&pred_t_dietcomp13 ^= 2) & (&pred_t_dietcomp13 ^= 3) %then %put Error:  Invalid value for pred_t_dietcomp13.;

  %if (%upcase(&dencalc)^=Y) & (%upcase(&dencalc)^=N) %then %put Error:  Invalid value for dencalc.;

  %if (%upcase(&dencalc)=N) & (&r24vars_dietcomp1 ^= %str()) %then %put Error:  Incompatible input for dencalc and r24vars_dietcomp1.;

  %if ((&model14 = 1) | (&model14 = 3)) %then %do;
    %if (&pred_t_xbeta1_prob ^= %str()) %then %put Error:  Incompatible input for model14 and pred_t_xbeta1_prob.;
    %if (&xbeta1vars ^= %str()) %then %put Error:  Incompatible input for model14 and xbeta1vars.;
    %if (&var_u1 ^= %str()) %then %put Error:  Incompatible input for model14 and var_u1.;
  %end;
  %else %if ((&model14 = 2) | (&model14 = 4)) %then %do;
    %if ((&pred_t_dietcomp13 = 1) | (&pred_t_dietcomp13 = 3)) & (&pred_t_xbeta1_prob = %str()) %then %put Error:  Incompatible input for model14, pred_t_dietcomp13, and pred_t_xbeta1_prob.;
    %else %if (&pred_t_dietcomp13 = 2) & (&pred_t_xbeta1_prob ^= %str()) %then %put Error:  Incompatible input for model14, pred_t_dietcomp13, and pred_t_xbeta1_prob.;
  %end;


  %if (&var_e2 = %str()) %then %put Error:  Specify value for var_e2.;
  %if (&lambda2 = %str()) %then %put Error:  Specify value for lambda2.;


  %if ((&model14 = 1) | (&model14 = 2)) %then %do;
    %if ((&pred_t_dietcomp13 = 2) | (&pred_t_dietcomp13 = 3)) %then %put Error:  Incompatible input for model14 and pred_t_dietcomp13.;
    %if (&r24vars_dietcomp2 ^= %str()) %then %put Error:  Incompatible input for model14 and r24vars_dietcomp2.;
    %if (&xbeta3vars ^= %str()) %then %put Error:  Incompatible input for model14 and xbeta3vars.;
    %if (&var_u3 ^= %str()) %then %put Error:  Incompatible input for model14 and var_u3.;
    %if (&var_e3 ^= %str()) %then %put Error:  Incompatible input for model14 and var_e3.;
    %if (&corr_e1e3 ^= %str()) %then %put Error:  Incompatible input for model14 and corr_e1e3.;
    %if (&corr_e2e3 ^= %str()) %then %put Error:  Incompatible input for model14 and corr_e2e3.;
    %if (&lambda3 ^= %str()) %then %put Error:  Incompatible input for model14 and lambda3.;
  %end;
  %else %if ((&model14=3) | (&model14=4)) %then %do;
    %if (%upcase(&dencalc)=N) & (&r24vars_dietcomp2 ^= %str()) %then %put Error:  Incompatible input for model14, dencalc, and r24vars_dietcomp2.;
    %if (&var_e3 = %str()) %then %put Error:  Incompatible input for model14 and var_e3.;
    %if (&corr_e2e3 = %str()) %then %put Error:  Incompatible input for model14 and corr_e2e3.;
    %if (&lambda3 = %str()) %then %put Error:  Incompatible input for model14 and lambda3.;
    %if (&model14 = 4 & %upcase(&link)=PROBIT) %then %do;
      %if (&corr_e1e3 = %str()) %then %put Error:  Incompatible input for model14 and corr_e1e3.;
    %end;
  %end;


  %if %upcase(&boxcox_t_lamt)=Y %then %do;
    %if (&lamt = %str()) %then %put Error:  Incompatible input for boxcox_t_lamt and lamt.;
  %end;
  %else %if %upcase(&boxcox_t_lamt)=N %then %do;
    %if (&lamt ^= %str()) %then %put Error:  Incompatible input for boxcox_t_lamt and lamt.;
  %end;
  %else %put Error:  Invalid value for boxcox_t_lamt.;

  %if (&subj1recdata = %str()) %then %put Error:  Specify value for subj1recdata.;
  %if (&recid = %str()) %then %put Error:  Specify value for recid.;

  %if %upcase(&notesprt)=Y %then %do;
    options notes;
  %end;
  %else %if %upcase(&notesprt)=N %then %do;
    options nonotes;
  %end;
  %else %put Error:  Unknown value for notesprt.;


  %if (&var_u1 ^= %str()) & (&var_u2 ^= %str()) & (&cov_u1u2 = %str()) %then %put Error:  Incompatible input for var_u1, var_u2, and cov_u1u2.;
  %else %if ((&var_u1 = %str()) | (&var_u2 = %str())) & (&cov_u1u2 ^= %str()) %then %put Error:  Incompatible input for var_u1, var_u2, and cov_u1u2.;

  %if (&var_u1 ^= %str()) & (&var_u3 ^= %str()) & (&cov_u1u3 = %str()) %then %put Error:  Incompatible input for var_u1, var_u3, and cov_u1u3.;
  %else %if ((&var_u1 = %str()) | (&var_u3 = %str())) & (&cov_u1u3 ^= %str()) %then %put Error:  Incompatible input for var_u1, var_u3, and cov_u1u3.;

  %if (&var_u2 ^= %str()) & (&var_u3 ^= %str()) & (&cov_u2u3 = %str()) %then %put Error:  Incompatible input for var_u2, var_u3, and cov_u2u3.;
  %else %if ((&var_u2 = %str()) | (&var_u3 = %str())) & (&cov_u2u3 ^= %str()) %then %put Error:  Incompatible input for var_u2, var_u3, and cov_u2u3.;


  /**********************************************/
  /*** Random effects macro variables for IML ***/
  /**********************************************/

  %if (&var_u1 = %str()) & (&var_u2 = %str()) & (&var_u3 = %str()) %then %do;

    %put Error:  Variance of the random effect(s) was not provided.;

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 = %str()) & (&var_u2 = %str()) & (&var_u3 ^= %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 1;

    %let imlu1 = 0;
    %let imlu2 = 0;
    %let imlu3 = u[1];

    /*** Covariance matrix for random effect ***/

    %let imlcovmatu = %str(
                           covmatu = &var_u3;
                           );

    %let dimcovu = 1;

    %let imlufinaln = %str(
      u1finaln = .;
      u2finaln = .;
      u3finaln = ufinaln[1];
      );

    %let imlufinald = %str(
      u1finald = .;
      u2finald = .;
      u3finald = ufinald[1];
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 = %str()) & (&var_u2 ^= %str()) & (&var_u3 = %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 2;

    %let imlu1 = 0;
    %let imlu2 = u[1];
    %let imlu3 = 0;

    /*** Covariance matrix for random effect ***/

    %let imlcovmatu = %str(
                           covmatu = &var_u2;
                           );

    %let dimcovu = 1;

    %let imlufinaln = %str(
      u1finaln = .;
      u2finaln = ufinaln[1];
      u3finaln = .;
      );

    %let imlufinald = %str(
      u1finald = .;
      u2finald = ufinald[1];
      u3finald = .;
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 = %str()) & (&var_u2 ^= %str()) & (&var_u3 ^= %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 3;

    %let imlu1 = 0;
    %let imlu2 = u[1];
    %let imlu3 = u[2];

    /*** Covariance matrix for random effects ***/

    %let imlcovmatu = %str(
      row1    = &var_u2 || &cov_u2u3;
      row2    = &cov_u2u3 || &var_u3;
      covmatu = row1 // row2;
      );

    %let dimcovu = 2;

    %let imlufinaln = %str(
      u1finaln = .;
      u2finaln = ufinaln[1];
      u3finaln = ufinaln[2];
      );

    %let imlufinald = %str(
      u1finald = .;
      u2finald = ufinald[1];
      u3finald = ufinald[2];
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 ^= %str()) & (&var_u2 = %str()) & (&var_u3 = %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 4;

    %let imlu1 = u[1];
    %let imlu2 = 0;
    %let imlu3 = 0;

    /*** Covariance matrix for random effect ***/

    %let imlcovmatu = %str(
                           covmatu = &var_u1;
                           );

    %let dimcovu = 1;

    %let imlufinaln = %str(
      u1finaln = ufinaln[1];
      u2finaln = .;
      u3finaln = .;
      );

    %let imlufinald = %str(
      u1finald = ufinald[1];
      u2finald = .;
      u3finald = .;
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 ^= %str()) & (&var_u2 = %str()) & (&var_u3 ^= %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 5;

    %let imlu1 = u[1];
    %let imlu2 = 0;
    %let imlu3 = u[2];

    /*** Covariance matrix for random effects ***/

    %let imlcovmatu = %str(
      row1    = &var_u1 || &cov_u1u3;
      row2    = &cov_u1u3 || &var_u3;
      covmatu = row1 // row2;
      );

    %let dimcovu = 2;

    %let imlufinaln = %str(
      u1finaln = ufinaln[1];
      u2finaln = .;
      u3finaln = ufinaln[2];
      );

    %let imlufinald = %str(
      u1finald = ufinald[1];
      u2finald = .;
      u3finald = ufinald[2];
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 ^= %str()) & (&var_u2 ^= %str()) & (&var_u3 = %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 6;

    %let imlu1 = u[1];
    %let imlu2 = u[2];
    %let imlu3 = 0;

    /*** Covariance matrix for random effects ***/

    %let imlcovmatu = %str(
      row1    = &var_u1 || &cov_u1u2;
      row2    = &cov_u1u2 || &var_u2;
      covmatu = row1 // row2;
      );

    %let dimcovu = 2;

    %let imlufinaln = %str(
      u1finaln = ufinaln[1];
      u2finaln = ufinaln[2];
      u3finaln = .;
      );

    %let imlufinald = %str(
      u1finald = ufinald[1];
      u2finald = ufinald[2];
      u3finald = .;
      );

  %end;
  /***********************************************************************************/
  %else %if (&var_u1 ^= %str()) & (&var_u2 ^= %str()) & (&var_u3 ^= %str()) %then %do;

    /***************************************/
    /*** Random effects code used in IML ***/
    /***************************************/

    %let randeffset = 7;

    %let imlu1 = u[1];
    %let imlu2 = u[2];
    %let imlu3 = u[3];

    /*** Covariance matrix for random effects ***/

    %let imlcovmatu = %str(
      row1    = &var_u1   || &cov_u1u2 || &cov_u1u3;
      row2    = &cov_u1u2 || &var_u2   || &cov_u2u3;
      row3    = &cov_u1u3 || &cov_u2u3 || &var_u3;
      covmatu = row1 // row2 // row3;
      );

    %let dimcovu = 3;

    %let imlufinaln = %str(
      u1finaln = ufinaln[1];
      u2finaln = ufinaln[2];
      u3finaln = ufinaln[3];
      );

    %let imlufinald = %str(
      u1finald = ufinald[1];
      u2finald = ufinald[2];
      u3finald = ufinald[3];
      );

  %end;




  /********************************************************************/
  /*** Covariance matrix for e2 or (e2,e3), macro variables for IML ***/
  /********************************************************************/
  %if (&model14 = 1) | (&model14 = 2) %then %do;

    /*** Covariance matrix for e2 ***/

    %let imlcovmate = %str(
                           covmate = &var_e2;
                           );

  %end;
  /********************************************************************/
  %else %if (&model14 = 3) | (&model14 = 4) %then %do;

    %let cov_e2e3 = cov_e2e3;

    /*** Covariance matrix for (e2,e3) ***/

    %let imlcovmate = %str(
      row1    = &var_e2   || &cov_e2e3;
      row2    = &cov_e2e3 || &var_e3;
      covmate = row1 // row2;
      );

  %end;
  /********************************************************************/






  data _null_;
    set &subj1recdata(obs=1);

    %if &lamt^=%str() %then %do;
      if &lamt<0 then put "Error:  Macro does not allow the lamt variable to have a negative value.";
    %end;

    /*** Determine the number of repeats and check input variables ***/

    %if &r24vars_dietcomp1=%str() %then %do;
      dc1numr24=0;
    %end;
    %else %do;
      array dc1rvec(*)  &r24vars_dietcomp1;
      dc1numr24 = dim(dc1rvec);
    %end;
    num_repeats = dc1numr24;

    %if &r24vars_dietcomp2=%str() %then %do;
      dc2numr24=0;
    %end;
    %else %do;
      array dc2rvec(*)  &r24vars_dietcomp2;
      dc2numr24 = dim(dc2rvec);
    %end;

    %if &xbeta1vars=%str() %then %do;
      numxbeta1=0;
    %end;
    %else %do;
      array xb1vec(*)  &xbeta1vars;
      numxbeta1 = dim(xb1vec);
    %end;

    %if &xbeta2vars=%str() %then %do;
      numxbeta2=0;
    %end;
    %else %do;
      array xb2vec(*)  &xbeta2vars;
      numxbeta2 = dim(xb2vec);
    %end;

    %if &xbeta3vars=%str() %then %do;
      numxbeta3=0;
    %end;
    %else %do;
      array xb3vec(*)  &xbeta3vars;
      numxbeta3 = dim(xb3vec);
    %end;


    %if (&model14 = 1) %then %do;
      if (numxbeta2 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta2.";
        num_repeats = .;
      end;
    %end;
    %else %if (&model14 = 2) %then %do;
      if (numxbeta1 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta1.";
        num_repeats = .;
      end;
      else if (numxbeta2 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta2.";
        num_repeats = .;
      end;
    %end;
    %else %if (&model14 = 3) %then %do;
      if (dc2numr24 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary components 1 and 2.";
        num_repeats = .;
      end;
      else if (numxbeta2 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta2.";
        num_repeats = .;
      end;
      else if (numxbeta3 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta3.";
        num_repeats = .;
      end;
    %end;
    %else %if (&model14 = 4) %then %do;
      if (dc2numr24 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary components 1 and 2.";
        num_repeats = .;
      end;
      else if (numxbeta1 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta1.";
        num_repeats = .;
      end;
      else if (numxbeta2 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta2.";
        num_repeats = .;
      end;
      else if (numxbeta3 ^= num_repeats) then do;
        put "Error:  Number of repeat variables unequal for 24-hour recall, dietary component 1, and xbeta3.";
        num_repeats = .;
      end;
    %end;

    call symput("num_repeatvars",trim(left(put(num_repeats, 3.))));
    call symput("numr24vars_dietcomp1",trim(left(put(dc1numr24, 3.))));
    call symput("numr24vars_dietcomp2",trim(left(put(dc2numr24, 3.))));
    call symput("numxbeta1vars",trim(left(put(numxbeta1, 3.))));
    call symput("numxbeta2vars",trim(left(put(numxbeta2, 3.))));
    call symput("numxbeta3vars",trim(left(put(numxbeta3, 3.))));

  run;

  %put Note:  Number of 24-hour recall variables for dietary component 1 is &numr24vars_dietcomp1;
  %put Note:  Number of 24-hour recall variables for dietary component 2 is &numr24vars_dietcomp2;
  %put Note:  Number of xbeta1 variables is &numxbeta1vars;
  %put Note:  Number of xbeta2 variables is &numxbeta2vars;
  %put Note:  Number of xbeta3 variables is &numxbeta3vars;




  /*******************************************************************************************/
  /*******************************************************************************************/


  proc format;
    value r24errzf
      0         = "Error: 0 Invalid"
      0< - high = "GT Zero"
      ;
    value r24gtzf
      0< - high = "GT Zero"
      ;
  run;


  data _subj1rec(drop =
                    %if &r24vars_dietcomp1^=%str() %then %str( kr24vars_dietcomp1 );
                    %if &r24vars_dietcomp2^=%str() %then %str( kr24vars_dietcomp2 );
                    dc1miss24message dc2miss24message);
    set &subj1recdata end=lastrec;
    by &recid;

    retain dc1miss24message dc2miss24message 1;


    %if %upcase(&boxcox_t_lamt)=N %then %do;
      %let lamt=lamt;
      &lamt=1;  /*** No transformation.  Therefore, lamt must be 1 because T##lamt/lamt is used. ***/
    %end;



    %if ((&model14=3) | (&model14=4)) %then %do;
      cov_e2e3 = sqrt(&var_e2) * sqrt(&var_e3) * &corr_e2e3;
    %end;



    %if &r24vars_dietcomp1^=%str() %then %do;

      array dc1rvecin(&num_repeatvars)  &r24vars_dietcomp1;
      array dc1rvec(&num_repeatvars)    dc1r1-dc1r&num_repeatvars;

      do kr24vars_dietcomp1 = 1 to &num_repeatvars;
        dc1rvec(kr24vars_dietcomp1) = dc1rvecin(kr24vars_dietcomp1);
        if dc1rvecin(kr24vars_dietcomp1) >= 0 then dc1miss24message=0;  /*** Set flag so missing value warning message is not printed ***/
      end;

    %end;

    %if &r24vars_dietcomp2^=%str() %then %do;

      array dc2rvecin(&num_repeatvars)  &r24vars_dietcomp2;
      array dc2rvec(&num_repeatvars)    dc2r1-dc2r&num_repeatvars;

      do kr24vars_dietcomp2 = 1 to &num_repeatvars;
        dc2rvec(kr24vars_dietcomp2) = dc2rvecin(kr24vars_dietcomp2);
        if dc2rvecin(kr24vars_dietcomp2) >= 0 then dc2miss24message=0;  /*** Set flag so missing value warning message is not printed ***/
      end;

    %end;

    %if &xbeta1vars^=%str() %then %do;

      array xb1vecin(&num_repeatvars)  &xbeta1vars;
      array xb1vec(&num_repeatvars)    xbeta1_rep1-xbeta1_rep&num_repeatvars;

      do kxbeta1vars = 1 to &num_repeatvars;
        xb1vec(kxbeta1vars) = xb1vecin(kxbeta1vars);
      end;

    %end;

    %if &xbeta2vars^=%str() %then %do;

      array xb2vecin(&num_repeatvars)  &xbeta2vars;
      array xb2vec(&num_repeatvars)    xbeta2_rep1-xbeta2_rep&num_repeatvars;

      do kxbeta2vars = 1 to &num_repeatvars;
        xb2vec(kxbeta2vars) = xb2vecin(kxbeta2vars);
      end;

    %end;

    %if &xbeta3vars^=%str() %then %do;

      array xb3vecin(&num_repeatvars)  &xbeta3vars;
      array xb3vec(&num_repeatvars)    xbeta3_rep1-xbeta3_rep&num_repeatvars;

      do kxbeta3vars = 1 to &num_repeatvars;
        xb3vec(kxbeta3vars) = xb3vecin(kxbeta3vars);
      end;

    %end;



    if (lastrec) then do;
      if (dc1miss24message = 1) & (dc2miss24message = 1) then do;
        put "Note:  If all 24-hour recalls are missing then denominator integration is not needed (i.e. dencalc=n).";
      end;
    end;

    output _subj1rec;
  run;






  %if &r24vars_dietcomp1^=%str() %then %do;

    proc freq data=_subj1rec;
      tables &r24vars_dietcomp1 / missing;

      %if (&model14=1 | &model14=3) %then %do;
        format &r24vars_dietcomp1 r24errzf.;
        title%eval(&titles+1) "Frequency of 24-Hour Recall Variables for Dietary Component 1.  Selected Model Does Not Allow Zeros for Dietary Component 1.";
      %end;
      %else %if (&model14=2 | &model14=4) %then %do;
        format &r24vars_dietcomp1 r24gtzf.;
        title%eval(&titles+1) "Frequency of 24-Hour Recall Variables for Dietary Component 1";
      %end;

    run;
    title%eval(&titles+1);

  %end;

  %if &r24vars_dietcomp2^=%str() & (&model14=3 | &model14=4) %then %do;

    proc freq data=_subj1rec;
      tables &r24vars_dietcomp2 / missing;

      %if (&model14=3) %then %do;
        format &r24vars_dietcomp2 r24errzf.;
        title%eval(&titles+1) "Frequency of 24-Hour Recall Variables for Dietary Component 2.  Selected Model Does Not Allow Zeros for Dietary Component 2.";
      %end;
      %else %do;
        format &r24vars_dietcomp2 r24gtzf.;
        title%eval(&titles+1) "Frequency of 24-Hour Recall Variables for Dietary Component 2";
      %end;

    run;
    title%eval(&titles+1);

  %end;




  /****************************************************************************/
  /****************************************************************************/




  proc iml;

    use _subj1rec  var{ &recid,

                        %if &r24vars_dietcomp1^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc1r&hrec,
                          %end;
                        %end;

                        %if &r24vars_dietcomp2^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc2r&hrec,
                          %end;
                        %end;

                        %if &xbeta1vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta1_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta2vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta2_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta3vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta3_rep&hrec,
                          %end;
                        %end;

                        %if (&randeffset = 1) %then %do;
                          &var_u3,
                        %end;
                        %else %if (&randeffset = 2) %then %do;
                          &var_u2,
                        %end;
                        %else %if (&randeffset = 3) %then %do;
                          &var_u2, &var_u3, &cov_u2u3,
                        %end;
                        %else %if (&randeffset = 4) %then %do;
                          &var_u1,
                        %end;
                        %else %if (&randeffset = 5) %then %do;
                          &var_u1, &var_u3, &cov_u1u3,
                        %end;
                        %else %if (&randeffset = 6) %then %do;
                          &var_u1, &var_u2, &cov_u1u2,
                        %end;
                        %else %if (&randeffset = 7) %then %do;
                          &var_u1, &var_u2, &var_u3, &cov_u1u2, &cov_u1u3, &cov_u2u3,
                        %end;

                        %if (&model14 = 1) | (&model14 = 2) %then %do;
                          &var_e2, &lambda2,
                        %end;
                        %else %if (&model14 = 3) | (&model14 = 4) %then %do;
                          &var_e2, &var_e3, &cov_e2e3, &lambda2, &lambda3,
                        %end;
                        
                        %if (&model14 = 4 & %upcase(&link)=PROBIT) %then %do;
                          &corr_e1e3, &corr_e2e3,
                        %end;

                        %if ((&pred_t_dietcomp13=1) | (&pred_t_dietcomp13=3)) %then %do;

                          %if ((&model14=2) | (&model14=4)) %then %do;
                            &pred_t_xbeta1_prob,
                          %end;

                          &pred_t_xbeta2_amt, &pred_t_min_amt2,

                        %end;

                        %if ((&pred_t_dietcomp13=2) | (&pred_t_dietcomp13=3)) %then %do;

                          &pred_t_xbeta3_amt, &pred_t_min_amt3,

                        %end;

                        &lamt
                        };



    start gbc_inverse(z_vector, bc_lambda_scalar, minamt_scalar);

      /*** Computes the inverse of the Box-Cox transformation and ensures result >= minimum amount ***/

      if bc_lambda_scalar = 0 then do;
        notneg_gbc_inv = exp(z_vector);
      end;
      else do;
        temp_notneg = (1 + bc_lambda_scalar # z_vector) <> 0;
        notneg_gbc_inv = temp_notneg ## (1 / bc_lambda_scalar);
      end;

      gbc_inv = notneg_gbc_inv <> minamt_scalar;

      return (gbc_inv);

    finish gbc_inverse;



    start backtransform(n, tran_lambda_j, tran_center_j, tran_scale_j, xbeta_u_j, sigmae_jj, minamt_j);

      /*** When lambda is zero, compute the exact backtransformation for the dietary component ***/
      /*** When lambda is not zero, compute the 9 point backtransformation for the dietary component ***/

      /*** Abscissas and weights ***/

      x9pt = {-2.1, -1.3, -0.8, -0.5, 0.0, 0.5, 0.8, 1.3, 2.1};
      w9pt_given = {0.063345, 0.080255, 0.070458, 0.159698, 0.252489, 0.159698, 0.070458, 0.080255, 0.063345};
      w9pt = w9pt_given / sum(w9pt_given);


      gstar = j(n, 1, 0);

      if tran_lambda_j = 0 then do;
        notneg_gstar = exp( tran_center_j + tran_scale_j # xbeta_u_j + ((tran_scale_j##2) # sigmae_jj / 2) );
        gstar = notneg_gstar <> minamt_j;
      end;
      else do;
        do qq = 1 to 9;

          bc_qq9pt = tran_center_j + tran_scale_j # (xbeta_u_j + x9pt[qq, 1] # (sigmae_jj)##0.5);
          g_inv_qq9pt = gbc_inverse(bc_qq9pt, tran_lambda_j, minamt_j);

          gstar = gstar + w9pt[qq, 1] # g_inv_qq9pt;

        end;
      end;

      return (gstar);

    finish backtransform;



    start mlogintn(u) global(covmatu, nrowcovu, covmate, nrowcove, inflate_constant,

                        %if &r24vars_dietcomp1^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc1r&hrec, dc1z&hrec,
                          %end;
                        %end;

                        %if &r24vars_dietcomp2^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc2r&hrec, dc2z&hrec,
                          %end;
                        %end;

                        %if &xbeta1vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta1_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta2vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta2_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta3vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta3_rep&hrec,
                          %end;
                          &corr_e1e3, &corr_e2e3,
                        %end;

                        %if ((&pred_t_dietcomp13=1) | (&pred_t_dietcomp13=3)) %then %do;

                          %if ((&model14=2) | (&model14=4)) %then %do;
                            &pred_t_xbeta1_prob,
                          %end;

                          &lambda2, &pred_t_xbeta2_amt, &pred_t_min_amt2,

                        %end;

                        %if ((&pred_t_dietcomp13=2) | (&pred_t_dietcomp13=3)) %then %do;

                          &lambda3, &pred_t_xbeta3_amt, &pred_t_min_amt3,

                        %end;

                        &lamt, pi);



      /****************************************************************************/
      /*** Calculate gt                                                         ***/
      /****************************************************************************/

      nsampsize_in_backtr   = 1;
      tran_center_in_backtr = 0;
      tran_scale_in_backtr  = 1;

      %if ((&pred_t_dietcomp13=1) | (&pred_t_dietcomp13=3)) %then %do;

        /*** Calculate probability for gt calculation ***/

        %if (&model14=1 | &model14=3) %then %do;

          prob_for_gt = 1;  /*** Selected model assumes dietary component 1 is consumed daily  ***/

        %end;
        %else %if (&model14=2 | &model14=4) %then %do;

          xbetaprobu1 = &pred_t_xbeta1_prob+&imlu1;

          %if (%upcase(&link)^=PROBIT) %then %do;
            if xbetaprobu1 > 0 then do;
              exp_neg1_xbetaprobu1    = exp(-xbetaprobu1);
              if exp_neg1_xbetaprobu1 = 0 then exp_neg1_xbetaprobu1=0.000000000001;
              prob_for_gt = 1/(1+exp_neg1_xbetaprobu1);
            end;
            else do;
              exp_pos1_xbetaprobu1    = exp(xbetaprobu1);
              if exp_pos1_xbetaprobu1 = 0 then exp_pos1_xbetaprobu1=0.000000000001;
              prob_for_gt = exp_pos1_xbetaprobu1/(1+exp_pos1_xbetaprobu1);
            end;
          %end;
          %else %do;
            prob_for_gt = cdf('normal', xbetaprobu1);
            if prob_for_gt > 1/1.000000000001
              then prob_for_gt = 1/1.000000000001;
            if prob_for_gt < 0.000000000001
              then prob_for_gt = 0.000000000001/1.000000000001;
          %end;

        %end;


        /**************************************************************************/

        covmate11 = covmate[1,1];

        pred_t_xbeta2_amt_u2  = &pred_t_xbeta2_amt + &imlu2;
        half_min_amt2         = 0.5 # &pred_t_min_amt2;

        g_star2 = backtransform(nsampsize_in_backtr, &lambda2, tran_center_in_backtr,
                                tran_scale_in_backtr, pred_t_xbeta2_amt_u2, covmate11, half_min_amt2);


        %if (&pred_t_dietcomp13=1) %then %do;

          if &lamt=0 then do;
            inflated_prob_amt = inflate_constant # (prob_for_gt # g_star2);
            gt_temp = log(inflated_prob_amt);
            gt = max(gt_temp, 0.000000000001);
          end;
          else if &lamt>0 then do;
            gt_temp = ((prob_for_gt # g_star2)##&lamt)/&lamt;
            gt = max(gt_temp, 0.000000000001);
          end;

        %end;

      %end;

      %if ((&pred_t_dietcomp13=2) | (&pred_t_dietcomp13=3)) %then %do;

        covmate22 = covmate[2,2];

        pred_t_xbeta3_amt_u3  = &pred_t_xbeta3_amt + &imlu3;
        half_min_amt3         = 0.5 # &pred_t_min_amt3;

        g_star3 = backtransform(nsampsize_in_backtr, &lambda3, tran_center_in_backtr,
                                tran_scale_in_backtr, pred_t_xbeta3_amt_u3, covmate22, half_min_amt3);


        %if (&pred_t_dietcomp13=2) %then %do;

          if &lamt=0 then do;
            inflated_prob_amt = inflate_constant # (g_star3);
            gt_temp = log(inflated_prob_amt);
            gt = max(gt_temp, 0.000000000001);
          end;
          else if &lamt>0 then do;
            gt_temp = ((g_star3)##&lamt)/&lamt;
            gt = max(gt_temp, 0.000000000001);
          end;

        %end;
        %else %if (&pred_t_dietcomp13=3) %then %do;

          if g_star3=0 then g_star3=0.000000000001;
          
          if &lamt=0 then do;
            inflated_prob_amt = inflate_constant # ((&pred_t_constant # prob_for_gt # g_star2 / g_star3));
            gt_temp = log(inflated_prob_amt);
            gt = max(gt_temp, 0.000000000001);
          end;
          else if &lamt>0 then do;
            gt_temp = ((&pred_t_constant # prob_for_gt # g_star2 / g_star3)##&lamt)/&lamt;
            gt = max(gt_temp, 0.000000000001);
          end;

        %end;

      %end;



      /****************************************************************************/
      /*** Calculate log_f_gr_givu                                              ***/
      /****************************************************************************/


      %if &model14=1 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp ***/

            logp = 0;  /*** For selected model, log(p) = log(1) = 0 ***/


            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            logp =           .;
            lognorme2 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=2 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp and log1_p ***/

            if dc1r&krec >= 0 then do;
              xbeta1u1    = xbeta1_rep&krec + &imlu1;

              %if (%upcase(&link)^=PROBIT) %then %do;
                if xbeta1u1 > 0 then do;
                  exp_neg1_xbeta1u1    = exp(-xbeta1u1);
                  if exp_neg1_xbeta1u1 = 0 then exp_neg1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = -log(1+exp_neg1_xbeta1u1);
                  else if dc1r&krec = 0 then log1_p = log( exp_neg1_xbeta1u1/(1+exp_neg1_xbeta1u1) );
                end;
                else do;
                  exp_pos1_xbeta1u1    = exp(xbeta1u1);
                  if exp_pos1_xbeta1u1 = 0 then exp_pos1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = log( exp_pos1_xbeta1u1/(1+exp_pos1_xbeta1u1) );
                  else if dc1r&krec = 0 then log1_p = -log(1+exp_pos1_xbeta1u1);
                end;
              %end;
              %else %do;
                temp = cdf('normal', xbeta1u1);
                if temp > 1/1.000000000001
                  then temp = 1/1.000000000001;
                if temp < 0.000000000001/1.000000000001
                  then temp = 0.000000000001/1.000000000001;
                logp=log(temp);
                log1_p=log(1-temp);
              %end;
            end;


            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc1r&krec = 0 then daycontrib  = log1_p;
            else daycontrib = 0;

            log_f_gr_givu = log_f_gr_givu + daycontrib;

            xbeta1u1 =       .;
            logp =           .;
            log1_p =         .;
            lognorme2 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=3 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) & (&r24vars_dietcomp2^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp ***/

            logp = 0;  /*** For selected model, log(p) = log(1) = 0 ***/


            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 & dc2r&krec > 0 then do;

              evect = ( dc1z&krec - xbeta2_rep&krec - &imlu2 ) || ( dc2z&krec - xbeta3_rep&krec - &imlu3 );
              logbivnorme2e3 = log( 1 / sqrt(((2#pi)##nrowcove) # det(covmate)) ) - 0.5#( evect * inv(covmate) * (evect`) );
              daycontrib  = logp + logbivnorme2e3;

            end;
            else if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = lognorme3;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            logp =           .;
            logbivnorme2e3 = .;
            lognorme2 =      .;
            lognorme3 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=4 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) & (&r24vars_dietcomp2^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp and log1_p ***/
            if dc1r&krec >= 0 then do;
              xbeta1u1    = xbeta1_rep&krec + &imlu1;

              %if (%upcase(&link)^=PROBIT) %then %do;
                if xbeta1u1 > 0 then do;
                  exp_neg1_xbeta1u1    = exp(-xbeta1u1);
                  if exp_neg1_xbeta1u1 = 0 then exp_neg1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = -log(1+exp_neg1_xbeta1u1);
                  else if dc1r&krec = 0 then log1_p = log( exp_neg1_xbeta1u1/(1+exp_neg1_xbeta1u1) );
                end;
                else do;
                  exp_pos1_xbeta1u1    = exp(xbeta1u1);
                  if exp_pos1_xbeta1u1 = 0 then exp_pos1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = log( exp_pos1_xbeta1u1/(1+exp_pos1_xbeta1u1) );
                  else if dc1r&krec = 0 then log1_p = -log(1+exp_pos1_xbeta1u1);
                end;
              %end;
              %else %do;
                if (dc1r&krec = 0 & dc2r&krec > 0) then do;
                  _z2 = (dc2z&krec - xbeta3_rep&krec - &imlu3)/sqrt(covmate[2,2]);
                  _m = &corr_e1e3 * _z2;
                  _s = sqrt(1 - &corr_e1e3**2);
                  _w = (-xbeta1u1 - _m) / _s;
                  l_b = probnorm(_w);
                  log1_p = log(l_b);
                  logp = log(1-l_b);
                end;
                else if (dc1r&krec > 0 & dc2r&krec > 0) then do;
                  _z1 = (dc1z&krec - xbeta2_rep&krec - &imlu2)/sqrt(covmate[1,1]);
                  _z2 = (dc2z&krec - xbeta3_rep&krec - &imlu3)/sqrt(covmate[2,2]);
                  _m = &corr_e1e3 / (1 - &corr_e2e3**2) * (_z2 - &corr_e2e3 * _z1);
                  _s = sqrt((1 - &corr_e1e3**2 - &corr_e2e3**2) / (1 - &corr_e2e3**2));
                  _w = (-xbeta1u1 - _m) / _s;
                  l_b = 1 - probnorm(_w);
                  logp = log(l_b);
                  log1_p = log(1-l_b);
                end;
                else do;
                  logp = .;
                  log1_p = .;
                end;
              %end;
            end;

            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 & dc2r&krec > 0 then do;

              evect = ( dc1z&krec - xbeta2_rep&krec - &imlu2 ) || ( dc2z&krec - xbeta3_rep&krec - &imlu3 );
              logbivnorme2e3 = log( 1 / sqrt(((2#pi)##nrowcove) # det(covmate)) ) - 0.5#( evect * inv(covmate) * (evect`) );
              daycontrib  = logp + logbivnorme2e3;

            end;
            else if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc1r&krec = 0 & dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = log1_p + lognorme3;

            end;
            else if dc1r&krec = 0 then do;

              daycontrib  = log1_p;

            end;
            else if dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = lognorme3;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            xbeta1u1 =       .;
            logp =           .;
            log1_p =         .;
            logbivnorme2e3 = .;
            lognorme2 =      .;
            lognorme3 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;



      /****************************************************************************/
      /*** Calculate the negative log of the numerator integrand                ***/
      /****************************************************************************/


      fn      =  - log( gt )
                 - log_f_gr_givu
                 - log( 1 / sqrt(((2#pi)##nrowcovu) # det(covmatu)) ) + 0.5#( u * inv(covmatu) * (u`) ) ;

      return(fn);
    finish mlogintn;



    start mlogintd(u) global(covmatu, nrowcovu, covmate, nrowcove,

                        %if &r24vars_dietcomp1^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc1r&hrec, dc1z&hrec,
                          %end;
                        %end;

                        %if &r24vars_dietcomp2^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            dc2r&hrec, dc2z&hrec,
                          %end;
                        %end;

                        %if &xbeta1vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta1_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta2vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta2_rep&hrec,
                          %end;
                        %end;

                        %if &xbeta3vars^=%str() %then %do;
                          %do hrec = 1 %to &num_repeatvars;
                            xbeta3_rep&hrec,
                          %end;
                          &corr_e1e3, &corr_e2e3,
                        %end;

                        pi);



      /****************************************************************************/
      /*** Calculate log_f_gr_givu                                              ***/
      /****************************************************************************/


      %if &model14=1 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp ***/

            logp = 0;  /*** For selected model, log(p) = log(1) = 0 ***/


            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            logp =           .;
            lognorme2 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=2 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp and log1_p ***/

            if dc1r&krec >= 0 then do;
              xbeta1u1    = xbeta1_rep&krec + &imlu1;

              %if (%upcase(&link)^=PROBIT) %then %do;
                if xbeta1u1 > 0 then do;
                  exp_neg1_xbeta1u1    = exp(-xbeta1u1);
                  if exp_neg1_xbeta1u1 = 0 then exp_neg1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = -log(1+exp_neg1_xbeta1u1);
                  else if dc1r&krec = 0 then log1_p = log( exp_neg1_xbeta1u1/(1+exp_neg1_xbeta1u1) );
                end;
                else do;
                  exp_pos1_xbeta1u1    = exp(xbeta1u1);
                  if exp_pos1_xbeta1u1 = 0 then exp_pos1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = log( exp_pos1_xbeta1u1/(1+exp_pos1_xbeta1u1) );
                  else if dc1r&krec = 0 then log1_p = -log(1+exp_pos1_xbeta1u1);
                end;
              %end;
              %else %do;
                temp = cdf('normal', xbeta1u1);
                if temp > 1/1.000000000001
                  then temp = 1/1.000000000001;
                if temp < 0.000000000001/1.000000000001
                  then temp = 0.000000000001/1.000000000001;
                logp=log(temp);
                log1_p=log(1-temp);
              %end;
            end;

            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc1r&krec = 0 then daycontrib  = log1_p;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            xbeta1u1 =       .;
            logp =           .;
            log1_p =         .;
            lognorme2 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=3 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) & (&r24vars_dietcomp2^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp ***/

            logp = 0;  /*** For selected model, log(p) = log(1) = 0 ***/


            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 & dc2r&krec > 0 then do;

              evect = ( dc1z&krec - xbeta2_rep&krec - &imlu2 ) || ( dc2z&krec - xbeta3_rep&krec - &imlu3 );
              logbivnorme2e3 = log( 1 / sqrt(((2#pi)##nrowcove) # det(covmate)) ) - 0.5#( evect * inv(covmate) * (evect`) );
              daycontrib  = logp + logbivnorme2e3;

            end;
            else if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = lognorme3;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            logp =           .;
            logbivnorme2e3 = .;
            lognorme2 =      .;
            lognorme3 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;
      %else %if &model14=4 %then %do;


        log_f_gr_givu = 0;
        %if (&r24vars_dietcomp1^=%str()) & (&r24vars_dietcomp2^=%str()) %then %do;
          %do krec = 1 %to &num_repeatvars;


            /*** Calculate logp and log1_p ***/

            if dc1r&krec >= 0 then do;
              xbeta1u1    = xbeta1_rep&krec + &imlu1;

              %if (%upcase(&link)^=PROBIT) %then %do;
                if xbeta1u1 > 0 then do;
                  exp_neg1_xbeta1u1    = exp(-xbeta1u1);
                  if exp_neg1_xbeta1u1 = 0 then exp_neg1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = -log(1+exp_neg1_xbeta1u1);
                  else if dc1r&krec = 0 then log1_p = log( exp_neg1_xbeta1u1/(1+exp_neg1_xbeta1u1) );
                end;
                else do;
                  exp_pos1_xbeta1u1    = exp(xbeta1u1);
                  if exp_pos1_xbeta1u1 = 0 then exp_pos1_xbeta1u1=0.000000000001;
                  if dc1r&krec > 0 then logp = log( exp_pos1_xbeta1u1/(1+exp_pos1_xbeta1u1) );
                  else if dc1r&krec = 0 then log1_p = -log(1+exp_pos1_xbeta1u1);
                end;
              %end;
              %else %do;
                if (dc1r&krec = 0 & dc2r&krec > 0) then do;
                  _z2 = (dc2z&krec - xbeta3_rep&krec - &imlu3)/sqrt(covmate[2,2]);
                  _m = &corr_e1e3 * _z2;
                  _s = sqrt(1 - &corr_e1e3**2);
                  _w = (-xbeta1u1 - _m) / _s;
                  l_b = probnorm(_w);
                  log1_p = log(l_b);
                  logp = log(1-l_b);
                end;
                else if (dc1r&krec > 0 & dc2r&krec > 0) then do;
                  _z1 = (dc1z&krec - xbeta2_rep&krec - &imlu2)/sqrt(covmate[1,1]);
                  _z2 = (dc2z&krec - xbeta3_rep&krec - &imlu3)/sqrt(covmate[2,2]);
                  _m = &corr_e1e3 / (1 - &corr_e2e3**2) * (_z2 - &corr_e2e3 * _z1);
                  _s = sqrt((1 - &corr_e1e3**2 - &corr_e2e3**2) / (1 - &corr_e2e3**2));
                  _w = (-xbeta1u1 - _m) / _s;
                  l_b = 1 - probnorm(_w);
                  logp = log(l_b);
                  log1_p = log(1-l_b);
                end;
                else do;
                  logp = .;
                  log1_p = .;
                end;
              %end;
            end;

            /*** Calculate contribution for this day ***/

            if dc1r&krec > 0 & dc2r&krec > 0 then do;

              evect = ( dc1z&krec - xbeta2_rep&krec - &imlu2 ) || ( dc2z&krec - xbeta3_rep&krec - &imlu3 );
              logbivnorme2e3 = log( 1 / sqrt(((2#pi)##nrowcove) # det(covmate)) ) - 0.5#( evect * inv(covmate) * (evect`) );
              daycontrib  = logp + logbivnorme2e3;

            end;
            else if dc1r&krec > 0 then do;

              lognorme2 = log( 1 / sqrt( 2 # pi # covmate[1,1])) - 0.5#( ( dc1z&krec - xbeta2_rep&krec - &imlu2 )##2 )/covmate[1,1];
              daycontrib  = logp + lognorme2;

            end;
            else if dc1r&krec = 0 & dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = log1_p + lognorme3;

            end;
            else if dc1r&krec = 0 then do;

              daycontrib  = log1_p;

            end;
            else if dc2r&krec > 0 then do;

              lognorme3 = log( 1 / sqrt( 2 # pi # covmate[2,2])) - 0.5#( ( dc2z&krec - xbeta3_rep&krec - &imlu3 )##2 )/covmate[2,2];
              daycontrib  = lognorme3;

            end;
            else daycontrib = 0;


            log_f_gr_givu = log_f_gr_givu + daycontrib;

            xbeta1u1 =       .;
            logp =           .;
            log1_p =         .;
            logbivnorme2e3 = .;
            lognorme2 =      .;
            lognorme3 =      .;
            daycontrib =     .;

          %end;
        %end;


      %end;



      /****************************************************************************/
      /*** Calculate the negative log of the denominator integrand              ***/
      /****************************************************************************/


      fd      =  - log_f_gr_givu
                 - log( 1 / sqrt(((2#pi)##nrowcovu) # det(covmatu)) ) + 0.5#( u * inv(covmatu) * (u`) ) ;

      return(fd);
    finish mlogintd;




    options nonotes;


    /*************************************************************************************/
    /*************************************************************************************/


    do data;

      read next;

      pi = arcos(-1);
      inflate_constant = 10##100;   *** constant to inflate the amount*probability value when lamt=0 ***;

      /*** 24-hour recall observed if recall >= 0. ***/

      dc1totcons = 0;
      dc1nirep = 0;

      %if &r24vars_dietcomp1^=%str() %then %do;
        %do krecdc1 = 1 %to &num_repeatvars;
          dc1z&krecdc1 = .;
          dc1consume&krecdc1 = .;
          if dc1r&krecdc1>=0 then do;
            dc1z&krecdc1 = ( dc1r&krecdc1##&lambda2 - 1 ) / &lambda2;
            dc1consume&krecdc1 = (dc1r&krecdc1>0);
            dc1totcons = dc1totcons + dc1consume&krecdc1;
            dc1nirep = dc1nirep + 1;
          end;
        %end;
      %end;


      dc2totcons = 0;
      dc2nirep = 0;

      %if &r24vars_dietcomp2^=%str() %then %do;
        %do krecdc2 = 1 %to &num_repeatvars;
          dc2z&krecdc2 = .;
          dc2consume&krecdc2 = .;
          if dc2r&krecdc2>=0 then do;
            dc2z&krecdc2 = ( dc2r&krecdc2##&lambda3 - 1 ) / &lambda3;
            dc2consume&krecdc2 = (dc2r&krecdc2>0);
            dc2totcons = dc2totcons + dc2consume&krecdc2;
            dc2nirep = dc2nirep + 1;
          end;
        %end;
      %end;



      /*** Covariance matrix for random effect(s) ***/

      &imlcovmatu;
      nrowcovu = nrow(covmatu);


      /*** Covariance matrix for e2 or (e2,e3) ***/

      &imlcovmate;
      nrowcove = nrow(covmate);


      /***************************************************************************************/
      /*** Optimization                                                                    ***/
      /***************************************************************************************/

      optn = {0 0};

      /*** Optimization for denominator integral ***/

      %if %upcase(&dencalc)=Y %then %do;

          ustart  = j(1, nrowcovu, 0.0);
          call nlpqn(rcd, ufinald, "mlogintd", ustart, optn);
          call nlpfdd(fufinald, gufinald, hufinald, "mlogintd", ufinald);
          ustartn  = ufinald;

      %end;
      %else %do;
        ustartn  = j(1, nrowcovu, 0.0);
      %end;

      /*** Optimization for numerator integral ***/

      call nlpqn(rcn, ufinaln, "mlogintn", ustartn, optn);
      call nlpfdd(fufinaln, gufinaln, hufinaln, "mlogintn", ufinaln);




      /****************************************************/
      /*** Evaluate numerator and denominator integrals ***/
      /****************************************************/

      sumresn = 0;
      %if %upcase(&dencalc)=Y %then %do;
        sumresd = 0;
      %end;

      /*** See Table 25.10 of Abramowitz and Stegun (1972) ***/


      zn6gh = { - 2.350604973674492
                - 1.335849074013697
                - 0.436077411927617
                  0.436077411927617
                  1.335849074013697
                  2.350604973674492 };

      ezn6sq = exp( zn6gh##2 );

      wn6gh = { 0.004530009905509
                0.1570673203229
                0.7246295952244
                0.7246295952244
                0.1570673203229
                0.004530009905509 };





      call eigen(eivalsn, eivecsn, hufinaln);

      eimhalfn = eivalsn ## (-0.5);
      scalmatn = sqrt(2) # ( eivecsn * diag(eimhalfn) * eivecsn` );



      %if %upcase(&dencalc)=Y %then %do;

        call eigen(eivalsd, eivecsd, hufinald);

        eimhalfd = eivalsd ## (-0.5);
        scalmatd = sqrt(2) # ( eivecsd * diag(eimhalfd) * eivecsd` );

      %end;







      %if &dimcovu = 1 %then %do;

        do j1 = 1 to 6;
          zvect = zn6gh[j1];
          wezsq = wn6gh[j1] # ezn6sq[j1];

      %end;
      %else %if &dimcovu = 2 %then %do;

        do j1 = 1 to 6;
          do j2 = 1 to 6;
            zvect = zn6gh[j1] // zn6gh[j2];
            wezsq = wn6gh[j1] # ezn6sq[j1] # wn6gh[j2] # ezn6sq[j2];

      %end;
      %else %if &dimcovu = 3 %then %do;

        do j1 = 1 to 6;
          do j2 = 1 to 6;
            do j3 = 1 to 6;
              zvect = zn6gh[j1] // zn6gh[j2] // zn6gh[j3];
              wezsq = wn6gh[j1] # ezn6sq[j1] # wn6gh[j2] # ezn6sq[j2] # wn6gh[j3] # ezn6sq[j3];

      %end;

              avectn = ufinaln + (scalmatn * zvect)`;
              sumresn = sumresn + exp(-mlogintn(avectn)) # wezsq;

      %if %upcase(&dencalc)=Y %then %do;

              avectd = ufinald + (scalmatd * zvect)`;
              sumresd = sumresd + exp(-mlogintd(avectd)) # wezsq;

      %end;

      %if &dimcovu = 1 %then %do;

        end; /*** End of j1 loop ***/

      %end;
      %else %if &dimcovu = 2 %then %do;

          end; /*** End of j2 loop ***/
        end; /*** End of j1 loop ***/

      %end;
      %else %if &dimcovu = 3 %then %do;

            end; /*** End of j3 loop ***/
          end; /*** End of j2 loop ***/
        end; /*** End of j1 loop ***/

      %end;




      intresn = (2##(nrowcovu/2)) # (det(hufinaln))##(-0.5) # sumresn;

      &imlufinaln;


      %if %upcase(&dencalc)=Y %then %do;

        intresd = (2##(nrowcovu/2)) # (det(hufinald))##(-0.5) # sumresd;

        &imlufinald;


        /******************************************************************/
        /*** E[(T##lamt-1)/lamt|Cond]=E[T##lamt/lamt|Cond]-1/lamt       ***/
        /******************************************************************/
        %if %upcase(&boxcox_t_lamt)=N %then %do;
          indusint = intresn/intresd;  /*** No transformation.  Macro sets lamt to 1 because T##lamt/lamt is used. ***/
        %end;
        %else %if %upcase(&boxcox_t_lamt)=Y %then %do;
          if &lamt=0 then do;
            indusint = intresn/intresd - log(inflate_constant);
          end;
          else if &lamt>0 then do;
            indusint = intresn/intresd - 1/&lamt;
          end;
        %end;

      %end;
      %else %do;

        /******************************************************************/
        /*** E[(T##lamt-1)/lamt|Cond]=E[T##lamt/lamt|Cond]-1/lamt       ***/
        /******************************************************************/
        %if %upcase(&boxcox_t_lamt)=N %then %do;
          indusint = intresn;  /*** No transformation.  Macro sets lamt to 1 because T##lamt/lamt is used. ***/
        %end;
        %else %if %upcase(&boxcox_t_lamt)=Y %then %do;
          if &lamt=0 then do;
            indusint = intresn - log(inflate_constant);
          end;
          else if &lamt>0 then do;
            indusint = intresn - 1/&lamt;
          end;
        %end;

      %end;



      resvals   =  &recid || indusint || intresn
                   %if %upcase(&dencalc)=Y %then %str(|| intresd);
                   || dc1nirep || dc1totcons
                   %if ((&model14=3) | (&model14=4)) %then %str(|| dc2nirep || dc2totcons);
                   || rcn  || u1finaln  || u2finaln  || u3finaln
                   %if %upcase(&dencalc)=Y %then %str(|| rcd  || u1finald  || u2finald || u3finald);
                   ;

      resmat    =  resmat // resvals ;

    end;

    /*************************************************************************************/
    /*************************************************************************************/

    %if %upcase(&notesprt)=Y %then %do;
      options notes;
    %end;




    /*** Print resmat ***/

    varnames = {"&recid" "indusint" "intresn"
                %if %upcase(&dencalc)=Y %then %str("intresd");
                "dc1nirep" "dc1totcons"
                %if ((&model14=3) | (&model14=4)) %then %str("dc2nirep" "dc2totcons");
                "rcn" "u1finaln" "u2finaln" "u3finaln"
                %if %upcase(&dencalc)=Y %then %str("rcd" "u1finald" "u2finald" "u3finald");
                };

    create _resdata from resmat [colname=varnames];
    append from resmat;

  quit;



  /****************************************************************************/
  /****************************************************************************/

  proc freq data=_resdata;
    tables dc1nirep*dc1totcons
    %if ((&model14=3) | (&model14=4)) %then %str(dc2nirep*dc2totcons);
    / norow nocol nopercent missing;
    title%eval(&titles+1) "Number of Observed 24-Hour Recalls and Number of 24-Hour Recalls Greater Than 0";
  run;
  title%eval(&titles+1);


  proc freq data=_resdata;
    tables rcn
           %if %upcase(&dencalc)=Y %then %str(rcd);
           ;
    title%eval(&titles+1) "Return Code from Optimization Step of Adaptive Gaussian Quadrature.  Positive Value Indicates Successful Termination.";
  run;
  title%eval(&titles+1);



%mend predict_intake_density;  /*** End macro ***/


