data demoadv;
  set nh.demoadv;
if w0304_0 ne . ;   
run ;

First, select only those people with dietary data by selecting those without missing BRR weights.

data adultf;
set demoadv;
if riagendr= 2 and ridageyr>= 19 ;
run ;

Because the INDIVINT macro also requires a dataset with 1 record per person, a dataset called adultf for adult females ages 19 years and older is created.

data day1;
set demoadv;
if riagendr= 2 and ridageyr>= 19 ;
d_milk=d_milk_d1;
day= 1 ;
run ;

data day2;
set demoadv;
if riagendr= 2 and ridageyr>= 19 ;
d_milk=d_milk_d2;
day= 2 ;
run ;

The variables d_milk_d1 and d_milk_d2 are calculated variables representing total milk consumed on days 1 and 2 respectively from all foods and beverages (see Module 4, Task 2 and Module 9, Task 4 for more details).  To create a dataset with 2 records per person, the demoadv dataset is set 2 times to create 2 datasets, one where day=1 and one where day=2.  The same variable name, d_milk, is used for calcium on both days.  It is created by setting it equal to d_milk_d1 for day 1 and d_milk_d2  for day 2. 

Finally, these data sets are appended.

%include   'C:\NHANES\Macros\mixtran_macro_v1.1.sas' ;
%include   'C:\NHANES\Macros\indivint_macro_v1.1.sas' ;
%include   'C:\NHANES\Macros\distrib_macro_v1.1.sas' ;

This code reads the MIXTRAN, DISTRIB, and INDIVINT macros into SAS so that these macros may be called.

title2 "Task 3 - Milk and Osteoporosis Treatment" ;

This title should be printed on all output.

%macro BRR212(data, response, datam, var12, regresponse, foodtype, subject, repeat, covars_prob, covars_amt, outlib, modeltype, lambda, covars_reg, seq, weekend, vargroup, numvargroups, subgroup, start_val1, start_val2, start_val3, vcontrol, nloptions, titles, printlevel, cutpts, ncutpts, nsim_mc, byvar, t0, t1, final);

The start of the BRR212 macro is defined.  All of the terms inside the parentheses are the macro variables.

%MIXTRAN (data=&data, response=&response, foodtype=&foodtype, subject=&subject, repeat=&repeat, covars_prob=&covars_prob, covars_amt=&covars_amt, outlib=&outlib, modeltype=&modeltype,                 lambda=&lambda, replicate_var=w0304_0, seq=&seq, weekend=&weekend, vargroup=&vargroup, numvargroups= &numvargroups, subgroup=&subgroup, start_val1=&start_val1, start_val2=&start_val2, start_val3= &start_val3, vcontrol= &vcontrol, nloptions=&nloptions, titles= &titles, printlevel= &printlevel)

Within the BRR212 macro, the MIXTRAN macro is called.  All of the variables preceded by & will be defined by the BRR212 macro call.  The only variable without an & is the replicate_var macro variable; it is set to w0304_0 for the first run.

data parampred(keep = &subject a_var_u2 a_var_e a_lambda x2b2);
if (_n_ = 1 ) then set & outlib.._ param_unc_&foodtype;
set & outlib.._ pred_unc_&foodtype;
run;

Input datasets for the INDIVINT macro are created.  The datasets _pred_unc_&foodtype and _param_unc_&foodtype are created in the MIXTRAN run.  The syntax “if (_n_=1) then set” is a way of doing a 1-to-many merge, where the dataset _param_unc_&foodtype is merged with each row of the _pred_unc_&foodtype dataset.

data parsubj1rec;

  merge parampred &datam(keep = &subject &var12);

  by &subject;

The INDIVINT macro requires a dataset with only one line per person called parsubj1rec.  This is defined by the macro variable datam in the BRR212 macro, and the response variables are defined by the macro variable var12.

%DISTRIB (seed= 0 , nsim_mc= 100 ,modeltype=&modeltype, pred=& outlib.._ pred_&foodtype, param= & outlib.._ param_&foodtype, outlib=&outlib, subject=&subject,food=&foodtype);

The DISTRIB macro is used to find the lambda for the transformation of the estimate of T (true usual intake) that will be used in the regression model.

title&titles; 

title %eval (&titles+ 1 );  

This code clears the titles created in the DISTRIB macro and resets the title for the example.

data _dist;

set & outlib..d escript_&foodtype._w0304_0;

keep  tpercentile1-tpercentile99;

run;

The dataset descript_&foodtype_w0304_0 is from the DISTRIB macro.  The percentiles are saved from this dataset.

data _lambda (keep=lambda sse);
set _dist;
  array Z  (*) Z1-Z99;
  array T  (*) T1-T99;
  array CP (*) CP1-CP99;
  array P  (*) tpercentile1-tpercentile99;

  do i= 1 to dim(z);
    Z(i) = probit((i - 0.375 ) / 99.25 );
    end;

  MZ  = mean(of Z1-Z99);
  SZZ = CSS(of Z1-Z99);
  best_SSE    = SZZ;
  best_lambda = 0 ;
  do lambda = 0 to 1 by 0.01 ;
    do i = 1 to 99 ;
     if (lambda = 0 ) then T(i) = log(P(i));
     else T(i) = (P(i)**lambda - 1 ) / lambda;
      end;

    MT  = mean(of T1-T99);
    STT = CSS(of T1-T99);
    do i = 1 to 99 ;
      CP(i) = (T(i) - MT) * (Z(i) - MZ);
      end;

    STZ = sum(of CP1-CP99);
    SSE = SZZ - ((STZ** 2 )/STT);
    if (SSE < best_SSE) then do;
      best_SSE     = SSE;
      best_lambda  = lambda;
      end;
    end;

  lambda = best_lambda;
  sse    = best_sse;
  call symput( 'lamt' ,lambda);
run;

This code finds the best lambda to transform the estimate of T. This is done by finding the lambda that minimizes the error sum of squares.  Z is the inverse of the normal cumulative distribution function (probit).  At the end of the code, the best lambda is saved as a macro variable called lamt.

proc univariate data=&data noprint;
  where &response> 0 ;
  var &response;
  output out=outmin_amt min=min_amt;
run;

The minimum amount consumed on a consumption day is calculated and saved in a dataset called outmin_amt.

data paramsubj1rec;
  if _n_= 1 then set outmin_amt;
  set parsubj1rec;
  lamt = &lamt;
  array zero( 2 ) &var12;
  do i= 1 to 2 ;
  if zero[i]= 0 then zero[i]=(min_amt/ 2 );
  end;
run;

proc sort; by &subject; run;

The minimum amount is added to the dataset.  Zero values are set to ½ of the mimimum amount.  The dataset is sorted by subject.  The variable lamt is created from the macro value lamt that was created previously.

%indivint(model12= 2 , subj1recdata=paramsubj1rec, recid=&subject, r24vars=&var12, min_amt=min_amt, var_u1=p_var_u1, var_u2=a_var_u2, cov_u1u2=cov_u1u2, var_e=a_var_e, lambda=a_lambda, xbeta1=x1b1, xbeta2=x2b2, boxcox_t_lamt=y, lamt=lamt, dencalc=y, denopt=y, u1nlmix=, u2nlmix=, titles=&titles, notesprt=y);

The INDIVINT macro is called.  Because we are fitting a two-part model, the macro variable model12 is set equal to 2. This macro uses the dataset with 1 record per person (paramsubj1rec), and specifies the two 24-hour recall variables (var12).  The minimum amount that was saved in the paramsubj1rec dataset is specified for the min_amt macro variable, as are the parameters that are saved in that dataset, including p_var_u1, a_var_u2, cov_u1u2, a_var_e, and a_lambda. The predicted values x1b1 and x2b2 from the MIXTRAN run are specified. The macro variable boxcox_t_lamt is set to y for ‘Yes’, indicating that a Box-Cox parameter is specified.  This parameter is lamt, which was defined in the paramsubj1rec dataset in the last step.  The macro variables dencalc, denopt, and notesprt are set to y for ‘Yes’, indicating that denominator integration is perfomed, denominator optimization is performed, and notes are printed in the SAS log, respectively.

data subj1recres;
  merge &datam _resdata(keep = &subject indusint);
  by &subject;
run; 

The dataset _resdata is created in the INDIVINT macro.  It is merged with the dataset adultf.

proc reg data=subj1recres outest=regout;
  freq w0304_0;
  model &regresponse = indusint &covars_reg;
title %eval (&titles+ 1 ) "Regression Model:  Adjusted for Measurement Error" ;

The regression model is run using the expected usual intake from the INDIVINT macro, called indusint.  Because boxcox_t_lamt was set to y, this variable was already transformed in the INDIVINT macro, using the value of lamt in a Box-Cox transformation. Title3 is reset.

data diethealthout0;
  set regout(rename = (indusint=rcbeta));
  lamt = &lamt;
  t0boxcox = (&t0**lamt - 1 )/lamt;
  t1boxcox = (&t1**lamt - 1 )/lamt;
  tdiffrcbeta = (t1boxcox - t0boxcox)*rcbeta;
run; 

The output dataset from the regression above (regout) is set, and the parameter for the expected true usual intake is renamed to rcbeta.  In addition, lamt is saved, and the variables t0boxcox and t1boxcox and their difference multipled by rcbeta are computed.  The macro variables t0 and t1 are set to a range that is relevant to the data, e.g., the 10th and 90th percentile.

proc datasets nolist; delete parampred parsubj1rec paramsubj1rec subj1recres _resdata  regout; run;

Unneeded datasets are deleted.

%do run= 1 %to 16 ;

This code starts a loop to run the 16 BRR runs.

%put ~~~~~~~~~~~~~~~~~~~ Run &run ~~~~~~~~~~~~~~~~~~~~;

The run number is printed to the log.

%MIXTRAN (data=&data, response=&response, foodtype= &foodtype, subject=&subject, repeat=&repeat,     covars_prob=&covars_prob, covars_amt=&covars_amt, outlib=&outlib, modeltype=&modeltype, lambda=&lambda, replicate_var=w0304_&run, seq=&seq, weekend=&weekend, vargroup=&vargroup, numvargroups= &numvargroups, subgroup=&subgroup,  start_val1=&start_val1, start_val2=&start_val2, start_val3=&start_val3, vcontrol=&vcontrol, nloptions=&nloptions, titles= &titles, printlevel=&printlevel)

Within the BRR212 macro, the MIXTRAN macro is called.  All of the variables preceded by & will be defined by the BRR212 macro call.  The only variable without an & is the replicate_var macro variable; it is set to w0304_&run where &run equals 1 to 16.

data parampred(keep = &subject a_var_u2 a_var_e a_lambda x2b2);

if (_n_ = 1 ) then set & outlib.._ param_unc_&foodtype;  set & outlib.._ pred_unc_&foodtype;

run;

Input datasets for the INDIVINT macro are created.  The datasets _pred_unc_&foodtype and _param_unc_&foodtype are created in the MIXTRAN run.  The syntax “if (_n_=1) then set” is a way of doing a 1-to-many merge, where the dataset _param_unc_&foodtype is merged with each row of the _pred_unc_&foodtype dataset.

data parsubj1rec;

  merge parampred &datam(keep = &subject &var12);

  by &subject;

The INDIVINT macro requires a dataset with only one line per person called parsubj1rec.  This is defined by the macro variable datam in the BRR212 macro, and the response variables are defined by the macro variable var12.

proc univariate data=&data noprint;
  where &response> 0 ;
  var &response;
  output out=outmin_amt min=min_amt;
run;

The minimum amount consumed on a consumption day is calculated and saved in a dataset called outmin_amt.

data paramsubj1rec;
  if _n_= 1 then set outmin_amt;
  set parsubj1rec;
  array zero( 2 ) &var12;
  do i= 1 to 2 ;
  if zero[i]= 0 then zero[i]=(min_amt/ 2 );
  end;
  lamt = &lamt;
run;

The minimum amount is added to the dataset.  Zero values are set to ½ of the mimimum amount.  The dataset is sorted by subject.  The variable lamt is created from the macro value lamt that was created previously.

%indivint(model12= 2 , subj1recdata=paramsubj1rec, recid=&subject, r24vars=&var12, min_amt=min_amt, var_u1=, var_u2=a_var_u2, cov_u1u2=, var_e=a_var_e, lambda=a_lambda, xbeta1=, xbeta2=x2b2, boxcox_t_lamt=y, lamt=lamt, dencalc=y, denopt=y, u1nlmix=, u2nlmix=, titles=&titles, notesprt=y);

The INDIVINT macro is called.  Because we are fitting a one-part model, the macro variable model12 is set equal to 2. This macro uses the dataset with 1 record per person (paramsubj1rec), and specifies the two 24-hour recall variables (var12).  The minimum amount that was saved in the paramsubj1rec dataset is specified for the min_amt macro variable, as are the parameters that are saved in that dataset, including p_var_u1, a_var_u2, cov_u1u2, a_var_e, and a_lambda. The predicted values x1b1and x2b2 from the MIXTRAN run are specified. The macro variable boxcox_t_lamt is set to y for ‘Yes’, indicating that a Box-Cox parameter is specified.  This parameter is lamt, which was defined in the paramsubj1rec dataset in the last step.  The macro variables dencalc, denopt, and notesprt are set to y for ‘Yes’, indicating that denominator integration is perfomed, denominator optimization is performed, and notes are printed in the SAS log, respectively.

data subj1recres;
  merge &datam _resdata(keep = &subject indusint);
  by &subject;
run; 

The dataset _resdata is created in the INDIVINT macro.  It is merged with the dataset adultf.

proc reg data=subj1recres outest=regout;
  freq w0304_&run;
  model &regresponse = indusint &covars_reg;
  title %eval (&titles+ 1 ) "Regression Model:  Adjusted for Measurement Error" ;

The regression model is run using the expected usual intake from the INDIVINT macro, called indusint.  Because boxcox_t_lamt was set to y, this variable was already transformed in the INDIVINT macro, using the value of lamt in a Box-Cox transformation. Title3 is reset.

data diethealthout;
  set regout(rename = (indusint=rcbeta));
  lamt = &lamt;
  t0boxcox = (&t0**lamt - 1 )/lamt;
  t1boxcox = (&t1**lamt - 1 )/lamt;
  tdiffrcbeta = (t1boxcox - t0boxcox)*rcbeta;
run; 

The output dataset from the regression above (regout) is set, and the parameter for the expected true usual intake is renamed to rcbeta.  In addition, lamt is saved, and the variables t0boxcox and t1boxcox and their difference multipled by rcbeta are computed.  The macro variables t0 and t1 are set to a range that is relevant to the data, e.g., the 10th and 90th percentile.

proc append base=brr_runs data=diethealthout;
run;

The BRR datasets are appended into a dataset called brr_runs.

proc datasets nolist; delete parampred parsubj1rec paramsubj1rec subj1recres _resdata  regout diethealthout; run;

After appending the information to brr_runs, diethealthout and the other datasets created in the BRR runs can be deleted.

%end ;

The BRR runs end.

data _null_;
format list listd listsum $255. ;
%let I=1;
%let varamtu= %upcase (&covars_reg);
%do %until ( %qscan (&varamtu,&I, %str ( ))= %str ());
  %let varb&I= %qscan (&varamtu,&I, %str ( ));
      num= %eval (&i+ 1 );
      varA= strip( "&&varb&i." );
      list =  trim(list)|| ' ' || 'b' ||trim(varA);
      listd =  trim(listd)|| ' ' || 'db' ||trim(varA);
      listsum =  trim(listsum)|| ' ' || 'sum_db' ||trim(varA);
%let I= %eval (&I+1);
%end ;
%let cnt= %eval (&I-1);
%if &covars_reg= %str () %then %let cnt=0;
call symput( 'list' ,list);
call symput( 'listd' ,listd);
call symput( 'listsum' ,listsum);
run;

This code creates macro calls that are a list of the variables with the prefix ‘b’ and ‘db’.

data parambrr;
set brr_runs;
rename rcbeta=brcbeta tdiffrcbeta=btdiffrcbeta intercept=bintercept;
array covs (&cnt) &covars_reg;
array newname (&cnt) &list;
do i= 1 to &cnt;
  newname[i]=covs[i];
 end;
data parambrr;
set parambrr;
keep bintercept brcbeta btdiffrcbeta &list run;
run;

The brr_runs dataset is set, and a new dataset is created.  This dataset renames the BRR parameters to have a prefix of ‘b’ using a rename statement and the macro calls created above.

data ss;
  if _n_= 1 then set diethealthout0;
  set parambrr;
  array bvar (*) bintercept brcbeta btdiffrcbeta &list;
  array varo  (*)intercept rcbeta tdiffrcbeta &covars_reg;
  array dsqr (*) dbintercept dbrcbeta dbtdiffrcbeta &listd;
  do i= 1 to dim(bvar);
   dsqr[i]=(bvar[i]-varo[i])** 2 ;
  end;
run;

The estimates from the first run are merged with the BRR estimates, and the squared difference is computed for each BRR run.

proc means data=ss sum;
 var  dbintercept dbrcbeta dbtdiffrcbeta &listd;
 output out=sums sum=sum_dbintercept sum_dbrcbeta  sum_dbtdiffrcbeta &listsum;
run;

The sum of squares is computed.

data brr;
 set sums;
 array sumt (*) sum_dbintercept sum_dbrcbeta  sum_dbtdiffrcbeta &listsum;
 array se  (*) intercept rcbeta tdiffrcbeta &covars_reg;
 do j= 1 to dim(sumt);
 se[j]=sqrt((sumt[j])/( 16 * .49 ));
 end;
 keep intercept rcbeta tdiffrcbeta &covars_reg;
run;

The standard errors are computed.

data wparms;
 set diethealthout0 brr;
keep intercept rcbeta tdiffrcbeta &covars_reg;

The point estimates are merged with the standard errors.

proc transpose data=wparms out=wparmst;
run;

The data are transposed.

data &final;
  format pvalue 6.4 ;
  set wparmst(rename=(col1=estimate col2=brrse));
  waldchisq = (estimate/brrse)** 2 ;
  pvalue = 1 - probchi(waldchisq, 1 );
  lower95lim = estimate - 1.96 *brrse;
  upper95lim = estimate + 1.96 *brrse;
  drop _label_;
run; 

The final dataset is created, and p-values and 95% confidence intervals are computed.

proc print;
run;

The final dataset is printed.

%mend BRR212;

The end of the BRR212 macro is indicated.