/****************************************************************************/
/****************************************************************************/
/* The following files are required to run this example:                    */
/*   brr_pvalue_ci_macro_v1.1.sas                                           */
/*   param_unc_add_sug_male0-add_sug_male32 (parameter estimates from the   */
/*     univariate measurement error model fit for the original data and for */
/*     each of the replicates)                                              */
/*   descript_add_sug_male0-add_sug_male32 (distribution estimates for      */
/*     the original data and for each of the replicates)                    */
/****************************************************************************/
/*                                                                          */
/****************************************************************************/
/* This example fits a univariate measurement model and estimates the       */
/* distribution of usual intake for a dietary component consumed almost     */
/* daily.  For this analysis, the example data includes males age 9+ from   */
/* the 2001-2004 public use NHANES, and the 24-hour dietary recall is the   */
/* main instrument.  The dietary component is added sugars.  In addition    */
/* recommended intakes (EAR and TUL) for several components have been added */
/* to the data for each individual.  The example data include repeated      */
/* intake measurements from 24-hour dietary recalls.                        */
/*                                                                          */
/* This analysis uses balanced repeated replication (BRR) variance          */
/* estimation, so the univariate measurement error model and distribution   */
/* of usual intake use the original weight and 32 BRR weights.  This step   */
/* of the example uses estimates calculated from the previous step and uses */
/* the BRR_PVALUE_CI macro to perform balanced repeated replication (BRR)   */
/* estimation of the variance.                                              */
/****************************************************************************/



title1 "Fit Univariate Measurement Error Model Using MLE with 24-Hour Recall as Main Instrument";
title2 "Estimate Distribution of Usual Intake for a Dietary Component Consumed Nearly Every Day";
title3 "Using 2001-2004 NHANES Public Use Data for Males Age 9+";
title4 "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/univar_surveillance_example1_mle_main24hr/outdata";



******************************************************************************;
***              Confidence Intervals of Parameter Estimates               ***;
******************************************************************************;

******************************************************************************;
*** 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 add_sug_male_0;
  set outlib._param_unc_add_sug_male0;
  drop FreqName min_amt numvargrps weekendflag;  *** variables for which SEs are not wanted ***;
run;

data add_sug_male_1to32;
  set outlib._param_unc_add_sug_male1-outlib._param_unc_add_sug_male32;
  drop FreqName min_amt numvargrps weekendflag;
run;


*********************************************************;
*** Put the names of the parameter variables into a   ***;
*** macro variable for use by the BRR_PVALUE_CI macro ***;
*********************************************************;

proc sql noprint;
  select name
  into :parm_vars separated by ' '
  from dictionary.columns
  where libname='WORK' and memname = 'ADD_SUG_MALE_0'
  order by name;
quit;



/*****************************************************************************/
/*****************************************************************************/
/*                                                                           */
/* 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.                                         ***;
******************************************************************************;

%brr_pvalue_ci(data_orig_estimates_1rec = add_sug_male_0,
               data_brr_estimates_brecs = add_sug_male_1to32,
               param_estimate_names     = &parm_vars,
               set_f_method_fay         = 0.3,
               set_confidence_level     = 0.95,
               set_degrees_freedom      = 30,
               print                    = y,
               titles                   = 4
               );


***************************************************************************************;
*** Save the output data set of standard errors, p-values, and confidence intervals ***;
***************************************************************************************;

data outlib.brr_pvalue_ci_out_param_male;
  set brr_pvalue_ci_out;
run;



******************************************************************************;
***           Confidence Intervals of Distribution Estimates               ***;
******************************************************************************;

******************************************************************************;
*** Calculate the confidence intervals for the percentiles etc. output     ***;
*** into the "descript" data files by the DISTRIB macro.  In this example  ***;
*** we used a subgroup named agegrp, so unlike the parameter files the     ***;
*** output data sets contain more than one record per replicate.  To       ***;
*** handle this in the BRR_PVALUE_CI macro  we will need to process the    ***;
*** data one subgroup at a time. To accomplish this a utility macro        ***;
*** BRR_REPLICATE_LOOP will be used.                                       ***;
***                                                                        ***;
*** NOTE: If there are no subgroups in the "descript" data sets it is not  ***;
*** necessary to use the subgroupforCI macro.  Simply follow the example   ***;
*** code for calculating the parameter CIs, substituting the "descript"    ***;
*** data for the "_param" data sets                                        ***;
******************************************************************************;

data add_sug_maleD_0;
  set outlib.descript_add_sug_male_rndw0;
  drop numsubjects;
run;

data add_sug_maleD_1to32;
  set outlib.descript_add_sug_male_rndw1-outlib.descript_add_sug_male_rndw32;
  drop numsubjects;
run;



****************************************************;
*** Get the number of categories in the subgroup ***;
****************************************************;

proc sql noprint;
   select count(*) into :numgrps from add_sug_maleD_0;
quit;



**************************************************************;
*** Create the utility macro BRR_REPLICATE_LOOP to prepare ***;
*** the data and call the macro BRR_PVALUE_CI              ***;
**************************************************************;

%macro brr_replicate_loop(subgrp=);

  /*************************************************************************/
  /* The macro BRR_REPLICATE_LOOP selects the data into separate data sets */
  /* per subgroup category, identifies the variables to be used in the     */
  /* BRR_PVALUE_CI macro, and calls that macro. Additionally all the CI    */
  /* data sets output by each iteration are concatenated, the subgroup     */
  /* added and the final data set is saved to disk.                        */
  /*                                                                       */
  /* The syntax for calling the macro is:                                  */
  /*   %brr_replicate_loop(subgrp=);                                       */
  /*                                                                       */
  /* where:                                                                */
  /*   subgroup is used to identify the name of the subgroup variable      */
  /*                                                                       */
  /*  For example:                                                         */
  /*    %brr_replicate_loop(subgrp=agegrp);                                */
  /*************************************************************************/

  **********************************************************************;
  *** Create an identifying variable based on the number of subgroup ***;
  *** categories to make it easy to select the records.              ***;
  *** Save this on the base data set and create a data set with only ***;
  *** the subgroup information and the record_id for later use.      ***;
  **********************************************************************;

  *********************;
  *** Base data set ***;
  *********************;

  data add_sug_maleD_0 subgrps(keep = record_id &subgrp);
    set add_sug_maleD_0;
      record_id=_n_;
    output add_sug_maleD_0 subgrps;
  run;



  *********************************************************************************;
  *** Replicate data sets                                                       ***;
  *** Add the record_id associated with each subgrp value to the replicate data ***;
  *********************************************************************************;

  proc sort data=subgrps;
    by &subgrp;
  run;

  proc sort data=add_sug_maleD_1to32;
    by &subgrp;
  run;

  data add_sug_maleD_1to32 ;
    merge subgrps add_sug_maleD_1to32;
    by &subgrp;
  run;



  ***********************************************************************************;
  *** Create data sets for the CI macro, each data set can have only one subgroup ***;
  ***********************************************************************************;

  %do a = 1 %to &numgrps;

    *********************;
    *** Base data set ***;
    *********************;

    data base (drop=record_id &subgrp);
      set add_sug_maleD_0;
      if &a = record_id;
    run;


    ***************************;
    *** Replicate data sets ***;
    ***************************;

    data replicate (drop=record_id &subgrp);
      set add_sug_maleD_1to32;
      if &a = record_id;
    run;


    *********************************************************************************;
    *** Get the names of the variables to be processed in the BRR_PVALUE_CI macro ***;
    *********************************************************************************;

    %if &a = 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;


    ************************************;
    *** Call the BRR_PVALUE_CI macro ***;
    ************************************;

    title5 "Output for Subgroup &a";

    %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 sets ***;
    ***************************************************;

    data brr_pvalue_ci_out_pct_male&a;
      set brr_pvalue_ci_out;
      record_id = &a;
    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 set of CIs for the "descript" data, add the       ***;
    *** subgroup values appropriately, and save the final data set ***;
    ******************************************************************;

    %if &a = &numgrps %then %do ;

      data brr_pvalue_ci_out_pct_male_all ;
        set brr_pvalue_ci_out_pct_male1-brr_pvalue_ci_out_pct_male&a;
      run;

      proc sort data=brr_pvalue_ci_out_pct_male_all;
        by record_id;
      run;

      proc sort data=subgrps;
        by record_id;
      run;

      data outlib.brr_pvalue_ci_out_pct_male_all (drop=record_id);
        merge brr_pvalue_ci_out_pct_male_all subgrps;
        by record_id;
      run;

    %end;


  %end; /* of iterations for each subgroup */


%mend; /* of BRR_REPLICATE_LOOP macro */



*************************************************;
*** Call the utility macro brr_replicate_loop ***;
*************************************************;

%brr_replicate_loop(subgrp=agegrp);

