The PSMATCH Procedure

(View the complete code for this example.)

This example illustrates how you can create observation weights that are appropriate for estimating the average treatment effect (ATE) in a subsequent outcome analysis (the outcome analysis itself is not shown here).

The data for this example are observations on patients in a nonrandomized clinical trial. The trial and the Drugs data set that contains the patient information are described in the section Getting Started: PSMATCH Procedure.

The following statements specify a logistic regression model for obtaining propensity scores, compute observation weights from the propensity scores, request statistics and plots for balance assessment, and save the weights in an output data set:

ods graphics on;
proc psmatch data=drugs region=allobs;
   class Drug Gender;
   psmodel Drug(Treated='Drug_X')= Gender Age BMI;
   psweight weight=atewgt nlargestwgt=6;
   assess lps var=(Gender Age BMI)
          / varinfo plots=(barchart boxplot(display=(lps BMI)) wgtcloud);
   id BMI;
   output out(obs=all)=OutEx1 weight=_ATEWgt_;
run;

The PSMODEL statement specifies the logistic regression model that creates the propensity score for each observation, which is the probability that the patient receives Drug_X. The CLASS statement specifies the classification variables in the model. The Drug variable is the binary treatment indicator variable, and TREATED='Drug_X' identifies Drug_X as the treated group. The Gender, Age, and BMI variables are included in the model because they are believed to be related to the assignment. The REGION=ALLOBS option specifies that the support region contains all observations. Weights are computed for all observations, regardless of the REGION= option.

The "Data Information" table in Output 101.1.1 displays the numbers of observations in the treated and control groups, the lower and upper limits of the propensity scores for observations in the support region, and the numbers of observations in the treated and control groups that fall within the support region. Because REGION=ALLOBS is specified, the lower and upper limits for of the propensity scores for observations in the support region are the minimum and maximum of the propensity scores for all observations. Consequently, all 373 observations in the control group fall within the support region, and all 133 observations in the treated group fall within the support region.

The "Propensity Score Information" table in Output 101.1.2 displays summary statistics by treatment group for all observations (labeled "All"), for observations in the support region (labeled "Region"), and for weighted observations in the support region (labeled "Weighted"). Because the support region consists of all observations, the first two rows in the table are identical. The WEIGHT=ATEWGT option in the PSWEIGHT statement displays summary statistics for ATE weighted observations.

The ASSESS statement produces tables and plots, shown in Output 101.1.3 through Output 101.1.5 and in Output 101.1.7 through Output 101.1.10, that summarize differences in the distributions of specified variables between treated and control groups. As requested by the LPS and VAR= options, these variables are the logit of the propensity score and the data variables Gender, Age, and BMI. Differences are summarized for all observations and for observations in the support region. Again, these two sets of differences are identical because REGION=ALLOBS is specified. The WEIGHT=ATEWGT option in the PSWEIGHT statement requests that differences in the variables also be summarized for the weighted observations. By comparing the differences for weighted observations to the differences for observations in the support region, you can assess how well weighting improves the balance for each variable.

The VARINFO option requests the "Variable Information" table, shown in Output 101.1.3, which displays variable summary statistics and differences between the treated and control groups for all observations (labeled "All"), for observations in the support region (labeled "Region"), and for weighted observations (labeled "Weighted"). For the binary classification variable (Gender), the difference is in the proportion of the first ordered level (Female).

The statistics in Output 101.1.3 are identical for all observations and for observations in the support region because REGION=ALLOBS is specified.

As indicated in the column labeled Weight, the total weight of the treated units is 460.45 and the total weight of the control units is 489.59, which are close to 486, the total number of units. The weights are ATE weights because WEIGHT=ATEWGT is specified in the PSWEIGHT statement. For information about ATE weights, see the section Inverse Probability of Treatment Weighting.

Note that in comparison to the unweighted means, the weighted means for the control units are closer in absolute value to the corresponding weighted means for the treated units.

The "Standardized Mean Differences" table, shown in Output 101.1.4, displays standardized mean differences in the variables between the treated and control groups, based on all observations, on observations in the support region, and on weighted observations.

The standardized mean differences based on weighted observations are significantly reduced; the largest of these differences is 0.12637 in absolute value, which is less than the upper limit of 0.25 that is recommended by Rubin (2001, p. 174) and Stuart (2010, p. 11). The treated-to-control variance ratios between the two groups are within the recommended range of 0.5 to 2. The percentage of reduction in variable mean difference is 0 for observations in the support region because REGION=ALLOBS is specified.

The PSMATCH procedure displays a standardized mean differences plot, shown in Output 101.1.5, for the variables that are specified in the ASSESS statement.

The "Standardized Mean Differences Plot" displays the differences that are shown in the "Standardized Mean Differences" table in Output 101.1.4. All differences for the weighted observations are within the recommended limits of –0.25 and 0.25, which are indicated by the shaded area.

The NLARGESTWGT=6 option displays the "Observations with Largest Weights" table, shown in Output 101.1.6, which lists the observations that have the six largest weights in the treated and control groups.

In the table, the scaled weights (which are the weights divided by their expected weights) are also displayed for ease of comparison. For more information about the expected weights in the treated and control group, see the section Propensity Score Weighting.

The PLOTS=WGTCLOUD option displays a cloud plot for the stabilized weights, which is shown in Output 101.1.7. This plot is called a cloud plot because the points are jittered in the vertical direction in order to avoid overplotting.

By default, the plot displays reference lines that represent 10 times the expected ATE weights in the treated and control groups. For information about these average weights, see the section Inverse Probability of Treatment Weighting.

The PLOTS=BARCHART option displays a bar chart for each classification variable that is specified in the ASSESS statement. As shown in Output 101.1.8, the bar chart shows the distributions of Gender based on all observations, on observations in the support region, and on weighted observations. By default, the bar chart displays the proportions of levels of Gender. Weighting the observations makes a slight improvement in the balance between males and females.

The PLOTS=BOXPLOT(DISPLAY=(LPS BMI)) option displays box plots for LPS and BMI, as shown in Output 101.1.9 and Output 101.1.10, respectively. These plots compare the distributions of the variables for the treated and control groups. Weighting the observations makes a good improvement in the balance between males and females.

Because there is good balance in the weighted distributions of the variables Gender, Age, and BMI, the observations and their weights can be saved in an output data set for use in a subsequent outcome analysis.

In situations where you are not satisfied with the variable balance, you can do one or more of the following to improve the balance: you can select another set of variables to fit the propensity score model, you can modify the specification of the propensity score model by using nonlinear terms for the continuous variables or by adding interactions (Rosenbaum and Rubin 1984), or you can choose another propensity score method (such as matching).

The OUT(OBS=ALL)=OutEx1 option in the OUTPUT statement creates an output data set, OutEx1, that contains all available observations. The following statements list the first 10 observations in OutEx1, as shown in Output 101.1.11.

proc print data=OutEx1(obs=10);
   var PatientID Drug Gender Age BMI _ps_ _AteWgt_;
run;

By default, the output data set includes the variable _PS_, which provides the propensity score. The weight for each treated unit is computed as 1 / p and the weight for each control unit is computed as 1 / (1 – p), where p is the propensity score.

After the responses for the trial are observed, they can be added to the data set OutEx1 as the starting point for an outcome analysis. Assuming that no other confounding variables are associated with both the response variable and the treatment group indicator Drug, you can estimate the ATE by performing a weighted version of the outcome analysis that you would have used to estimate the treatment effect if the original data set had resulted from a randomized trial.