options mergenoby=nowarn;

/***************************************************************************/
/* The following files are required to run this example:                   */
/*   brr_pvalue_ci_macro_v1.1.sas                                          */
/*   pctl_b_sfat_energy0.sas7bdat - pctl_b_sfat_energy32.sas7bdat          */
/***************************************************************************/
/*                                                                         */
/***************************************************************************/
/* This example analysis fits univariate and bivariate measurement error   */
/* models for two dietary components that are consumed almost daily.  It   */
/* also estimates the percentiles of the population distribution of the    */
/* ratio of usual intakes of two dietary components.  For this analysis,   */
/* the example data include data from the public use 2001-2004 NHANES,     */
/* children aged 1-8, males aged 9+ years, and females aged 9+ years, and  */
/* the 24-hour dietary recall is the main instrument.  The two dietary     */
/* components are saturated fat and energy.  The example data include      */
/* repeated intake measurements from 24-hour dietary recalls.              */
/*                                                                         */
/* This analysis uses balanced repeated replication (BRR) variance         */
/* estimation, so the bivariate measurement error model and distribution   */
/* of usual intake use the original weight and 32 BRR weights.
/***************************************************************************/



title1 "Fit Bivariate Measurement Error Model Using MLE with 24-Hour Recall as Main Instrument";
title2 "Estimate Percentiles of the Population Distribution of the Ratio of Usual Intakes of Two Dietary Components";
title3 "Balanced Repeated Replication (BRR) Variance Estimation";



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

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

*** Include the required macros ***;
%include "&home/include.files.macros/brr_pvalue_ci_macro_v1.1.sas";

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



***************;
*** Formats ***;
***************;

proc format;
  value subfmt
    1 = "Children 1-3"
    2 = "Children 4-8"
    3 = "Males 9-13"
    4 = "Males 14-18"
    5 = "Males 19-30"
    6 = "Males 31-50"
    7 = "Males 51-70"
    8 = "Males 71+"
    9 = "Males 19+"
   10 = "Females 9-13"
   11 = "Females 14-18"
   12 = "Females 19-30"
   13 = "Females 31-50"
   14 = "Females 51-70"
   15 = "Females 71+"
   16 = "Females 19+"
   17 = "All 1+";
run;



*****************************************************************************;
*** Import the data sets of intermediate results needed for this example  ***;
*** and create a data set with one record including parameter estimates   ***;
*** from the original data and a data set with records 1, 2, ..., 32      ***;
*** including parameter estimates from replicate analyses performed for   ***;
*** BRR variance estimation.                                              ***;
*****************************************************************************;

data pctl_b_sfat_energy0(keep = mean subpopulation Pctile5 Pctile10 Pctile25 Pctile50
                                Pctile75 Pctile90 Pctile95 Prob1-Prob3);
  set outlib.pctl_b_sfat_energy0;
run;

data pctl_b_sfat_energy1to32(keep = mean subpopulation Pctile5 Pctile10 Pctile25 Pctile50
                                    Pctile75 Pctile90 Pctile95 Prob1-Prob3);
  set outlib.pctl_b_sfat_energy1-outlib.pctl_b_sfat_energy32;
run;



****************************************************;
*** Get the number of categories in the subgroup ***;
****************************************************;
proc sql noprint;
  select count(*) into :num_subgroups from pctl_b_sfat_energy0;
quit;



******************************************************************************;
*** The utility macro BRR_REPLICATE_LOOP selects the data into separate    ***;
*** data sets for each subgroup category, identifies the variables to be   ***;
*** used in the BRR_PVALUE_CI macro, and calls that macro. Additionally    ***;
*** all the data sets output by each iteration of the macro are            ***;
*** concatenated and saved into one data set.                              ***;
***                                                                        ***;
*** The syntax for calling the brr_replicate_loop macro is:                ***;
***                                                                        ***;
*** %brr_replicate_loop(subgroup = );                                      ***;
***                                                                        ***;
*** where                                                                  ***;
***                                                                        ***;
*** subgroup  =    specifies the name of the subgroup variable             ***;
******************************************************************************;

%macro brr_replicate_loop(subgroup=);

  ******************************************************************;
  *** Create an identifying variable based on the number of      ***;
  *** subgroup cateogories to make it easy to select the records ***;
  ******************************************************************;
  
  data pctl_b_sfat_energy0 subgroups(keep=record_id &subgroup);
    set pctl_b_sfat_energy0;
    record_id = _n_;
    output pctl_b_sfat_energy0 subgroups;
  run;



  *********************************************************************************;
  *** Add the record_id associated with each subgrp value to the replicate data ***;
  *********************************************************************************;
  
  proc sort data=subgroups;
    by &subgroup;
  run;
  
  proc sort data=pctl_b_sfat_energy1to32;
    by &subgroup;
  run;
  
  data pctl_b_sfat_energy1to32 ;
    merge subgroups pctl_b_sfat_energy1to32;
    by &subgroup;
  run;
  
  
  
  ****************************************************************;
  *** Loop through all the categories in the subgroup variable ***;
  ****************************************************************;

  %do i = 1 %to &num_subgroups;
    
    ******************************************************************;
    *** Create a data sets for the BRR_PVALUE_CI with one subgroup ***;
    ******************************************************************;
      
    data base(drop=record_id &subgroup);
      set pctl_b_sfat_energy0;
      if &i = record_id;
    run;
    
    data replicate(drop=record_id &subgroup);
      set pctl_b_sfat_energy1to32;
      if &i = record_id;
    run; 
    
    
    
    *********************************************************************************;
    *** Get the names of the variables to be processed in the BRR_PVALUE_CI macro ***;
    *********************************************************************************;
    
    %if &i = 1 %then %do;
      proc sql noprint;
        select name
        into :desc_vars separated by ' '
        from dictionary.columns
        where libname='WORK' and memname = 'BASE'
        order by name;
      quit;
    %end;



    /*****************************************************************************/
    /*****************************************************************************/
    /*                                                                           */
    /* Description of the BRR_PVALUE_CI macro                                    */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /*                                                                           */
    /* The BRR_PVALUE_CI macro performs balanced repeated replication (BRR)      */
    /* variance estimation, calculates a p-value for each user specified         */
    /* parameter to test the null hypothesis that the parameter value is zero,   */
    /* and calculates a confidence interval for each user specified parameter.   */
    /* The macro produces an output data set including the parameter estimates,  */
    /* the BRR standard error estimates, the p-value used to test the null       */
    /* hypothesis that the parameter value is zero versus the alternative        */
    /* hypothesis that the parameter value is not zero, and the lower and upper  */
    /* confidence limits.                                                        */
    /*                                                                           */
    /* The macro input includes a data set with one record including parameter   */
    /* estimates from the original data and a data set with records 1, 2, ...,   */
    /* including parameter estimates from replicate analyses performed for BRR   */
    /* variance estimation.  Optionally, Fay's BRR method can be specified.      */
    /*                                                                           */
    /* Reference:                                                                */
    /*                                                                           */
    /*   Korn EL, Graubard BI. Analysis of Health Surveys, John Wiley & Sons,    */
    /*   Inc., New York, 1999.                                                   */
    /*                                                                           */
    /*                                                                           */
    /* The syntax for calling the BRR_PVALUE_CI macro is:                        */
    /*                                                                           */
    /* %brr_pvalue_ci(data_orig_estimates_1rec =,                                */
    /*                data_brr_estimates_brecs =,                                */
    /*                param_estimate_names     =,                                */
    /*                set_f_method_fay         =,                                */
    /*                set_confidence_level     =,                                */
    /*                set_degrees_freedom      =,                                */
    /*                print                    =,                                */
    /*                titles                   =                                 */
    /*                );                                                         */
    /*                                                                           */
    /*  where                                                                    */
    /*                                                                           */
    /*  "data_orig_estimates_1rec"    Specifies an input data set that includes  */
    /*                                one record with parameter estimates from   */
    /*                                the original data.                         */
    /*                                                                           */
    /*  "data_brr_estimates_brecs"    Specifies an input data set that includes  */
    /*                                records 1, 2, ..., B with parameter        */
    /*                                estimates from replicate analyses          */
    /*                                performed for balanced repeated            */
    /*                                replication (BRR) variance estimation.     */
    /*                                For example, record 1 includes the         */
    /*                                parameter estimates obtained from          */
    /*                                replicate 1, and record B includes the     */
    /*                                parameter estimates obtained from          */
    /*                                replicate B, i.e. the final replicate.     */
    /*                                                                           */
    /*  "param_estimate_names"        Specifies a list of parameter estimates    */
    /*                                that are included as variables in the      */
    /*                                input data sets specified for the          */
    /*                                "data_orig_estimates_1rec" and             */
    /*                                "data_brr_estimates_brecs" macro           */
    /*                                parameters.                                */
    /*                                                                           */
    /*  "set_f_method_fay"            Specifies a value f that will be used to   */
    /*                                define an adjustment factor 1/[(1-f)**2]   */
    /*                                used to perform an appropriate adjustment  */
    /*                                if Fay's BRR method was used to create the */
    /*                                replicate weights.  The specified value    */
    /*                                should be nonnegative and less than 1.     */
    /*                                The default value is "0" so by default,    */
    /*                                the BRR calculations are based on the      */
    /*                                assumption that Fay's method was not used  */
    /*                                to create the replicate weights.           */
    /*                                                                           */
    /*  "set_confidence_level"        Specifies the confidence level.  The       */
    /*                                default value is "0.95".                   */
    /*                                                                           */
    /*  "set_degrees_freedom"         Specifies the degrees of freedom for the t */
    /*                                distribution and for the denominator of    */
    /*                                the F distribution.                        */
    /*                                                                           */
    /*  "print"                       If "print=n" or "print=N" then macro       */
    /*                                results are not printed.  The default      */
    /*                                value is "y".                              */
    /*                                                                           */
    /*  "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".          */
    /*                                                                           */
    /*****************************************************************************/
    /*                                                                           */
    /* Macro Output Data Set:  brr_pvalue_ci_out                                 */
    /*                                                                           */
    /* The BRR_PVALUE_CI macro produces a SAS data set "brr_pvalue_ci_out" that  */
    /* includes the following variables:                                         */
    /*               _name_, estimate, brr_se, pvalue_param0_probf_df_1_r,       */
    /*               lowerlim_t_df_r, and upperlim_t_df_r                        */
    /* where r is the value that was specified using the "set_degrees_freedom"   */
    /* macro parameter.  The variable "_name_" identifies each parameter, and    */
    /* the variable "estimate" includes the parameter estimates from the data    */
    /* set specified using the "data_orig_estimates_1rec" macro parameter.  The  */
    /* variable "brr_se" includes the BRR standard error estimates.  The         */
    /* variable "pvalue_param0_probf_df_1_r" is the p-value used to test the     */
    /* null hypothesis that the parameter value is zero versus the alternative   */
    /* hypothesis that the parameter value is not zero.  The p-value is          */
    /* calculated using the F distribution with 1 numerator degree of freedom    */
    /* and r denominator degrees of freedom.  An equivalent p-value could be     */
    /* obtained using the t distribution with r degrees of freedom.  The         */
    /* variables "lowerlim_t_df_r" and "upperlim_t_df_r" are the lower and upper */
    /* confidence limits which are calculated using the t distribution with r    */
    /* degrees of freedom.                                                       */
    /*                                                                           */
    /*****************************************************************************/
    /*****************************************************************************/

    ******************************************************************************;
    ***       Call the BRR_PVALUE_CI macro to calculate standard errors,       ***;
    ***                   p-values, and confidence intervals                   ***;
    ******************************************************************************;
    *** Note that the degrees of freedom is specified as 30 (i.e. the number   ***;
    *** of sampled primary sampling units (PSUs) minus the number of strata is ***;
    *** 60-30=30 for NHANES 2001-2004). The degrees of freedom associated with ***;
    *** replication variance estimators is discussed by Korn and Graubard      ***;
    *** (1999, Analysis of Health Surveys, pg. 34) and in the documentation    ***;
    *** for the SAS survey procedures.                                         ***;
    ******************************************************************************;

    title5 "Output for Subgroup = &i";

    %brr_pvalue_ci(data_orig_estimates_1rec = base,
                   data_brr_estimates_brecs = replicate,
                   param_estimate_names     = &desc_vars,
                   set_f_method_fay         = 0.3,
                   set_confidence_level     = 0.95,
                   set_degrees_freedom      = 30,
                   print                    = y,
                   titles                   = 5
                   );



    **************************************************;
    *** Keep each iteration of the output data set ***;
    **************************************************;
       
    data brr_pvalue_ci_out_pctl&i;
       set brr_pvalue_ci_out;
       record_id = &i;
    run;
    
    
    
    *********************************************************************;
    *** Delete data sets that will be recreated in the next iteration ***;
    *********************************************************************;
    
    proc datasets nolist lib=work;
      delete brr_pvalue_ci_out base replicate;
    quit;
    
    
    
    *****************************************************************;
    *** When all subgroups categories have been processed combine ***;
    *** the data sets of standard errors, p-values, and CIs for   ***;
    *** each subgroup and save the final data set                 ***;
    *****************************************************************;
    
    %if &i = &num_subgroups %then %do ;
    
      data brr_pvalue_ci_out_pctl_all;
        set brr_pvalue_ci_out_pctl1-brr_pvalue_ci_out_pctl&i;
      run;
   
      proc sort data=brr_pvalue_ci_out_pctl_all;
        by record_id;
      run;
      
      proc sort data=subgroups;
        by record_id;
      run;
  
      data outlib.brr_pvalue_ci_out_pctl_all (drop=record_id);
        merge brr_pvalue_ci_out_pctl_all subgroups;
        by record_id;
        format &subgroup subfmt.;
      run;
      
    %end;


  %end;


%mend;



%brr_replicate_loop(subgroup=Subpopulation);
