/**************************************************************************** * * * SAS macro NLMixed_Bivariate fits a bivariate nonlinear mixed 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. * * * * The model for two foods/nutrients that are consumed every day is * * described in: * * Freedman LS, et al. The population distribution of ratios of usual * * intake of dietary components that are consume every day can be * * estimated from repeated 24-hour recalls. J Nutr. 2010;140:111-116. * * * * 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 = LOGIT. * * 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 = N. * * 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. * * The verbose option is useful for debugging. * * 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_Lambda = Box-Cox transformation parameter for the * * first food/nutrient. * * A2_Lambda = Box-Cox transformation parameter 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 (intital estimates for these * * parameters will be set to zero): * * z_e2e3 * * Cov_e2e3 * * Corr_e2e3 * * Cov_u1u2 Cov_u1u3 Cov_u2u3 * * Corr_u1u2 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 = , init_parms = , covars_prob1 = , covars_amt1 = , covars_amt2 = , link = LOGIT, lambda1 = , lambda2 = , var_u1 = , var_u2 = , var_u3 = , indep_u1 = N, indep_u2 = N, indep_u3 = N, threshold = N, replicate_var = , nloptions = , print = Y, ntitle = 2); %let print = %upcase(%substr(&print,1,1)); %let modeltype = %upcase(&modeltype); %let link = %upcase(&link); %let threshold = %upcase(%substr(&threshold,1,1)); %let indep_u1 = %upcase(%substr(&indep_u1,1,1)); %let indep_u2 = %upcase(%substr(&indep_u2,1,1)); %let indep_u3 = %upcase(%substr(&indep_u3,1,1)); %if (&link ^= PROBIT) %then %let link = LOGIT; %if (&modeltype ^= ONEPART) %then %let modeltype = TWOPART; %if (&modeltype = ONEPART) %then %let link = NONE; %if (&modeltype = ONEPART) %then %let var_u1 = 0; %if (&modeltype = ONEPART) %then %let threshold = N; %if (&threshold = Y) %then %let link = PROBIT; %if (&var_u1 ^= %str()) %then %let indep_u1 = Y; %if (&var_u2 ^= %str()) %then %let indep_u2 = Y; %if (&var_u3 ^= %str()) %then %let indep_u3 = Y; %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. ***/ proc iml; 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; close &init_parms; stddev_u = exp(logsd_u)`; cov_u = diag(stddev_u##2); dim_u = nrow(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; z_e2e3 = 0; z_gamma = 0; z_theta = 0; initcov = %if (&modeltype = ONEPART | &threshold ^= Y) %then z_e2e3 ||; %else z_gamma || z_theta ||; log_d` %if (&dim_u > 1) %then || s_lower`; %str(;) names = {%if (&modeltype = ONEPART | &threshold ^= Y) %then "Z_E2E3"; %else "Z_GAMMA" "Z_THETA"; %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;}; 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. */ drop LogSDu2 LogSDu3 %if (&modeltype ^= ONEPART) %then LogSDu1; %do i = 1 %to &dim_u; %let ii = %scan(&index_u,&i,%str( )); %if (&&var_u&ii ^= %str()) %then LOG_D&i; %if (&&var_u&ii ^= %str() | &&indep_u&ii = Y) %then %do j = 1 %to &dim_u; %if (&j < &i) %then S&i&j; %if (&j > &i) %then S&j&i; %end; %end; ; run; %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(;); /* 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. */ %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() | &&indep_u&ii = Y) %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; /* 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 * ξ %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 * ξ %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 * ξ %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. */ 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; /* estimate additional parameters. */ %do i = 1 %to &dim_u; %let ii = %scan(&index_u,&i,%str( )); estimate "VAR_U&ii" VAR_Uⅈ %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; /*** end of macro NLMixed_Bivariate ****************************************/