
/*****************************************************************************/
/*****************************************************************************/
/*                                                                           */
/*                            BRR_PVALUE_CI MACRO                            */
/*                                                                           */
/*****************************************************************************/
/*                         VERSION 1.1      1/19/2018                        */
/*                                                                           */
/* 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.                                                       */
/*                                                                           */
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/


%macro brr_pvalue_ci(data_orig_estimates_1rec =,
                     data_brr_estimates_brecs =,
                     param_estimate_names     =,
                     set_f_method_fay         = 0,
                     set_confidence_level     = 0.95,
                     set_degrees_freedom      =,
                     print                    = y,
                     titles                   = 0
                     );




  /**********************************************************************************/
  /*** Ensure user macro call does not overwrite default values with a null value ***/
  /**********************************************************************************/

  %if (&set_f_method_fay=%str()) %then %let set_f_method_fay = 0;
  %if (&set_confidence_level=%str()) %then %let set_confidence_level = 0.95;
  %if (&print=%str()) %then %let print = y;
  %if (&titles=%str()) %then %let titles = 0;


  /****************************************************/
  /*** Process macro input and perform calculations ***/
  /****************************************************/

  data _null_;
    set &data_orig_estimates_1rec;

    array estnames(*) &param_estimate_names;

    num_param_estimate_names = dim(estnames);

    call symputx("k_est_names", num_param_estimate_names);

  run;


  data keep_num_brr_replic(keep = num_brr_replic);
    set &data_brr_estimates_brecs end=lastbrr;

    num_brr_replic = _n_;
    if (lastbrr = 1) then output keep_num_brr_replic;

  run;


  data renamed_orig_estimates(keep = num_brr_replic origest1 - origest&k_est_names);
    if _n_=1 then set keep_num_brr_replic;
    set &data_orig_estimates_1rec(keep = &param_estimate_names);

    array estnames(&k_est_names) &param_estimate_names;
    array origest(&k_est_names)  origest1 - origest&k_est_names;

    do j = 1 to &k_est_names;
      origest(j) = estnames(j);
    end;

  run;


  data presumdata(drop = k);
    if _n_=1 then set renamed_orig_estimates;
    set &data_brr_estimates_brecs(keep = &param_estimate_names);

    array origest(&k_est_names) origest1 - origest&k_est_names;
    array brrest(&k_est_names)  &param_estimate_names;
    array presum(&k_est_names)  presum1 - presum&k_est_names;

    do k = 1 to &k_est_names;
      presum(k) =  ( (brrest(k) - origest(k))**2 )/(num_brr_replic * (1 - &set_f_method_fay)**2);
    end;

  run;


  proc means data=presumdata(keep=presum1-presum&k_est_names) sum noprint;
    var presum1-presum&k_est_names;
    output out=brrvar sum=&param_estimate_names;  /*** Sum is BRR variance ***/
  run;


  data brrvar;
    set brrvar;
    sort_order_1est_2brrvar = 2;
  run;


  data orig_est_1rec_addsortvar(keep = sort_order_1est_2brrvar &param_estimate_names);
    set &data_orig_estimates_1rec;
    sort_order_1est_2brrvar = 1;
  run;


  data parambrrvar;
    /*** Set brrvar first since variable names are in controlled order. ***/
    /*** Records are then sorted before transpose.                      ***/
    set brrvar(keep = sort_order_1est_2brrvar &param_estimate_names)
        orig_est_1rec_addsortvar
        ;
  run;


  proc sort data=parambrrvar;
    by sort_order_1est_2brrvar;
  run;


  data parambrrvar;
    set parambrrvar(drop = sort_order_1est_2brrvar);
  run;


  proc transpose data=parambrrvar out=transposed_data;
  run;


  data brr_pvalue_ci_out(drop = brr_variance sq_estimate_div_brr_se alpha one_minus_half_alpha t_quant);
    set transposed_data(rename=(col1=estimate col2=brr_variance));

    brr_se = sqrt(brr_variance);

    if (brr_se > 0) then do;
      sq_estimate_div_brr_se = (estimate/brr_se)**2;
      pvalue_param0_probf_df_1_&set_degrees_freedom = 1 - cdf('f', sq_estimate_div_brr_se, 1, &set_degrees_freedom);
    end;

    alpha = 1 - &set_confidence_level;
    one_minus_half_alpha = 1 - alpha/2;

    t_quant = quantile('t', one_minus_half_alpha, &set_degrees_freedom);

    lowerlim_t_df_&set_degrees_freedom = estimate - t_quant * brr_se;
    upperlim_t_df_&set_degrees_freedom = estimate + t_quant * brr_se;

  run;


  %if %upcase(&print)^=N %then %do;

    proc print data=brr_pvalue_ci_out;
      id _name_;
      title%eval(&titles+1) "Listing of Output Dataset: brr_pvalue_ci_out, Macro Input Includes: set_confidence_level &set_confidence_level, set_degrees_freedom &set_degrees_freedom";
      format pvalue_param0_probf_df_1_&set_degrees_freedom pvalue8.5;
    run;
    title%eval(&titles+1);

  %end;
  

%mend brr_pvalue_ci;  /*** End macro ***/
