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

%macro NLMixed_Bivariate (data           = ,
                          subject        = ,
                          repeat         = ,
                          response1      = ,
                          response2      = ,
                          modeltype      = TWOPART,
                          init_parms     = ,
                          covars_prob1   = ,
                          covars_amt1    = ,
                          covars_amt2    = ,
                          link           = PROBIT,
                          lambda1        = ,
                          lambda2        = ,
                          var_u1         = ,
                          var_u2         = ,
                          var_u3         = ,
                          corr_u1u2      = ,
                          corr_u1u3      = ,
                          corr_u2u3      = ,
                          threshold      = Y,
                          replicate_var  = ,
                          init_cov_zero  = Y,
                          cov_u_sds      = Y,
                          nloptions      = ,
                          print          = Y,
                          ntitle         = 2);

%let print      = %upcase(%substr(&print,1,1));
%let modeltype  = %upcase(&modeltype);
%let link       = %upcase(&link);
%let init_cov_zero = %upcase(%substr(&init_cov_zero,1,1));
%let threshold  = %upcase(%substr(&threshold,1,1));

%if (&link      ^= LOGIT)   %then %let link = PROBIT;
%if (&link       = LOGIT)   %then %let threshold = N;
%if (&modeltype ^= ONEPART) %then %let modeltype = TWOPART;

%if (&modeltype = ONEPART) %then %do;
  %let link = NONE;
  %let var_u1 = 0;
  %let threshold = N;
  %end;

%if (&print = V) %then %let nloptions = &nloptions itdetails;

%if (&replicate_var ^= %str()) %then %do;
  footnote "Note: Standard Errors are not valid if data are not from a simple random sample";
  %end;

proc sort data=&data;  by &subject &repeat;  run;

 /*** determine number and order of random effects u. ***/

%let dim_u   = 0;
%let index_u = ;

%if (&modeltype ^= ONEPART & &var_u1 ^= 0) %then %do;
  %let dim_u   = %eval(&dim_u + 1);
  %let index_u = &index_u 1;
  %end;

%if (&var_u2 ^= 0) %then %do;
  %let dim_u   = %eval(&dim_u + 1);
  %let index_u = &index_u 2;
  %end;

%let dim_u_f1 = &dim_u;

%if (&var_u3 ^= 0) %then %do;
  %let dim_u   = %eval(&dim_u + 1);
  %let index_u = &index_u 3;
  %end;

 /*** determine number of covariates in each part of the model. ***/

data _null_;
  set &data;

  %if (&covars_prob1 = %str()) %then nx_p1 = 0%str(;);
  %else %do;
    array _p1 (*) &covars_prob1;
    nx_p1 = dim(_p1);
    %end;

  %if (&covars_amt1 = %str()) %then nx_a1 = 0%str(;);
  %else %do;
    array _a1 (*) &covars_amt1;
    nx_a1 = dim(_a1);
    %end;

  %if (&covars_amt2 = %str()) %then nx_a2 = 0%str(;);
  %else %do;
    array _a2 (*) &covars_amt2;
    nx_a2 = dim(_a2);
    %end;

  call symput("nx_p1",trim(left(put(nx_p1, 3.))));
  call symput("nx_a1",trim(left(put(nx_a1, 3.))));
  call symput("nx_a2",trim(left(put(nx_a2, 3.))));
  run;

 /*** calculate initial estimates of covariance parameters ***/
 /*** for random effects u and within-person errors e.     ***/

%if (&init_cov_zero = Y) %then %do;

proc iml;

  z_e2e3  = 0;
  z_gamma = 0;
  z_theta = 0;

  initcov = %if (&modeltype = ONEPART | &threshold ^= Y) %then z_e2e3;
            %else  z_gamma || z_theta; %str(;)

  names   = {%if (&modeltype = ONEPART | &threshold ^= Y) %then "Z_E2E3";
             %else "Z_GAMMA" "Z_THETA"; };

  %if (&cov_u_sds = Y & &dim_u > 0) %then %do;
    use &init_parms;
      read all var{%do i = 1 %to &dim_u;
                   %let ii = %scan(&index_u,&i,%str( ));
                   LogSDu&ii
                   %end;
                } into logsd_u;
    %if (&modeltype ^= ONEPART) %then %do;
      read all var{corr_u1u2} into corr_u1u2;
      %end;
    close &init_parms;

    stddev_u = exp(logsd_u)`;
    dim_u    = nrow(stddev_u);

    corr_u = i(dim_u);  
    %if (&modeltype ^= ONEPART) %then %do;
      corr_u[1,2] = corr_u1u2;
      if (corr_u[1,2] < -0.8) then corr_u[1,2] = -0.8;
      if (corr_u[1,2] >  0.8) then corr_u[1,2] =  0.8;
      corr_u[2,1] = corr_u[1,2];
      %end;

    cov_u  = corr_u # (stddev_u * stddev_u`);

 /* parameterize cov(u) as cov(u) = SDS', where S is a lower triangular matrix with 1's  */
 /* on the diagonal, and D is a diagonal matrix with positive values on the diagonal.    */

    root_u = root(cov_u)`;

    sqrt_d  = vecdiag(root_u);
    log_d   = log(sqrt_d##2);
    %if (&dim_u > 1) %then %do;
      s = root_u / repeat(sqrt_d`,dim_u,1);
      s_lower = symsqr(s[2:dim_u,1:(dim_u-1)]);
      %end;

  initcov = initcov || log_d`
            %if (&dim_u > 1) %then || s_lower`; %str(;)

  names   = names ||
            {%do i = 1 %to &dim_u;
               "LOG_D&i"
               %end;
             %do i = 2 %to &dim_u;
               %do j = 1 %to %eval(&i-1);
                 "S&i&j"
                 %end;
               %end; };

  %end;  /* %if (&cov_u_sds = Y & &dim_u > 0) %then %do */

  create _initcov from initcov [colname = names];
  append from initcov;
  close _initcov;
  quit;

 /* combine initial parameter estimates. */

  data _init_parms;
    merge &init_parms _initcov;

 /* drop uneeded parameters. */

    %if (&cov_u_sds = Y & &dim_u > 0) %then %do;
      drop LogSDu2 LogSDu3
           %if (&modeltype ^= ONEPART) %then LogSDu1 corr_u1u2;
           ;
      %end;
    run;

  %end;   /* %if (&init_cov_u = Y | &init_cov_e = Y) %then %do */

%else %do;
  data _init_parms;
    set &init_parms;
    run;
  %end;

%if (&print = V) %then %do;
  proc print data=_init_parms noobs;
    title%eval(&ntitle+1) "Initial Parameter Estimates for Bivariate Foods/Nutrients Model";
    run;
  %end;

 /*** fit model using proc NLMixed. ***/

      %if (&link = NONE)  %then title%eval(&ntitle+1) "Bivariate Foods/Nutrients Model"%str(;);
%else %if (&link = LOGIT) %then title%eval(&ntitle+1) "Bivariate Foods/Nutrients Model (Link = Logit)"%str(;);
%else                           title%eval(&ntitle+1) "Bivariate Foods/Nutrients Model (Link = Probit)"%str(;);

data parms_b;
  run;

data pred_x_b;
  run;

data conv_b;
  status = 1;
  run;

%if (&print = N) %then ods exclude all %str(;);

proc nlmixed  data=&data &nloptions;
  parms / data=_init_parms;

  %if (&replicate_var ^= %str()) %then replicate &replicate_var%str(;);    /* weight variable */

  %if (&lambda1 ^= %str()) %then A1_Lambda = &lambda1%str(;);
  %if (&lambda2 ^= %str()) %then A2_Lambda = &lambda2%str(;);

  %if (&var_u1 = 0) %then u1 = 0 %str(;);
  %if (&var_u2 = 0) %then u2 = 0 %str(;);
  %if (&var_u3 = 0) %then u3 = 0 %str(;);

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

  %if (%quote(&corr_u1u2) ^= %str()) %then corr_u1u2 = &corr_u1u2 %str(;);
  %if (%quote(&corr_u1u3) ^= %str()) %then corr_u1u3 = &corr_u1u3 %str(;);
  %if (%quote(&corr_u2u3) ^= %str()) %then corr_u2u3 = &corr_u2u3 %str(;);

 /* reparameterize the variance/covariance of e. */

  VAR_E2 = exp(2*A1_LogSDe);
  VAR_E3 = exp(2*A2_LogSDe);

  %if (&modeltype = ONEPART | &threshold ^= Y) %then %do;
    CORR_E1E3 = 0;
    CORR_E2E3 = (exp(2*Z_E2E3) - 1) / (exp(2*Z_E2E3) + 1);
    %end;
  %else %do;
    GAMMA = (exp(2*Z_GAMMA) - 1) / (exp(2*Z_GAMMA) + 1);
    THETA = (exp(2*Z_THETA) - 1) / (exp(2*Z_THETA) + 1);
    CORR_E1E3 = GAMMA * THETA;
    CORR_E2E3 = GAMMA * sqrt(1 - THETA**2);
    %end;

  COV_E1E3  = CORR_E1E3 * sqrt(VAR_E3);
  COV_E2E3  = CORR_E2E3 * sqrt(VAR_E2 * VAR_E3);

 /* reparameterize the variance/covariance of u. */

  %if (&cov_u_sds = Y & &dim_u > 0) %then %do;
    %do i = 1 %to &dim_u;
      %let ii = %scan(&index_u,&i,%str( ));
      %if (&&var_u&ii ^= %str()) %then D&i = &&var_u&ii%str(;);
      %else D&i = exp(LOG_D&i)%str(;);
      %if (&&var_u&ii ^= %str()) %then %do j = 1 %to &dim_u;
        %if (&j < &i) %then S&i&j = 0%str(;);
        %if (&j > &i) %then S&j&i = 0%str(;);
        %end;
     %end;

    %do i = 1 %to &dim_u;
      %let ii = %scan(&index_u,&i,%str( ));
      VAR_U&ii = D&i
           %do k = 1 %to %eval(&i-1);
             + D&k * S&i&k**2
             %end;
           ;
        %do j = 1 %to %eval(&i-1);
          %let jj = %scan(&index_u,&j,%str( ));
          COV_U&jj.U&ii = D&j * S&i&j
               %do k = 1 %to %eval(&j-1);
                 + D&k * S&i&k * S&j&k
                 %end;
               ;
          CORR_U&jj.U&ii = COV_U&jj.U&ii / sqrt(VAR_U&jj * VAR_U&ii);
          %end;
      %end;
    %end;  /* %if (&cov_u_sds = Y & &dim_u > 0) %then %do */

  %else %do;
    %do i = 1 %to &dim_u;
      %let ii = %scan(&index_u,&i,%str( ));
      %if (&&var_u&ii ^= %str()) %then VAR_U&ii = &&var_u&ii %str(;);
      %do j = 1 %to %eval(&i-1);
        %let jj = %scan(&index_u,&j,%str( ));
        %if (%quote(&&corr_u&jj.u&ii) ^= %str()) %then %do;
          CORR_U&jj.U&ii = %quote(&&corr_u&jj.u&ii);
    COV_U&jj.U&ii = sqrt(var_u&jj * var_u&ii) * corr_u&jj.u&ii;
          %end;
        %else %do;
          CORR_U&jj.U&ii = COV_U&jj.U&ii / sqrt(VAR_U&jj * VAR_U&ii);
          %end;
        %end;
      %end;
    %end;  /* %else %do */


 /* get response data. */

  %if (&modeltype = ONEPART) %then %do;
    if (&response1 > 0) then do;
      if (A1_Lambda = 0) then _boxcoxy1 = log(&response1);
      else _boxcoxy1 = (&response1**A1_Lambda - 1) / A1_Lambda;
      end;
    else do;
      _boxcoxy1 = .;
      end;
    %end;

  %else %do;
    if (&response1 > 0) then do;
      _y1 = 1;
      if (A1_Lambda = 0) then _boxcoxy1 = log(&response1);
      else _boxcoxy1 = (&response1**A1_Lambda - 1) / A1_Lambda;
      end;
    else if (&response1 = 0) then do;
      _y1 = 0;
      _boxcoxy1 = .;
      end;
    else do;
      _y1 = .;
      _boxcoxy1 = .;
      end;
    %end;

  if (&response2 > 0) then do;
    if (A2_Lambda = 0) then _boxcoxy2 = log(&response2);
    else _boxcoxy2 = (&response2**A2_Lambda - 1) / A2_Lambda;
    end;
  else do;
    _boxcoxy2 = .;
    end;

 /* calculate linear predictor for probability of consumption. */

  %if (&modeltype ^= ONEPART) %then %do;
    x_p1 = p1_intercept;
    %do i = 1 %to &nx_p1;
      %let xi = %scan(&covars_prob1,&i,%str( ));
      x_p1 = x_p1 + p1_&xi * &xi;
      %end;
    xu_p1 = x_p1 + u1;

 /* logit or probit link. */

    %if (&link = LOGIT) %then %do;
      _p1 = 1 /(1 + exp(-xu_p1));
      _c1 = probit(1 - _p1);
      %end;
    %else %do;
      _p1 = probnorm(xu_p1);
      _c1 = -xu_p1;
      %end;
    %end;

 /* calculate linear predictors for amount consumed. */

  x_a1 = a1_intercept;
  %do i = 1 %to &nx_a1;
    %let xi = %scan(&covars_amt1,&i,%str( ));
    x_a1 = x_a1 + a1_&xi * &xi;
    %end;
  xu_a1 = x_a1 + u2;

  x_a2 = a2_intercept;
  %do i = 1 %to &nx_a2;
    %let xi = %scan(&covars_amt2,&i,%str( ));
    x_a2 = x_a2 + a2_&xi * &xi;
    %end;
  xu_a2 = x_a2 + u3;

  if (_boxcoxy1 ^= .) then _z1 = (_boxcoxy1 - xu_a1) / sqrt(var_e2);
  else _z1 = .;

  if (_boxcoxy2 ^= .) then _z2 = (_boxcoxy2 - xu_a2) / sqrt(var_e3);
  else _z2 = .;

 /* calculate likelihood for probability of consumption. */

  %if (&modeltype = ONEPART) %then %do;
    l_b = 1;
    %end;

  %else %do;
    %if (&threshold ^= Y) %then %do;
      if (_y1 = 0) then l_b = 1 - _p1;
      else if (_y1 = 1) then l_b = _p1;
      else if (_boxcoxy2 ^= .) then l_b = 1;
      else l_b = .;
      %end;

 /* latent variable threshold model. */

    %else %do;
      if (_y1 = 0 & _boxcoxy2 ^= .) then do;
        _m = corr_e1e3 * _z2;
        _s = sqrt(1 - corr_e1e3**2);
        _w = (_c1 - _m) / _s;
        l_b = probnorm(_w);
        end;

      else if (_y1 = 1 & _boxcoxy2 ^= .) then do;
        _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 = (_c1 - _m) / _s;
        l_b = 1 - probnorm(_w);
        end;

      else if (_y1 = 0) then l_b = 1 - _p1;
      else if (_y1 = 1) then l_b = _p1;
      else if (_boxcoxy2 ^= .) then l_b = 1;
      else l_b = .;
      %end;
    %end;

  if (l_b ^= .) then ll_b = log(l_b);
  else ll_b = .;

 /* calculate likelihood for amount consumed. */

  pi = arcos(-1);

  if (_z1 ^= . & _z2 ^= .) then
    ll_n = -log(2 * pi * sqrt(VAR_E2 * VAR_E3 * (1-CORR_E2E3**2))) -
           (_z1**2 - 2 * CORR_E2E3 * _z1 * _z2 +_z2**2) / (2 * (1-CORR_E2E3**2)) +
            (A1_Lambda - 1)*log(&response1) + (A2_Lambda - 1)*log(&response2);

  else if (_z1 ^= .) then 
    ll_n = -log(sqrt(2 * pi * VAR_E2)) - _z1**2 / 2 +
            (A1_Lambda - 1)*log(&response1);

  else if (_z2 ^= .) then 
    ll_n = -log(sqrt(2 * pi * VAR_E3)) - _z2**2 / 2 +
            (A2_Lambda - 1)*log(&response2);

  %if (&modeltype ^= ONEPART) %then %do;
    else if (_y1 = 0) then ll_n = 0;
    %end;

  else ll_n = .;

  ll = ll_b + ll_n;

  if (&response1 > .Z) then _response = &response1;
  else _response = &response2;

  model _response ~ general(ll);

 /* specify random effects. */

%if (&dim_u > 0) %then %do;
  random %do i = 1 %to &dim_u;
           %let ii = %scan(&index_u,&i,%str( ));
           u&ii
           %end;
         ~ normal([%do i = 1 %to &dim_u;
                     0
                     %if (&i < &dim_u) %then ,;
                     %end;
                  ],
                  [%do i = 1 %to &dim_u;
                    %let ii = %scan(&index_u,&i,%str( ));
                    %do j = 1 %to %eval(&i-1);
                      %let jj = %scan(&index_u,&j,%str( ));
                       cov_u&jj.u&ii,
                        %end;
                     var_u&ii
                     %if (&i < &dim_u) %then ,;
                     %end;
                  ])
        subject=&subject;
  %end;

 /* estimate additional parameters. */

  %do i = 1 %to &dim_u;
    %let ii = %scan(&index_u,&i,%str( ));
     estimate "VAR_U&ii" VAR_U&ii;
     %end;

  %do i = 1 %to &dim_u;
    %let ii = %scan(&index_u,&i,%str( ));
    %do j = %eval(&i+1) %to &dim_u;
      %let jj = %scan(&index_u,&j,%str( ));
       estimate "COV_U&ii.U&jj" COV_U&ii.U&jj;
       %end;
     %end;

  %do i = 1 %to &dim_u;
    %let ii = %scan(&index_u,&i,%str( ));
    %do j = %eval(&i+1) %to &dim_u;
      %let jj = %scan(&index_u,&j,%str( ));
       estimate "CORR_U&ii.U&jj" CORR_U&ii.U&jj;
       %end;
     %end;

  estimate "VAR_E2" VAR_E2;
  estimate "VAR_E3" VAR_E3;

  *%if (&modeltype ^= ONEPART & &threshold = Y) %then estimate "GAMMA" GAMMA %str(;);
  *%if (&modeltype ^= ONEPART & &threshold = Y) %then estimate "THETA" THETA %str(;);

  %if (&modeltype ^= ONEPART & &threshold = Y) %then estimate "COV_E1E3" COV_E1E3 %str(;);
  estimate "COV_E2E3" COV_E2E3;

  %if (&modeltype ^= ONEPART & &threshold = Y) %then estimate "CORR_E1E3" CORR_E1E3 %str(;);
  estimate "CORR_E2E3" CORR_E2E3;

  %if (&modeltype ^= ONEPART) %then predict x_p1 out=pred_x_p1%str(;);
  predict x_a1 out=pred_x_a1;
  predict x_a2 out=pred_x_a2;

  ods output ParameterEstimates=_parms_b;
  ods output AdditionalEstimates=_adparms_b;
  ods output FitStatistics=fit_b;
  ods output ConvergenceStatus=conv_b;
  run;

title%eval(&ntitle+1);

%if (&print = N) %then ods exclude none %str(;);

data _null_;
  set conv_b;
  call symput("status",trim(left(put(status, 6.))));
  run;

%if (&status ^= 0) %then %goto exit;

 /*** save parameter estimates and predicted values. ***/

proc transpose data=_parms_b out=_parms_b(drop=_name_);
  id parameter;
  var estimate;
  run;

proc transpose data=_adparms_b out=_adparms_b(drop=_name_);
   id label;
   var estimate;
   run;

data parms_b;
  merge _parms_b _adparms_b;

 /* add fixed parameters. */

  %if (&lambda1 ^= %str()) %then A1_Lambda = &lambda1%str(;);
  %if (&lambda2 ^= %str()) %then A2_Lambda = &lambda2%str(;);

  %if (&modeltype ^= ONEPART & &var_u1 = 0) %then VAR_U1 = 0%str(;);
  %if (&var_u2 = 0) %then VAR_U2 = 0%str(;);
  %if (&var_u3 = 0) %then VAR_U3 = 0%str(;);

  %if (&modeltype ^= ONEPART) %then %do;
    %if (&var_u1 = 0 | &var_u2 = 0) %then COV_U1U2 = 0%str(;);
    %if (&var_u1 = 0 | &var_u3 = 0) %then COV_U1U3 = 0%str(;);
    %end;
  %if (&var_u2 = 0 | &var_u3 = 0) %then COV_U2U3 = 0%str(;);

  %if (&modeltype ^= ONEPART) %then %do;
    %if (&var_u1 = 0 | &var_u2 = 0) %then CORR_U1U2 = 0%str(;);
    %if (&var_u1 = 0 | &var_u3 = 0) %then CORR_U1U3 = 0%str(;);
    %end;
  %if (&var_u2 = 0 | &var_u3 = 0) %then CORR_U2U3 = 0%str(;);
  run;

data pred_x_b;
  merge &data
        %if (&modeltype ^= ONEPART) %then pred_x_p1(keep=&subject &repeat pred rename=(pred=pred_x_p1));
                                           pred_x_a1(keep=&subject &repeat pred rename=(pred=pred_x_a1))
                                           pred_x_a2(keep=&subject &repeat pred rename=(pred=pred_x_a2));
    by &subject &repeat;

  label %if (&modeltype ^= ONEPART) %then pred_x_p1 = " ";
        pred_x_a1 = " " pred_x_a2 = " ";
  run;

 /*** delete unneeded data sets. */

proc datasets lib=work nolist;
  delete _initcov _init_parms _parms_b _adparms_b
        %if (&modeltype ^= ONEPART) %then pred_x_p1;
        pred_x_a1 pred_x_a2;
  run;
  quit;

%exit:

title%eval(&ntitle+1);

%if (&replicate_var ^= %str()) %then footnote %str(;);

%mend NLMixed_Bivariate;
