/* This supplemental file contains SAS code used to generate estimates and tables corresponding to the paper: 
 “Assessing Alternative Precision Measures when Adjusting for Conditional Bias at the Subnational Level through 
 Calibration Weighting.” Due to disclosure concerns, some portions of the code were altered or are limited as 
 discussed in the program comments. */

/*** PROGRAM 1: RECALIBRATE NCVS WEIGHTS AND CALCULATE VICTIMIZATION RATES BASED ON TWO DIFFERENT METHODS ***/

/*** Note: Due to disclosure concerns, the control totals for recalibration cannot be released. However, this SAS program demonstrates how weights and estimates are calculated based on the conventional and corrected methods described in the paper. The control totals in this program are for demonstration purposes only and are NOT the true control totals used in the paper. These programs demonstrate the calculations for a single outcome (the overall violent victimization rate). The other personal victimization rates and household victimization rates are calculated similarly.  ***/

/***Calculate Conventionally Implemented Estimates ***/

/* Perform the weight calibration step */
proc wtadjust data=NCVSperson design=wr adjust=post maxiter=20;
        weight wt1; /* NCVS Weight - Uncalibrated */
        nest _one_ tract / missunit;
        idvar sudaanid; /* Unique Person Identifier */
        lowerbd 0;
        upperbd 1e10;
        postwgt 5000000  /* Total-Intercept */
                2700000 2300000 /* Gender */
                250000 300000 400000 550000 750000 1250000 1500000 /* Age */ 
                2500000 1000000 750000 400000 350000 /* Race */
                750000 800000 850000 1000000 900000 700000 /* FPL */
                3500000 1500000 /* Tenure */
                1000000 2500000 900000 600000 /* Education */
                3000000 600000 400000 1000000 /* Marital Status */
                2750000 175000 500000 /* Employment Status */
                ;
        class gender age raceeth fpl tenure educ marstat empwk / nofreq;
        model _one_=gender age raceeth fpl tenure educ marstat empwk;
        output / predicted=all filename=perwgt filetype=sas replace;
run;

/* Merge the recalibrated weight (WTFINAL) with the person-level analysis file */
proc sort data=NCVSperson; by sudaanid; run;
proc sort data=perwgt out=perwgt2(keep=sudaanid wtfinal); by sudaanid; run;

data person_with_wtfinal;
merge NCVSperson perwgt2;
by sudaanid;
run;

/* Calculate the rate of violent victimization using the recalibrated weight.*/
proc vargen data=person_with_wtfinal design=wr;
  nest _one_ tract/missunit;
  weight wtfinal; /* Use the recalibrated weight */
  xmean violent_rate: TOTVIOLENT*1000;
  output estim seestim/filename=violent_rate_conventional filetype=sas replace; /* Output the violent victimization rate and standard error */
run;


/***Calculate Corrected Estimates ***/

/* Calculate the violent victimization rate at the time of calibration */
proc wtadjust data=NCVSperson design=wr adjust=post maxiter=20;
        weight wt1; /* NCVS Weight - Uncalibrated */
        nest _one_ tract / missunit;
        idvar sudaanid; /* Unique Person Identifier */
        lowerbd 0;
        upperbd 1e10;
        postwgt 5000000  /* Total-Intercept */
                2700000 2300000 /* Gender */
                250000 300000 400000 550000 750000 1250000 1500000 /* Age */ 
                2500000 1000000 750000 400000 350000 /* Race */
                750000 800000 850000 1000000 900000 700000 /* FPL */
                3500000 1500000 /* Tenure */
                1000000 2500000 900000 600000 /* Education */
                3000000 600000 400000 1000000 /* Marital Status */
                2750000 175000 500000 /* Employment Status */
                ;
        class gender age raceeth fpl tenure educ marstat empwk / nofreq;
        model _one_=gender age raceeth fpl tenure educ marstat empwk;
        output / predicted=all filename=perwgt filetype=sas replace; /* Output the final calibrated weight */
        var TOTVIOLENT; 
        tables _one_;
        output mean se_mean / filename=violent_rate_corrected filetype=sas MEANFMT=F10.9 SE_MEANFMT=F10.9 replace; /* Output the violent victimization rate and standard error */
run;


/*** PROGRAM 2: COMPARE ESTIMATES ACROSS METHODS AND CREATE TABLES FOR THE PAPER ***/

/*** This SAS program compares the calibrated and recalibrated victimization rates and variance estimates between the conventional and corrected methods described in the paper. The input dataset contains one row per victimization estimate type, with variables for each method’s estimated victimization rate, standard error, and RSE. Note: In the code, variables with _VAR were calculated using the conventional method, variables with _WT were calculated using the corrected method, and variables with _uncal are the nationally-weighted estimates (i.e., they were calculated based on the original NCVS weights). ***/
data compare_methods;
 set estimates;

/*Calculate the absolute and percent change in victimization rates resulting from recalibration*/
 PTESTDIFF=abs(VR3_wt-VR3_uncal);
 ESTCHANGE=(abs(VR3_var-VR3_uncal)/VR3_uncal)*100; 

/*Calculate percent RSE variables*/
 PCTRSEVR3_uncal=RSEVR3_uncal*100;
 PCTRSEVR3_var=RSEVR3_var*100;
 PCTRSEVR3_wt=RSEVR3_wt*100;

/*Classify estimates based on RSE of nationally-weighted estimates */
 if RSEVR3_uncal LE .2 then RSEORIGCAT=1;
 else if RSEVR3_uncal LE .3 then RSEORIGCAT=2;
 else if RSEVR3_uncal LE .4 then RSEORIGCAT=3;
 else if RSEVR3_uncal LE .5 then RSEORIGCAT=4;
 else RSEORIGCAT=5;

/* Calculate absolute change in RSEs due to recalibration for each method (_TRAD = conventional, _ALT=corrected) */
 RSECHANGE_TRAD=PCTRSEVR3_var-PCTRSEVR3_uncal;
 RSECHANGE_ALT=PCTRSEVR3_wt-PCTRSEVR3_uncal;

/* Calculate the absolute difference in RSEs between the conventional and corrected methods */
 RSEDIFF=PCTRSEVR3_var-PCTRSEVR3_wt;

/* Classify estimates based on the percent change due to recalibration */
 if ESTCHANGE LE 2.5 then ESTCHANGECAT=1;
 else if ESTCHANGE LE 5 then ESTCHANGECAT=2;
 else if ESTCHANGE LE 7.5 then ESTCHANGECAT=3;
 else if ESTCHANGE LE 10 then ESTCHANGECAT=4;
 else ESTCHANGECAT=5;

run;

/* Calculate the distribution of absolute changes in victimization rates resulting from recalibration by subnational area type (Table 1) */
proc means data= compare_methods nway; 
 class area_type; /* State or MSA */
 var PTESTDIFF; 
 output out=distn_byarea(drop=_TYPE_ _FREQ_) min=min_dif p10=p10_dif p25=p25_dif median=med_dif p75=p75_dif p90=p90_dif max=max_dif;
run;

/* Calculate the distribution of differences in percent RSE between conventionally implemented and corrected variance methods by subnational area type (Table 2) */
proc means data= compare_methods nway; 
 class area_type; /* State or MSA */
 var RSEDIFF; 
 output out=distn_byarea(drop=_TYPE_ _FREQ_) min=min_RSE p10=p10_RSE p25=p25_RSE median=med_RSE p75=p75_RSE p90=p90_RSE max=max_RSE;
run;

/*Calculate the mean percent change in estimates due to subnational recalibration by nationally weighted RSE and subnational area type (Table 3) */
proc means data= compare_methods mean; 
 where area_type=1; /* States */
 class RSEORIGCAT; 
 var ESTCHANGE; 
 output out=BYRSE_states(drop=_TYPE_ _FREQ_) n=num_ests_states  
        mean=estchange_states; 
run;

proc means data= compare_methods mean; 
 where area_type=2; /* MSAs */
 class RSEORIGCAT; 
 var ESTCHANGE; 
 output out=BYRSE_MSAs(drop=_TYPE_ _FREQ_) n=num_ests_MSAs 
        mean=estchange_MSAs; 
run;

data BYRSE;
 merge BYRSE_states 
       BYRSE_MSAs;
 by RSEORIGCAT;
run;

/*Calculate the mean change in percent RSE for conventionally implemented and corrected estimation methods by percent change in estimate and subnational area type (Table 4) */

proc means data=compare_methods mean; 
 where area_type=1; /* States */
 class ESTCHANGECAT; 
 var RSECHANGE_TRAD RSECHANGE_ALT; 
 output out=BYEST_states(drop=_TYPE_ _FREQ_) n=num_ests_states 
        mean(RSECHANGE_TRAD)=TRAD_states mean(RSECHANGE_ALT)=ALT_states;
run;

proc means data=compare_methods mean; 
 where area_type=2; /* MSAs */
 class ESTCHANGECAT; 
 var RSECHANGE_TRAD RSECHANGE_ALT; 
 output out=BYEST_MSAs(drop=_TYPE_ _FREQ_) n=num_ests_MSAs 
        mean(RSECHANGE_TRAD)=TRAD_MSAs mean(RSECHANGE_ALT)=ALT_MSAs;
run;

data BYEST;
 merge BYEST_states
       BYEST_MSAs(in=b);
 by ESTCHANGECAT;
run;