Chapter 4

Doubly Robust Estimation of Treatment Effects

Michele Jonsson Funk
Daniel Westreich
Chris Weisen
Marie Davidian

Abstract

4.1 Introduction

4.2 Implementation with the DR Macro

4.3 Sample analysis

4.4 Summary

4.5 Conclusion

References

 

Abstract

Estimation of the effect of a treatment or exposure with a causal interpretation from studies where exposure is not randomized may be biased if confounding is not taken into appropriate account. Adjustment for confounding is often carried out through regression modeling of the relationships among treatment, confounders, and outcome. Doubly robust (DR) estimation produces a consistent effect estimator as long as one of two component regression models is correctly specified and assuming that there are no unmeasured confounders, giving the analyst two chances to correctly specify at least one of the regression models. In this chapter, we provide a brief introduction to DR estimators; present sample code using a SAS macro; illustrate the use of the macro with results from analyses of simulated data; and discuss issues including interpretation of estimates, assumptions, and limitations of this approach.

 

4.1 Introduction

Correct specification of the regression model is one of the most fundamental assumptions in statistical analysis. In an observational data analysis, it is common to estimate the causal effect of treatment using a regression model for the relationship between outcome, treatment, and confounders. Even when all relevant confounders have been measured, an unbiased estimator for the causal treatment effect can be obtained only if the model itself reflects the true relationship among treatment, confounders, and outcome. This is the case for a typical analysis in which the outcome is modeled as a function of exposure and covariates as well as propensity score-based methods, in which the exposure is modeled as a function of covariates. Outside of simulation studies, we can never know whether or not the model we have constructed accurately depicts those relationships. Thus, correct specification of the regression model is an unverifiable assumption. The DR estimator does not eliminate the need for such an assumption but does give the analyst two chances to satisfy it.

 

4.1.1 Conceptual Overview

Doubly robust (DR) estimation builds on the propensity score approach of Rosenbaum and Rubin (1983) and the inverse probability of weighting (IPW) approach of Robins and colleagues (Robins, 1998; Robins, 1999a; Robins, 1999b; Robins et al., 2000). DR estimation combines inverse probability weighting by a propensity score with regression modeling of the relationship between covariates and outcome for each treatment. It combines it in such a way that, as long as either the propensity score model or the outcome regression models are correctly specified, the effect of the exposure on the outcome will be correctly estimated, assuming that there are no unmeasured confounders (Robins et al., 1994; Robins, 2000; van der Laan and Robins, 2003; Bang and Robins, 2005). Specifically, one builds and fits a (binary) regression model for the probability that a particular patient received a given treatment as a function of that individual’s covariates (the propensity score). Maximum likelihood regression is conducted separately within the exposed and unexposed populations to predict the mean response (outcome) as a function of confounders and risk factors. (These two sets of models are visually represented in Figure 4.1.) Finally, each individual observation is given a weight equal to the inverse of the probability of the treatment he/she received based on baseline covariates (as in IPW analysis) to create two pseudopopulations of subjects that represent the expected response in the entire population under those two treatment conditions. Results from simulations confirm that the estimator is consistent when an important confounder is omitted entirely from one of the two models (Lunceford and Davidian, 2004; Bang and Robins, 2005) and, in more realistic scenarios, when one of the component models is misspecified by categorizing a continuous variable when the true relation with the outcome is a function of the continuous form (Jonsson Funk and Westreich, 2008).

Figure 4.1 Component Models of the DR Estimator

The DR estimator is an alternative to the usual approach of estimating the causal treatment effect based on a regression model for the relationships among the outcome and covariates and treatment (or using standard propensity scoring adjustment methods). If the outcome regression model is correctly specified, then the estimator for the causal effect will be at least as precise asymptotically as the DR estimator. However, if the outcome regression is misspecified, the resulting causal effect estimator need not be consistent for the true casual effect. The DR estimator is consistent if either the propensity score or treatment-specific outcome regression models are correct; thus, one trades a possible loss of precision in using the DR estimator for this additional protection.

The DR effect estimates have a marginal, rather than a conditional (on covariates), interpretation and are directly comparable to the effect estimates that one would obtain from a randomized trial in which a population is randomly assigned to receive treatment. Because the estimates from a standard outcome regression model have a conditional interpretation, the two estimates might not agree simply because they are averaging the effect in two different target populations. In particular, this could arise in the presence of effect measure modification (Kurth et al., 2006) or due to the non-collapsibility of the effect estimate (StÜrmer et al., 2006; Petersen et al., 2006).

 

4.1.2 Statistical Expression and Assumptions

We use the following notation: Y is the observed response or outcome, Z is a binary treatment (exposure) variable taking values 0 or 1, and X represents a vector of baseline covariates. Y1 and Y0 are the counterfactual responses under treatment and no treatment, respectively (Hernán, 2004). All of these variables are further subscripted by i for subjects i=1,…, n. In this example, the causal effect of interest is the difference in means if everyone in the population received treatment versus everyone not receiving treatment, or D=E(Y1)-E(Y0). In the following equation, e(X,β) is a postulated model for the true propensity score (from logistic regression), and m0(X,α0) and m1(X,α1) are postulated regression models for the true relationship between the vector of covariates (confounders plus other prognostic factors) and the outcome within each stratum of treatment. With these definitions, the estimator of the causal effect is:

The standard error for the DR estimator is estimated using the delta method (Casella and Berger, 2002). The sampling variance for the DR estimator is calculated as:

$$S{E}_{{\widehat{\mathrm{\Delta }}}_{DR}}=\sqrt{n^{-2}\sum _{i=1}^n{\widehat{I}}_i^2}$$

where Î_i is defined as follows:

$${\widehat{I}}_i=\left[\frac{Z_i{Y}_i}{e(X_i,\widehat{\beta })}-\frac{\left\{Z_i-e(X_i,\widehat{\beta })\right\}}{e(X_i,\widehat{\beta })}{m}_1(X_i,{\widehat{\alpha }}_1)\right]-\left[\frac{(1-Z_i)Y_i}{1-e(X_i,\widehat{\beta })}+\frac{\left\{Z_i-e(X_i,\widehat{\beta })\right\}}{1-e(X_i,\widehat{\beta })}{m}_0(X_i,{\widehat{\alpha }}_0)\right]-{\widehat{\mathrm{\Delta }}}_{DR}$$

 

The DR estimator is consistent for the true value of D if the following assumptions are satisfied: 1) no unmeasured confounding (also called the assumption of exchangeability); 2) correct specification of at least one of the two component models; 3) positivity (the true propensity score is bounded away from 0 and 1 so that there is always a positive probability of receiving both treatments under any combination of covariates); and 4) stable unit treatment value (SUTVA) (Rubin, 1980), comprising consistency and no interference (Cox, 1958). The delta method standard error is appropriate when the sample size is sufficiently large and when both the propensity score model and the outcome regression models have been correctly specified. For reporting estimates from an analysis of data where the true propensity score and outcome regression models are not known and may have been misspecified, standard errors should be obtained using the bootstrap.

 

4.2 Implementation with the DR Macro

While several estimators have been found to have the doubly robust property (Robins et al., 2007), we will describe the augmented IPW estimator identified by Robins and colleagues (1994) that was subsequently recognized to be doubly robust (Scharfstein et al., 1999). This DR estimator has been implemented in a macro developed at the University of North Carolina at Chapel Hill for use with Base SAS (validated with SAS 9.1.3). The DR macro can be downloaded from the Resources section of http://harryguess.unc.edu/sas.htm along with the sample data. This chapter reflects the features of the macro as of version 1.0. Please review the readme file (ReadMe.pdf) for important information regarding setup and installation as well as updated details of features prior to use.

 

4.2.1 Getting Started

The DR macro runs two sets of models: one for the probability of receiving a dichotomous treatment or exposure (the propensity score) and another to predict either the probability of the outcome (for a dichotomous outcome) or its mean value (for a continuous outcome) within strata of the exposure (treatment-specific outcome regression models). We describe the general syntax for each model statement, the use of optional commands, and the resulting output in the following sections.

 

4.2.2 Specifying the Weight Model

The general weight (or propensity score) model is specified in the first MODEL statement using the form:

wtmodel exposure = <covariates> / method=dr dist=bin <other options>

We model the main exposure or treatment on the left side of the equal sign as a function of the covariates on the right side. Method=dr indicates that the DR estimation method should be used, while dist=bin indicates that this is a binary exposure. Therefore, logistic regression will be used to model the relationship between the covariates and the exposure.

The weight model should be specified with the same care and rigor that you would use in the specification of any other IPW or propensity score model based on substantive knowledge. Brookhart and colleagues (2006) have found that the propensity score (exposure) model should include all confounders as well as those covariates that are risk factors for the outcome. Including covariates unrelated to the outcome can result in increased variance of the ultimate effect estimator of interest without reducing bias, although Greenland (2008) has recently suggested that the net benefit of variable reduction is limited and that it is preferable to report results from a full model.

In order to assess the positivity assumption, the analyst can invoke the SHOWCURVES option in the WTMODEL statement to produce graphs of the propensity score distribution in the two treatment groups. To limit the range of data for a sensitivity analysis, the analyst might enter an optional common_support=number, where the number is between 0 and 1 (inclusive). This limits the region of analysis to those observations for which there is common support for counterfactual inference. A value of 1 indicates that the program should use the entire region of common support; a value of 0.8 indicates that the program should use only the middle 80% of the region of common support. If the lack of common support is limited to one tail of the distribution, specifying common_support=1 excludes only those observations in the affected tail. There is no option at this time that allows for defining an asymmetric region within the range of propensity score values for which there is some overlap (5th percentile through 80th percentile, for instance). If common_support=number is not specified, the macro uses all available data, including any regions of non-overlapping propensity score distributions. In addition, the region of common support appears as vertical lines in the plot generated by the SHOWCURVES option.

 

4.2.3 Specifying the Outcome Regression Models

The outcome regression models are specified in the second MODEL statement using the following form:

model outcome = <covariates> </options>

This models the main outcome of interest (continuous or dichotomous) as a function of the covariates within each exposure group. The main exposure should not be specified here a second time. Options should be specified after a slash (/). Dist=n indicates that the model form is a linear regression, with a normal error distribution. This is appropriate for a continuous outcome. If the outcome is dichotomous, dist=bin is appropriate and logistic regression is used.

We know of no studies at this time that have specifically addressed which variables to include in the outcome regression models in the context of DR estimation, but we suggest that the analyst develop this model with the same rigor that would be used in a single multivariable regression model (Harrell, 2001).

 

4.2.4 Output

The DR macro outputs several results nodes in the results pane (the vertically oriented pane, to the left of main window). Each of these nodes is described below.

+ Logistic: Weight Model

This is the logistic regression for the propensity score model (wtmodel) portion of DR estimation. The related SAS code from the macro is as follows, where &desc corresponds to the descending option and &wtm is of the form exposure = <covariates>:

Program 4.1 Excerpt of SAS Code from DR Macro v1.0 for Weight Model

proc logistic data=_usedata_ &desc;
  ods exclude modelinfo;
 ods exclude nobs;
  model &wtm;
 output out=_ps_ p=__ps;
run;

+ Print: Descriptive statistics for Exposure Probability

This node provides the mean, standard deviation, and minimum and maximum of the estimated propensity scores stratified by exposure group.

+ Print: Extreme Values by Exposure Group

This node provides the five highest and five lowest propensity score values in each exposure group.

+ Univariate: Propensity Score Curves Stratified by <Exposure>

This is the graph of the propensity score curves stratified by exposure group generated by the WTMODEL statement only if the SHOWCURVES option is specified. If the COMMON_SUPPORT option is also invoked, this graph has lines indicating the region of common support, as well. Program 4.2 shows the SAS code that produces the histograms, where &expvar represents the exposure variable and &cmspt is a flag to indicate whether or not the common support option was invoked:

Program 4.2 Excerpt of SAS Code from DR Macro v1.0 for Propensity Score Curves Stratified by Exposure Group

proc univariate noprint data=_shcrvs_;
  var __ps;
  class &expvar;
  histogram __ps/kernel
  %if &cmspt=1 %then %do;
        href=(&lo_spt, &hi_spt)
  %end;
  ;
run;

+ Print: Observation Information for Exposure=0

This node reports for the first exposure group the total number of observations read from the original data set, the number with no missing covariate values, and the number used in the outcome regression model after applying any optional common support criteria.

+ Logistic | GLM: Outcome Model for Exposure=0

This node reports the results of the outcome regression model within the first exposure group. If the distribution (DIST=) option of the MODEL statement is set to bin, then this node is labeled Logistic; if the distribution option is set to normal (dist=n), then this node is labeled GLM. Program 4.3 shows the related SAS code for a dichotomous outcome from the macro, where &desc corresponds to the descending option and &outmod is of the form outcome = <covariates>:

Program 4.3 Excerpt of SAS Code from DR Macro v1.0 for Outcome Regression Models within Exposure Groups for a Dichotomous Outcome

title "Outcome Model for &expvar=&_exp0_";
    proc logistic data=_wgts0_ &desc;* noprint;
      ods exclude modelinfo;
     ods exclude nobs;
     model &outmod;
     output out=_modres0_(keep=__m0) p=__m0;
    run;

title "Observation Information for &expvar=&_exp1_";
    proc logistic data=_wgts1_ &desc;* noprint;
      ods exclude modelinfo;
     ods exclude nobs;
     model &outmod;
     output out=_modres1_(keep=__m1) p=__m1;
    run;

Program 4.4 shows the related SAS code for a continuous, normally distributed outcome from the macro, where &outmod is of the form outcome = <covariates>:

Program 4.4 Excerpt of SAS Code from DR Macro v1.0 for Outcome Regression Models within Exposure Groups for a Continuous Outcome

title "Outcome Model for &expvar=&_exp0_";
   proc glm data=_wgts0_;
     model &outmod;
     output out=_modres0_(keep=__m0) p=__m0;
   run;

title "Outcome Model for &expvar=&_exp1_";
   proc glm data=_wgts1_;
     model &outmod;
     output out=_modres1_(keep=__m1) p=__m1;
   run;

+ Print: Observation Information for Exposure=1

This node reports for the second exposure group the total number of observations read from the original data set, the number with no missing covariate values, and the number used in the outcome regression model after applying any optional common support criteria.

+ Logistic | GLM: Outcome Model for Exposure=1

This node reports the results of the outcome regression model within the second exposure group. If the distribution (DIST=) option of the MODEL statement is set to bin, then this node is labeled Logistic; if the distribution option is set to normal (dist=n), then this node is labeled GLM. The related SAS code from the macro is shown in Programs 4.3 and 4.4.

+ Print: DR Estimate

This results node produces output in the following format:

Output Displayed in Print DR Estimate Node

Total Obs

is the number of observations in the specified data set.

Used Obs

is the number of observations actually used in the analysis. Usedobs is always less than or equal to totalobs; usedobs is less than totalobs when there are missing values for some individuals for some covariates and/or when the COMMON_SUPPORT option is used in the WTMODEL statement.

Modeled Mean for Exposure=0

is an estimate of the average response under the assumption that no one in the population receives treatment. This is the mean of dr0.

Modeled Mean for Exposure=1

is an estimate of the average response under the assumption that everyone in the population receives treatment. This is the mean of dr1.

In the case of a continuous outcome, this is the mean value for that continuous variable (such as blood pressure, cholesterol, weight). In the case of a dichotomous outcome, this is the average risk of the outcome. This is the expected mean response in subjects rather than the expected value for an average subject; these two values are the same in a linear model but not so for a logistic regression model, which could lead to discrepancies.

Deltadr

is the difference between the mean responses (dr1–dr0). In the case of a continuous outcome, this is the mean difference due to treatment or exposure. In the case of a dichotomous outcome, this is the difference in the average predicted probability of the outcome, comparing the response rate in the pseudopopulation if everyone had been unexposed (or untreated) to the response rate in the pseudopopulation had everyone been exposed (or treated).

Sedeltadr

is the delta method standard error associated with the measure deltadr and should be reported only when the weight and outcome regression models are known to be correctly specified.

Z

is the Z score based on the delta method standard error.

ProbZ

is the two-tailed probability that deltadr=0 based on the Z score. Like the estimated standard error and Z score, this p value should be reported only when the weight and outcome regression models are known to be correctly specified.

Program 4.5 shows the related SAS code.

Program 4.5 Excerpt of SAS Code from DR Macro v1.0 for Calculation of DR Estimate of the Mean Difference Due to Exposure and Estimated Standard Error

/* combine M0 and M1 */
       data _modres01_(keep=&expvar &resvar __ps __m0 __m1 __exp01 __res01);
         merge _ps_ _modres0_ _modres1_;
         %if ("&desc"="") %then %do;
           __ps=1-__ps;
              __m0=1-__m0;
              __m1=1-__m1;
         %end;
       run;

/* create DR0 and DR1 and their difference DR1_DR0 statistics */
       data _dr01_;
         set _modres01_;
         dr0=((1-__exp01)*__res01+(__exp01-__ps)*__m0)/(1-__ps);
         dr1=(__exp01*__res01-(__exp01-__ps)*__m1)/__ps;
         dr1_dr0=dr1-dr0;
    run;

/* obtain mean, variance and n of difference DR1_DR0 */
/* and the means of DR0 and DR1 */
       proc means noprint data=_dr01_ vardef=n;
         var dr1_dr0 dr0 dr1;
         output out=_mdr01_(drop=_type_) mean=deltadr dr0 dr1 var=i2 vdr0 vdr1 n=__n;
       run;

/* get the SE of the difference */
/* and the two variance components */
       data _mdr01_;
         merge _mdr01_;
         SEdeltadr=sqrt(i2/__n);
         vdr0=vdr0/__n;
         vdr1=vdr1/__n;
       run;

In the case of a dichotomous outcome, a separate table with parameter estimates is displayed in the output. The following parameters, with their standard errors, Z scores and p values, are displayed. Program 4.6 shows the relevant SAS code.

LogRiskRatio

The natural log of the risk ratio. Exponentiating this value returns the risk ratio.

LogOddsRatio

The natural log of the odds ratio. Exponentiating this value returns the odds ratio.

Program 4.6 Excerpt of SAS Code from DR Macro v1.0 for Calculation of DR Estimates of the Mean Difference, Log Relative Risk, Log Odds Ratio, and Estimated Standard Errors

/* if the distribution of the response is binary */
/* get the ratio and the standard error of the ratio */
    %if (&mdist=BINOMIAL or &mdist=B or &mdist=BIN) %then %do;

      /* obtain the 2 variances and the covariance */
      ods listing close;
      proc corr cov data=_dr01_ vardef=n;
       var dr0 dr1;
        ods output cov=_cov_;
      run;
      ods listing;
      data _cov_;
        set _cov_;
        if _n_=1 then do;
          v0=dr0;
           v01=dr1;
        end;
        else do;
          v1=dr1;
         keep v0 v1 v01;
         output;
        end;
        retain v0 v01;
      run;

/* risk ratio, odds ratio and standard errors */

/* combine the covariance information with the means information */
/* get the variance estimates of some functions of the means */
      data  _mdr01_;
        merge _mdr01_ _cov_;
/* derivatives */
/* log odds */
        alo=(1/(dr1*(1-dr1)));
        blo=(-1/(dr0*(1-dr0)));
/* mean difference */
        am=1;
        bm=-1;
/* log risk ratio */
        alrr=1/dr1;
        blrr=-1/dr0;
/* log odds ratio and log risk ratio */
         LogOddsRatio=log((dr1*(1-dr0))/((1-dr1)*dr0));
         LogRiskRatio=log(dr1/dr0);
/* special for DESCENDING */
         %if "&desc"="" %then %do;
           LogOddsRatio=-1*logoddsratio;
           LogRiskRatio=-1*logriskratio;
         %end;
/* standard error estimates */
    SELogOddsRatio=sqrt((alo*alo*v1+blo*blo*v0+2*alo*blo*v01)/__n);       
        SEMeanDifference=sqrt((am*am*v1+bm*bm*v0+2*am*bm*v01)/__n);           

SELogRiskRatio=sqrt((alrr*alrr*v1+blrr*blrr*v0+2*alrr*blrr*v01)/__n); 
         drop alo blo am bm alrr blrr v0 v1 v01;
      run;
  %end;

 

4.3 Sample analysis

 

4.3.1 Introduction to Sample Data

The examples presented here are based on a simulated observational cohort with 10,000 individuals. The main exposure is statin initiation at baseline (statin, p[statin=1]=0.51). The main outcomes are rmi1a (mean=-10.7, sd=4.7), which represents the change in lipid levels between the baseline and a follow-up visit, and mi1 (p[mi1=1]=0.19), which represents the occurrence of an acute myocardial infarction during the follow-up period for the cohort. The data set is structured as one record or observation per person with the individuals represented in rows. Because these data are simulated, we know the true causal effect of the exposure on the outcomes as well as the true relationships between the covariates and the exposure and between the covariates and the outcomes. In both cases, the true effect of statin use on the outcomes is null. While the true mean response for the continuous outcome in each treatment group is negative (representing a decrease in lipid levels at follow-up relative to the baseline), these are simulated data and the methods described here apply equally to an outcome with a mean positive value. The association between the exposure and the outcomes is confounded by seven variables, four continuous (Age, BMI, Chol, and Exer) and three dichotomous (Hs, Smk, and Hxcvd). In addition, there is one variable that is a risk factor only for the outcome (Female) and two variables that are risk factors only for the exposure (Black and Income).

To run these example analyses, download the simulated study data set from http://harryguess.unc.edu/sas.htm and create a libname for SampleData that points to the appropriate directory on your computer.

 

4.3.2 DR Analysis of a Continuous Outcome

This example represents an analysis of simulated data where the exposure of interest is statin use (statin) and the outcome of interest is a continuous cardiovascular disease score (rmi1a). The true effect of statin use on the outcome is null. Both the weight (propensity score) model and the regression model are specified correctly in this analysis (see Program 4.7):

Program 4.7 Call to DR Macro v1.0 for Analysis of rmi1a Outcome

title 'Continuous Example';
%dr(%str(options data=sampledata.study descending;
       wtmodel statin=hs smk hxcvd black bmi age income chol exer
         / method=dr dist=bin showcurves common_support=.99;
model rmi1a=hs female smk hxcvd bmi bmi2 age age2 chol exer /dist=n;) );

The first component of the output is the usual SAS output from a logistic regression model (Program 4.1). From this, we can confirm the total number of observations and that the probability modeled is statin=1. These results also allow the analyst to identify covariates that are strongly associated with exposure and assess the fit of the model.

Output Node from rmi1a Analysis: Logistic Weight Model

The next node in the results pane presents the mean, standard deviation, and minimum and maximum predicted probabilities (or propensity scores) stratified by exposure group. This allows the analyst to check the assumption of positivity. This assumption states that for each combination of characteristics, there must be a non-zero probability of being exposed and unexposed (Cole and Hernán, 2008). In the event that positivity is violated, consider whether this is a case of structural nonpositivity, in which it is not possible for individuals with a particular combination of characteristics to receive one of the exposures. Instances of structural nonpositivity suggest that these observations should not be included in the analysis. There may also be cases of random nonpositivity, particularly when some covariates are continuous. The regression models smooth over these instances of nonpositivity, but it is helpful to investigate the sensitivity of the findings to violation of the positivity assumption.

Output Node from rmi1a Analysis: Descriptive Statistics for Exposure Probability

The SHOWCURVES option produces a histogram that compares the distributions of the propensity score for the two levels of exposure with a nonparametric smoothed curve overlaid (see Figure 4.2, produced by Program 4.2). This allows the analyst to visually assess the degree to which there are unexposed individuals who can serve as counterfactuals for those who were exposed and vice versa.

Output Node from rmi1a Analysis: Propensity Score Curves Stratified by Exposure

Figure 4.2 Estimated Propensity Score Distributions Stratified by Exposure with Nonparametric Smoothed Curve

The COMMON_SUPPORT option directs the outcome regression model to trim off observations that lie at the extreme ends of the propensity score distribution in order to support sensitivity analyses.¹ Using common_support=0.99, the regression models are limited to those observations with a propensity score between the 0.5th percentile and the 99.5th percentile. The vertical dashed lines in Figure 4.2 indicate the boundaries for this portion of the data.

The following two nodes present the results of linear regression within the unexposed (statin=0) and exposed (statin=1, not shown) groups, respectively (Program 4.3). The number of observations used may be less than the total number of observations in the data set (the number of observations read). In this example, there are 9,858 total observations rather than 10,000 because of the use of the COMMON_SUPPORT option, of which 4,857 (out of the original 4,927) were in the unexposed group and, therefore, contribute to the outcome regression among the unexposed (shown here). In addition, observations without complete data for all covariates are excluded in a case-wise deletion manner.

Output Nodes from rmi1a Analysis: Observation Information for Exposure=0 and GLM Outcome Model for Exposure=0

Finally, under the node labeled Print: DR Estimate, we find the results of interest (Program 4.5). Again, the number of observations under Used Obs may be less than the total number of observations in the event of missing data for some observations for some covariates and/or use of the COMMON_SUPPORT option in the MACRO statement. Under Modeled Mean for Statin=0, we find the estimated mean response if all subjects in this cohort had been unexposed (statin=0) and, likewise under Modeled Mean for Statin=1, the estimated mean response if all subjects in this cohort had been exposed (statin=1). The doubly robust estimate (deltadr) of the average treatment effect (also known as the mean difference due to treatment) when we have specified both models correctly is a difference of 0.00 with a standard error of 0.02, a p value of 0.87, and a 95% confidence interval of -0.04 to 0.05 (calculated using the estimated standard error for the delta DR).

Output Node from rmi1a Analysis: DR Estimate

Results from a standard linear outcome regression model including treatment and all covariates can be obtained, for comparison, using the following SAS code:

Program 4.8 Standard Linear Regression Model

proc reg data=sampledata.study;
model rmi1a=statin hs female smk hxcvd bmi bmi2 age age2 chol exer;
run;

Results from this model are shown here. The mean difference due to statin use conditional on all other covariates in a single linear regression model is -0.01, similar to the DR estimate. The standard error of the DR estimate (0.024) is slightly greater than that from the single linear regression model (0.022), as we would expect. Recall that the DR estimate used a sample in which some observations at the extreme were trimmed, so these are not strictly comparable estimates.

Output from Program 4.8—Standard Linear Regression Model

 

4.3.3 DR Analysis of a Dichotomous Outcome

The second example represents an analysis of simulated data where the exposure of interest is statin use (statin) and the outcome of interest is a dichotomous variable indicating whether the subject experienced a myocardial infarction within the follow-up period (mi1). Both the weight (propensity score) model and the regression model are specified correctly:

Program 4.9 Call to DR Macro v1.0 for Analysis of mi1 Outcome

title 'Dichotomous Example';
%dr(%str(options data=sampledata.study descending;
       wtmodel statin=hs smk hxcvd black bmi age income chol exer
         / method=dr dist=bin showcurves;
model mi1=hs female smk hxcvd bmi bmi2 age age2 chol exer /dist=bin;
) );

As before, the first node in the results pane includes the complete output from the logistic regression model for the propensity score model. Because we specified the SHOWCURVES option, a graph of the two propensity score distributions is also produced. The results are identical because the same exposure and covariates were specified as in the previous example. There are two additional logistic regression nodes: the first for the model of the predicted response among the exposed and the second among the unexposed, respectively. The results of interest appear under the Print: DR estimates node (Program 4.6). These results include DR estimates of the risk (or probability) of the outcome had everyone in the population been untreated (dr0) or treated (dr1), the risk difference (delta DR), the log risk ratio, and the log odds ratio. In addition, the standard errors, Z score, and p value for each parameter are provided.

Output Node from mi1 Analysis: DR Estimate

Due to the noncollapsibility of the odds ratio, it is not particularly meaningful to compare the DR estimate of the odds ratio to that from a typical logistic regression model. Instead, we compare the DR estimate of the risk ratio to the risk ratio estimated using the GENMOD procedure (Spiegelman and Hertzmark, 2005). (Poisson regression is used in this case because of problems with convergence of binomial regression.)

Program 4.10 Analysis of mi1 Using Poisson Regression Model

title 'Estimate of RR using Poisson Regression';
proc genmod data=sampledata.study descending;
    class i;
    model mi1=statin hs female smk hxcvd bmi bmi2 age age2 chol exer
/dist=poisson link=log;
    repeated subject=i/type=ind;
    estimate 'log RR statin' statin 1/exp;
run;

The estimate of logrRR (0.12) is further from the true null value than the DR estimator (0.03), and the standard error of the logRR from Poisson regression is less efficient (0.04 vs. 0.03).

Output from Program 4.10–Poisson Regression Model

In the examples provided, the true forms of the propensity score and outcome regression models are used, but one can also explore the robustness of this method by intentionally misspecifying one of the two models using these sample data.

 

4.4 Summary

Given that we rarely know the true relations among exposures, outcomes, and confounders, the DR estimator is a potentially valuable tool for obtaining more robust effect estimates in observational studies of the effects (intended and unintended) of drugs, devices, and other interventions. The SAS macro described here makes this method relatively straightforward to apply to real-world analyses. But given the relative novelty of this estimator, we suggest that researchers take an especially careful approach to this type of analysis, with particular consideration of the method’s limitations, the SAS macro’s implementation in particular, and a variety of issues that remain under active debate. Given that this is a relatively new approach, we suggest using this method along with a more traditional bias adjustment approach as a sensitivity analysis.

 

4.4.1 Limitations

While the DR estimator provides the analyst with two chances to specify a regression model correctly, it still requires assumptions common to most regression models for causal effect estimation. While some of these can be verified, the most critical assumption of no unmeasured confounders is unverifiable. One of the two models must be correctly specified for the estimator to remain consistent, which means that all confounders must have been measured in order to be included in at least one of the models. In the event that neither model is correct, the estimator is no longer consistent and is not necessarily closer to the true effect than that from a single misspecified outcome regression model (Jonsson Funk and Westreich, 2008; Kang and Schafer, 2008).

Because this method reweights observations, the results can be sensitive to observations that are given a very large weight, as is the case with IPW methods more generally (Robins et al., 2000; Cole and Hernán, 2008). This can arise when an individual has a particular combination of characteristics that is almost always associated with one of the two exposure conditions. For instance, if there are 100 individuals over the age of 80 but only one of them is unexposed, that single unexposed individual will be reweighted in the final calculation of the treatment effect to count as 100 observations, making his or her outcome potentially overinfluential. In this situation, the analyst should first seek to understand which combinations of characteristics define these unusual individuals. Next, the analyst needs to make a decision based on the research question and substantive knowledge regarding whether these individuals are properly part of the population of interest. If so, then further investigations of effect measure modification would be appropriate in combination with sensitivity analyses to assess the degree to which the estimated treatment effect is sensitive to omission of individuals where there is limited positivity.

 

4.4.2 Practical Considerations

While there are other DR estimators and means of computing these estimators, the SAS macro for DR estimation described here has specific advantages and limitations. One advantage is that it provides estimates of the absolute risks, risk difference, and risk ratio in addition to the usual odds ratio for analyses of dichotomous outcomes. The macro also includes some built-in diagnostics to aid the analyst in evaluating the appropriateness of the propensity score (weight) model.

In terms of limitations, the current version (v1.0) is designed only to handle binary exposures and binary or continuous outcomes. In the event that there are more than two exposure groups, the analyst would need to conduct pairwise comparisons to make use of the macro. Estimated standard errors are provided for reference but appropriate estimates of the standard error and confidence limits should be obtained by bootstrapping for purposes of reporting results from original research where either the propensity score or outcome regression models may be misspecified. With respect to missing data, if a given observation has a missing value for any covariate, the individual’s data will not be utilized, as is the case with most SAS procedures. While methods such as multiple imputation can be used to more adequately address missing data, the analyst will need to do so outside of the macro. Although the code used by the analyst to run the macro was designed to look much like a typical SAS procedure to improve its usability, there are some SAS conventions that are not currently recognized. Specifically, variables for interaction terms and higher order terms must be created in a DATA step–not within the MODEL statements. The CLASS statement is also not recognized and, therefore, all categorical variables should be coded using indicator variables. Other common statements such as WHERE and BY are also not recognized. Development of the macro is ongoing, so the analyst should review the documentation provided with the current version on the Web site to be aware of any changes subsequent to this publication.

 

4.4.3 Areas of Onoing Investigation

Doubly robust estimation is a relatively new method. As such, many questions remain to be answered about its performance in applied analyses and its optimal use. Key questions under active discussion include detecting and handling effect modification, selecting covariates for the propensity score and outcome regression models, using diagnostics, and addressing violations of the assumption of positivity.

 

4.5 Conclusion

Doubly robust estimation methods–which provide the analyst with two chances to correctly specify the true relationships among exposures, covariates and outcomes–are potentially valuable tools for epidemiologic research. In this chapter, we have presented a basic introduction to the method, provided specific instruction on the use of DR estimation as implemented with a SAS macro, illustrated the use of the macro with two sample analyses of a simulated observational cohort, and discussed several issues that the analyst should take into consideration when using this approach.

 

References

Bang, H., and J. M. Robins. 2005. “Doubly robust estimation in missing data and causal inference models.” Biometrics 61(4): 962–973.

Brookhart, M. A., S. Schneeweiss, K. J. Rothman, R. J. Glynn, J. Avorn, and T. StÜrmer. 2006. “Variable selection for propensity score models.” American Journal of Epidemiology 163(12): 1149–1156.

Casella, G., and R. L. Berger. 2001. Statistical Inference. Australia: Duxbury Press.

Cole, S. R., and M. A. Hernán. 2008. “Constructing inverse probability weights for marginal structural models.” American Journal of Epidemiology 168(6): 656–64.

Cox, D. R. 1958. Planning of Experiments. New York: John Wiley & Sons, Inc.

Greenland, S. 2008. “Invited commentary: variable selection versus shrinkage in the control of multiple confounders.” American Journal of Epidemiology 167(5): 523–529; discussion 530–531.

Harrell, F. E. 2001. Regression Modeling Strategies: With Applications to Linear Models, Logistic Regression, and Survival Analysis. New York: Springer-Verlag.

Hernán, M. A. 2004. “A definition of causal effect for epidemiological research.” Journal of Epidemiology and Community Health 58(4): 265–71.

Jonsson Funk, M. L., and D. Westreich. 2008. “Doubly robust estimation under realistic conditions of model misspecification.” Pharmacoepidemiology and Drug Safety 17(S1): S106.

Kang, J. D.Y., and J. L. Schafer. 2007. “Demystifying double robustness: a comparison of alternative strategies for estimating a population mean from incomplete data (with discussion).” Statistical Science 22(4): 523–580.

Kurth, T., A. M. Walker, R. J. Glynn, K. A. Chan, J. M. Gaziano, K. Berger, and J. M. Robins. 2006. “Results of multivariable logistic regression, propensity matching, propensity adjustment, and propensity-based weighting under conditions of nonuniform effect.” American Journal of Epidemiology 163(3): 262–270.

Lunceford, J. K., and M. Davidian. 2004. “Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study.” Statistics in Medicine 23(19): 2937–2960.

Petersen, M. L., Y. Wang, M. J. van der Laan, and D. R. Bangsberg. 2006. “Assessing the effectiveness of antiretroviral adherence interventions: Using marginal structural models to replicate the findings of randomized controlled trials.” Journal of Acquired Immune Deficiency Syndromes 43: S96–S103.

Robins, J. M. 1997. “Marginal structural models.” In ASA Proceedings in the Section on Bayesian Statistical Science. Alexandria, VA: American Statistical Association, 1–10.

Robins, J. M. 1998. “Correction for non-compliance in equivalence trials.” Statistics in Medicine 17(3): 269–302; discussion 387–389.

Robins, J. M. 1999a. “Association, causation, and marginal structural models.” Synthese 121: 151–179.

Robins, J. M. 1999b. “Marginal structural models versus structural nested models as tools for causal inference.” In Statistical Models in Epidemiology: The Environment and Clinical Trials. Edited by M. E. Halloran and D. Berry. IMA 116: 95–134. New York: Springer-Verlag.

Robins, J. M. 1999c. “Robust estimation in sequentially ignorable missing data and causal inference models.” In ASA Proceedings in the Section on Bayesian Statistical Science 6–10.

Robins, J. M., A. Rotnitzky, and L. P. Zhao. 1994. “Estimation of regression coefficients when some of the regressors are not always observed.” Journal of the American Statistical Association 89: 846–866.

Robins, J. M., M. A. Hernán, and B. Brumback. 2000. “Marginal structural models and causal inference in epidemiology.” Epidemiology 11(5): 550–60.

Robins, J. M., M. Sued, Q. Lei-Gomez, and A. Rotnitzky. 2007. “Comment: Performance of double-robust estimators when ‘inverse probability’ weights are highly variable.” Statistical Science 22(4): 544–559.

Rosenbaum, P. R., and D. B. Rubin. 1983. “The central role of the propensity score in observational studies for causal effects.” Biometrika 70: 41–55.

Rubin, D. B. 1980. Comments on “Randomization analysis of experimental data in the Fisher randomization test.” Journal of the American Statistical Association 75: 591–593.

Scharfstein, D. O., A. Rotnitzky, and J. Robins. 1999. “Adjusting for nonignorable drop-out using semiparametric nonresponse models.” Journal of the American Statistical Association 94: C/R)1121–1146.

Spiegelman, D., and E. Hertzmark. 2005. “Easy SAS calculations for risk or prevalence ratios and differences.” American Journal of Epidemiology 162: 199–200.

Stefanski, L. A., and D. D. Boos. 2002. “The calculus of M-estimation.” The American Statistician 56: 29–38.

StÜrmer, T., K. J. Rothman, and R. J. Glynn. 2006. “Insights into different results from different causal contrasts in the presence of effect-measure modification.” Pharmacoepidemiology and Drug Safety 15(10): 698–709.

van der Laan, M. J., and J. M. Robins. 2003. Unified Methods for Censored and Longitudinal Data and Causality. New York: Springer-Verlag.

 

¹We have used it in this example only to demonstrate its application. There is no indication that it is needed in these data based on the good overlap across the full range of the propensity score distributions.