/*****************************************************************************/
/*** The following files are required for this example:                    ***/
/***   brr_pvalue_ci_macro_v1.1.sas                                        ***/
/***   means.hei2010.nonsmok.smok.ag20p.xpt.                               ***/
/***                                                                       ***/
/*** This example includes 5 programs.  Programs 1a, 1b, 1c can be used to ***/
/*** fit a multivariate measurement error model, calculate mean usual      ***/
/*** Healthy Eating Index–2010 (HEI–2010) scores, and perform t-tests      ***/
/*** comparing mean usual HEI–2010 scores using NHANES 2003-2004 data.     ***/
/*** Programs 1d and 1e can be used to extend the analysis to also         ***/
/*** calculate percentiles and cutpoint probablities for both the HEI-2010 ***/
/*** component scores and the HEI-2010 densities, conditional on           ***/
/*** categories of HEI-2010 total scores.                                  ***/
/***                                                                       ***/
/*** This SAS program illustrates the use of the BRR_PVALUE_CI macro to    ***/
/*** perform balanced repeated replication (BRR) variance estimation and   ***/
/*** t-tests comparing mean usual HEI–2010 scores for nonsmokers versus    ***/
/*** smokers.  The input data set  includes mean usual HEI-2010 scores     ***/
/*** calculated using the original weights (i.e. weights for replicate 0)  ***/
/*** and using the 14 valid BRR replicates (i.e. 1, 2, ..., 13, 14, 16)    ***/
/*** and includes the differences between these means for smokers and      ***/
/*** nonsmokers.                                                           ***/
/*****************************************************************************/


/*****************************************/
/*** Specify a title for this example. ***/
/*****************************************/

title1 "Fit Multivariate Measurement Error Model Using MCMC with 24-Hour Recall as Main Instrument";
title2 "Calculate Mean Usual HEI–2010 Scores and Perform T-Tests Comparing Mean Usual HEI–2010 Scores";
title3 "Calculate Conditional Distributions and Cutpoint Probabilities for HEI-2010 Densities and Scores";


/*****************************************************************************/
/*** The path must be modified in the following lines to specify the       ***/
/*** location of the SAS macro and the SAS data file.                      ***/
/*****************************************************************************/

%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";

*** Input data ***;
filename tranmean "&home/multivar_surveillance_example1_mcmc_nnc_main24hr/data/means.hei2010.nonsmok.smok.ag20p.xpt";

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



/*****************************************************************************/
/*** Import the data file 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, ..., B       ***/
/*** including parameter estimates from replicate analyses performed for   ***/
/*** BRR variance estimation.                                              ***/
/*****************************************************************************/

proc cimport data=means_hei2010_nonsmok_smok_ag20p infile=tranmean;
run;

data origwt_paramest brr_paramest;
  set means_hei2010_nonsmok_smok_ag20p;
  if replicate=0 then output origwt_paramest;
  else output brr_paramest;
run;



/*****************************************************************************/
/*** Use the BRR_PVALUE_CI macro to perform t-tests comparing mean usual   ***/
/*** HEI–2010 scores for nonsmokers versus smokers                         ***/
/*****************************************************************************/

/*****************************************************************************/
/*                                                                           */
/* 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, ..., B */
/* 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.                                                       */
/*                                                                           */
/*****************************************************************************/
/*****************************************************************************/

%let meansmok = mean_smok_hei_tot_fruit
                mean_smok_hei_wholefrt
                mean_smok_hei_tot_veg
                mean_smok_hei_beangrn
                mean_smok_hei_g_whl
                mean_smok_hei_d_total
                mean_smok_hei_tot_prot
                mean_smok_hei_seaplant
                mean_smok_hei_fat_acids_ratio
                mean_smok_hei_g_nwhl
                mean_smok_hei_tsodi
                mean_smok_hei_emptycal
                mean_smok_hei_total
                ;

%let meannonsmok = mean_nonsmok_hei_tot_fruit
                   mean_nonsmok_hei_wholefrt
                   mean_nonsmok_hei_tot_veg
                   mean_nonsmok_hei_beangrn
                   mean_nonsmok_hei_g_whl
                   mean_nonsmok_hei_d_total
                   mean_nonsmok_hei_tot_prot
                   mean_nonsmok_hei_seaplant
                   mean_nonsmok_hei_fat_acids_ratio
                   mean_nonsmok_hei_g_nwhl
                   mean_nonsmok_hei_tsodi
                   mean_nonsmok_hei_emptycal
                   mean_nonsmok_hei_total
                   ;

%let param_names = &meansmok &meannonsmok hei_mean_diff1 - hei_mean_diff13;


/*****************************************************************************/
/*** Delete replicates 12 and 14 due to unusable parameter estimates of    ***/
/*** covariance matrix parameters in the multivariate measurement error    ***/
/*** model for smokers.  Fourteen replicates are available for the         ***/
/*** analysis, and the BRR_PVALUE_CI macro call specifies the value 14 for ***/
/*** the "set_degrees_freedom" macro parameter.                            ***/
/***                                                                       ***/
/*** If all 16 replicates had yielded usable results, the degrees of       ***/
/*** freedom would have been limited to 15 (i.e. the number of sampled     ***/
/*** primary sampling units (PSUs) minus the number of strata is 30-15=15  ***/
/*** for NHANES 2003-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.                                        ***/
/*****************************************************************************/

data del_brr_paramest keep_brr_paramest;
  set brr_paramest;
  if (_n_ = 12) or (_n_ = 14) then output del_brr_paramest;
  else output keep_brr_paramest;
run;

%brr_pvalue_ci(data_orig_estimates_1rec = origwt_paramest,
               data_brr_estimates_brecs = keep_brr_paramest,
               param_estimate_names     = &param_names,
               set_f_method_fay         = 0.3,
               set_confidence_level     = 0.99,
               set_degrees_freedom      = 14,
               print                    = y,
               titles                   = 3
               );


/***************************************************/
/*** Save the output data set from BRR_PVALUE_CI ***/
/***************************************************/

data outlib.brr_pvalue_ci_out_hei_1c;
  set brr_pvalue_ci_out;
run;
