PSweight: An R Package for Propensity Score Weighting Analysis

Propensity score weighting is an important tool for comparative effectiveness research. Besides the inverse probability of treatment weights (IPW), recent development has introduced a general class of balancing weights, corresponding to alternative target populations and estimands. In particular, the overlap weights (OW) lead to optimal covariate balance and estimation efficiency, and a target population of scientific and policy interest. We develop the R package PSweight to provide a comprehensive design and analysis platform for causal inference based on propensity score weighting. PSweight supports (i) a variety of balancing weights, (ii) binary and multiple treatments, (iii) simple and augmented weighting estimators, (iv) nuisance-adjusted sandwich variances, and (v) ratio estimands. PSweight also provides diagnostic tables and graphs for covariate balance assessment. We demonstrate the functionality of the package using a data example from the National Child Development Survey (NCDS), where we evaluate the causal effect of educational attainment on income.

Tianhui Zhou (Department of Biostatistics and Bioinformatics, Duke University School of Medicine) , Guangyu Tong (Department of Biostatistics, Yale University School of Public Health) , Fan Li (Department of Statistical Science, Duke University) , Laine E. Thomas (Department of Biostatistics and Bioinformatics, Duke University School of Medicine) , Fan Li (Department of Biostatistics, Yale University School of Public Health)
2022-06-22

1 Introduction

Propensity score is one of the most widely used causal inference methods for observational studies (Rosenbaum and Rubin 1983). Propensity score methods include weighting, matching, stratification, regression, and mixed methods such as the augmented weighting estimators. The PSweight package provides a design and analysis pipeline for causal inference with propensity score weighting (Robins et al. 1994; Hirano et al. 2003; Lunceford and Davidian 2004; Li et al. 2018). There are a number of existing R packages on propensity score weighting (see Table 1). Comparing to those, PSweight offers three major advantages: it incorporates (i) visualization and diagnostic tools of checking covariate overlap and balance, (ii) a general class of balancing weights, including overlap weights and inverse probability of treatment weights, and (iii) multiple treatments. More importantly, PSweight comprises a wide range of functionalities, whereas each of the competing packages only supports a subset of these functionalities. As such, PSweight is currently the most comprehensive platform for causal inference with propensity score weighting, offering analysts a one-stop shop for the design and analysis. Table 1 summarizes the key functionalities of PSweight in comparison to related existing R packages. We elaborate the main features of PSweight below.

PSweight facilitates better practices in the design stage of observational studies, an aspect that has not been sufficiently emphasized in related packages. Specifically, we provide a design module that facilitates visualizing overlap (also known as the positivity assumption) and evaluating covariate balance without access to the final outcome (Austin and Stuart 2015). When there is limited overlap, PSweight allows for symmetric propensity score trimming (Crump et al. 2009; Yoshida et al. 2018) and optimal trimming (Crump et al. 2009; Yang et al. 2016) to improve the internal validity. We extend the class of balance metrics suggested in Austin and Stuart (2015) and Li et al. (2019) for binary treatments, and those in McCaffrey et al. (2013) and Li and Li (2019) for multiple treatments. In addition, the design module helps describe the weighted target population by providing the information required in the standard “Table 1” of a clinical article.

In addition to the standard inverse probability of treatment weights (IPW), PSweight implements the average treatment effect among the treated (Treated) weights, overlap weights (OW), matching weights (MW) and entropy weights (EW) for both binary (Li and Greene 2013; Li et al. 2018; Mao et al. 2018; Zhou et al. 2020) and multiple treatments (Yoshida et al. 2017; Li and Li 2019). All weights are members of the family of balancing weights (Li et al. 2018); the last three types of weights target at the subpopulation with improved overlap in the covariates between (or across) treatment groups, similar to the target population in randomized controlled trials (Thomas et al. 2020a,b). Among them, OW achieves optimal balance and estimation efficiency (Li et al. 2018, 2019). We also implement the augmented weighting estimators corresponding to each of the above weighting schemes (Mao et al. 2018). By default, PSweight employs parametric regression models to estimate propensity scores and potential outcomes. Nonetheless, it also allows for propensity scores to be estimated by external machine learning methods including generalized boosted regression models (McCaffrey et al. 2013) and super learner (Van der Laan et al. 2007), or importing any other propensity or outcome model estimates of interest.

Table 1: Comparisons of existing R packages that implement propensity score weighting with discrete treatments. Binary treatments and additive estimands are implemented in all packages, and therefore those two columns are omitted.
Multiple Balance IPW/ATT OW/other Ratio Augmented Nuisance-adj Optimal
treatments diagnostics weights weights estimands weighting variance trimming
PSweight \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\)
twang \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\times\) \(\times\) \(\times\) \(\times\) \(\times\)
CBPS \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\times\) \(\times\) \(\checkmark\) \(\checkmark\) \(\times\)
PSW \(\times\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\times\)
optweight \(\checkmark\) \(\times\) \(\checkmark\) \(\times\) \(\times\) \(\times\) \(\times\) \(\times\)
ATE \(\checkmark\) \(\checkmark\) \(\checkmark\) \(\times\) \(\times\) \(\times\) \(\checkmark\) \(\times\)
WeightIt \(\checkmark\) \(\times\) \(\checkmark\) \(\checkmark\) \(\times\) \(\times\) \(\times\) \(\times\)
causalweight \(\checkmark\) \(\times\) \(\checkmark\) \(\times\) \(\times\) \(\checkmark\) \(\times\) \(\times\)
sbw \(\times\) \(\checkmark\) \(\checkmark\) \(\times\) \(\times\) \(\times\) \(\times\) \(\times\)

\(\checkmark\) indicates that the functionality is currently implemented in the package; \(\times\) indicates otherwise.

References: twang (Version 1.6): Ridgeway et al. (2020); CBPS (Version 0.21): Fong et al. (2019); PSW (Version 1.1-3): Mao and Li (2018); optweight (Version 0.2.5): Greifer (2019); ATE (Version 0.2.0): Haris and Chan (2015); WeightIt (Version 0.10.2): Greifer (2020); causalweight (Version 0.2.1): Bodory and Huber (2020); sbw (Version 1.1.1): Zubizarreta and Li (2020).

To our knowledge, PSweight is the first R package to accommodate a variety of balancing weighting schemes with multiple treatments. Existing R packages such as twang (Ridgeway et al. 2020), CBPS (Fong et al. 2019), optweight (Greifer 2019) have also implemented weighting-based estimation with multiple treatments, but focus on IPW. The PSW R package (Mao and Li 2018) implements both OW and MW and allows for nuisance-adjusted variance estimation, but it is only for binary treatments. To better assist applied researchers to perform propensity score weighting analysis, this article provides a full introduction of the PSweight package. In what follows, we explain the methodological foundation of PSweight and outline the main functions and their arguments. We further illustrate the use of these functions with a data example that studies the causal effect of educational attainment on income, and finally conclude with a short discussion.

2 Overview of propensity score weighting

Before diving into the implementation details of PSweight, we briefly introduce the basics of the propensity score weighting framework.

Binary treatments

Assume we have an observational study with \(N\) units. Each unit \(i\) (\(i=1,2,\ldots,N\)) has a binary treatment indicator \(Z_{i}\) (\(Z_{i}=0\) for control and \(Z_{i}=1\) for treated), a vector of \(p\) covariates \(\mathbf{X}_{i}=(X_{1i},\cdots, X_{pi})\). For each unit \(i\), we assume a pair of potential outcomes \(\{Y_{i}(1),Y_{i}(0)\}\) mapped to the treatment and control status, of which only the one corresponding to the observed treatment is observed, denoted by \(Y_i=Z_{i}Y_{i}(1)+(1-Z_{i})Y_{i}(0)\); the other potential outcome is counterfactual.

Causal effects are contrasts of the potential outcomes of the same units in a target population, which usually is the population of a scientific interest (Thomas et al. 2020b). PSweight incorporates a general class of weighted average treatment effect (WATE) estimands. Specifically, assume the observed sample is drawn from a probability density \(f(\mathbf{x})\), and let \(g(\mathbf{x})\) denote the covariate density of the target population. The ratio \(h(\mathbf{x})\propto g(\mathbf{x})/f(\mathbf{x})\) is called the tilting function, which adjusts the distribution of the observed sample to represent the target population. Denote the conditional expectation of the potential outcome by \(m_z(\mathbf{x})=\mathop{\mathrm{\mathbb E}}[Y(z)|\mathbf{X}=\mathbf{x}]\) for \(z=0,1\). Then, we can represent the average treatment effect over the target population by a WATE estimand: \[\label{eq:estimand1} \tau^h=\mathop{\mathrm{\mathbb E}}_g[Y(1)-Y(0)]=\frac{\mathop{\mathrm{\mathbb E}}\{h(\mathbf{x})(m_1(\mathbf{x})-m_0(\mathbf{x}))\}}{\mathop{\mathrm{\mathbb E}}\{h(\mathbf{x})\}}. \tag{1}\] To estimate (1), PSweight maintains two standard assumptions: (1) unconfoundedness: \(\{Y(1),Y(0)\} \perp Z \mid \mathbf{X}\); (2) overlap: \(0<P(Z=1|\mathbf{X})<1\). The propensity score is the probability of a unit being assigned to the treatment group given the covariates (Rosenbaum and Rubin 1983): \(e(\mathbf{x})=P(Z=1|\mathbf{X}=\mathbf{x})\). While assumption (1) is generally untestable and critically depends on substantive knowledge, assumption (2) can be checked visually from data by comparing the distribution of propensity scores between treatment and control groups.

For a given tilting function \(h(\mathbf{x})\) (and correspondingly a WATE estimand \(\tau^h\)), we can define the balancing weights \((w_1,w_0)\) for the treated and control units: \(w_1(\mathbf{x}) \propto h(\mathbf{x})/{e(\mathbf{x})}\) and \(w_0(\mathbf{x}) \propto h(\mathbf{x})/\{1-e(\mathbf{x})\}\). These weights balance the covariate distributions between the treated and control groups towards the target population (Li et al. 2018). PSweight implements the following Hájek estimator for WATE: \[\label{eq:sampleWATE} \hat{\tau}^h=\hat{\mu}^h_1-\hat{\mu}^h_0=\frac{\sum_{i=1}^Nw_1(\mathbf{x}_i)Z_i Y_i}{\sum_{i=1}^Nw_1(\mathbf{x}_i)Z_i} - \frac{\sum_{i=1}^Nw_0(\mathbf{x}_i)(1-Z_i) Y_i}{\sum_{i=1}^Nw_0(\mathbf{x}_i)(1-Z_i)}, \tag{2}\] where the weights are calculated based on the propensity scores estimated from the data. Clearly, specification of \(h(\mathbf{x})\) defines the target population and estimands. PSweight primarily implements the following three types of balancing weights (see Table 2 for a summary):

IPW has been the dominant weighting method in the literature, but has a well-known shortcoming of being sensitive to extreme propensity scores, which induces bias and large variance in estimating treatment effects. OW addresses the conceptual and operational problems of IPW. Among all balancing weights, OW leads to the smallest asymptotic (and often finite-sample) variance of the weighting estimator (2). (Li et al. 2018, 2019). Recent simulations also show that OW provides more stable causal estimates under limited overlap (Yoshida et al. 2017; Mao et al. 2018; Yoshida et al. 2018; Li et al. 2019), and is more robust to misspecification of the propensity score model (Zhou et al. 2020).

PSweight implements two additional types of balancing weights: matching weights (MW) (Li and Greene 2013), and entropy weights (EW) (Zhou et al. 2020). Similar to OW, MW and EW focus on target populations with substantial overlap between treatment groups. Though having similar operating characteristics, MW and EW do not possess the same theoretical optimality as OW, and are less used in practice. Therefore, we will not separately describe MW and EW hereafter.

Table 2: Target populations, tilting functions, estimands and the corresponding balancing weights for binary treatments in PSweight.
Target population Tilting function \(h(\mathbf{x})\) Estimand Balancing weights \((w_1, w_0)\)
Combined 1 ATE \(\left(\frac{1}{e(\mathbf{x})}, \frac{1}{1-e(\mathbf{x})}\right)\)
Treated \(e(\mathbf{x})\) ATT \(\left(1, \frac{e(\mathbf{x})}{1-e(\mathbf{x})}\right)\)
Overlap \(e(\mathbf{x})(1-e(\mathbf{x}))\) ATO \((1-e(\mathbf{x}), e(\mathbf{x}))\)
Matching \(\xi_1(\mathbf{x})\) ATM \(\left( \frac{\xi_1(\mathbf{x})}{e(\mathbf{x})}, \frac{\xi_1(\mathbf{x})}{1-e(\mathbf{x})} \right)\)
Entropy \(\xi_2(\mathbf{x})\) ATEN \(\left( \frac{ \xi_2(\mathbf{x})}{e(\mathbf{x})}, \frac{ \xi_2(\mathbf{x})}{1-e(\mathbf{x})} \right)\)

Notes: \(\xi_1(\mathbf{x}) = \min \{e(\mathbf{x}),1-e(\mathbf{x})\}\) and \(\xi_2(\mathbf{x})=-\{e(\mathbf{x})\log(e(\mathbf{x}))+(1-e(\mathbf{x}))\log(1-e(\mathbf{x}))\}\).

In observational studies, propensity scores are generally unknown and need to be estimated. Therefore, propensity score analysis usually involves two steps: (1) estimating the propensity scores, and (2) estimating the causal effects based on the estimated propensity scores. In PSweight, the default model for estimating propensity scores with binary treatments is a logistic regression model. Spline or polynomial models can be easily incorporated by adding bs(), ns() or poly() terms into the model formula. PSweight also allows for importing propensity scores estimated from external routines, such as boosted models or super learner.

Goodness-of-fit of the propensity score model is usually assessed based on the resulting covariate balance. In the context of propensity score weighting, this is measured based on either the absolute standardized difference (ASD): \[\label{eq:ASD1} \text{ASD} = \left| {\frac{\sum_{i=1}^N w_1(\mathbf{x}_i)Z_iX_{pi} }{\sum_{i=1}^N w_1(\mathbf{x}_i)Z_i } - \frac{\sum_{i=1}^N w_0(\mathbf{x}_i)(1- Z_i)X_{pi} }{\sum_{i=1}^N w_0(\mathbf{x}_i)(1-Z_i) }}\right| \Bigg /{\sqrt{\frac{s_{1}^2 + s_{0}^2}{2}}}, \tag{3}\] or the target population standardized difference (PSD), \(\max\{\text{PSD}_0,\text{PSD}_1\}\), where \[\label{eq:PSD1} \text{PSD}_z = \left|{\frac{\sum_{i=1}^N w_z(\mathbf{x}_i)\mathbb{1}\{Z_i=z\}X_{pi} }{\sum_{i=1}^N w_z(\mathbf{x}_i)\mathbb{1}\{Z_i=z\} } - \frac{\sum_{i=1}^N h(\mathbf{x}_i)X_{pi} }{\sum_{i=1}^N h(\mathbf{x}_i)}}\right|\Bigg /{\sqrt{\frac{s_{1}^2 + s_{0}^2}{2}}}. \tag{4}\] In (3) and (4), \(s_z^2\) is the variance (either unweighted or weighted, depending on user specification) of the \(p\)th covariate in group \(z\), and \((w_0,w_1)\) are the specified balancing weights. Setting \(w_0=w_1=1\) corresponds to the unweighted mean differences. ASD and PSD are often displayed as column in the baseline characteristics table (known as the “Table 1”) and visualized via a Love plot (also known as a forest plot) (Greifer 2018). A rule of thumb for determining adequate balance is when ASD of all covariates is controlled within \(0.1\) (Austin and Stuart 2015).

Multiple treatments

(Li and Li 2019) extend the framework of balancing weights to multiple treatments. Assume that we have \(J\) \((J\geq 3)\) treatment groups, and let \(Z_i\) stand for the treatment received by unit \(i\), \(Z_i\in \{1,\ldots,J\}\). We further define \(D_{ij}=\mathbb{1}\{Z_i=j\}\) as a set of multinomial indicator, satisfying \(\sum_{i=1}^J D_{ij}=1\) for all \(j\). Denote the potential outcome for unit \(i\) under treatment \(j\) as \(Y_{i}(j)\), of which only the one corresponding to the received treatment, \(Y_i=Y_i(Z_i)\), is observed. The generalized propensity score is the probability of receiving a potential treatment \(j\) given \(\mathbf{X}\) (Imbens 2000): \(e_j(\mathbf{x})=P(Z=j|\mathbf{X}=\mathbf{x})\), with the constraint that \(\sum_{j=1}^J e_j(\mathbf{x})=1\).

To define the target estimand, let \(m_j(\mathbf{x})=\mathop{\mathrm{\mathbb E}}[Y(j)|\mathbf{X}=\mathbf{x}]\) be the conditional expectation of the potential outcome in group \(j\). For specified tilting function \(h(\mathbf{x})\) and target density \(g(\mathbf{x})\propto f(\mathbf{x})h(\mathbf{x})\), the \(j\)th average potential outcome among the target population is \[\label{eq:meanpo} \mu_j^h=\mathop{\mathrm{\mathbb E}}_g[Y(j)]=\frac{\mathop{\mathrm{\mathbb E}}\{h(\mathbf{x})m_j(\mathbf{x})\}}{\mathop{\mathrm{\mathbb E}}\{h(\mathbf{x})\}}. \tag{5}\] Causal estimands can then be constructed in a general manner as contrasts based on \(\mu_j^h\). For example, the most commonly seen estimands in multiple treatments are the pairwise average treatment effects between groups \(j\) and \(j'\): \(\tau_{j,j'}^h=\mu_j^h-\mu_{j'}^h\). This definition can be generalized to arbitrary linear contrasts. Denote \(\pmb{a}=(a_{i},\cdots, a_{J})\) as a contrast vector of length \(J\). A general class of additive estimands is \[\label{eq:meantau} \tau^h(\pmb{a})=\sum\limits_{j=1}^{J}a_j\mu_j^h. \tag{6}\] Specific choices for \(\mathbf{a}\) with nominal and ordinal treatments can be found in Li and Li (2019). Similar to before, propensity score weighting analysis with multiple treatments rests on two assumptions: (1) weak unconfoundedness: \(Y(j)\perp \mathbb{1}\{Z=j\} |\mathbf{X},\) for all \(j\), and (2) Overlap: the generalized propensity score is bounded away from 0 and 1: \(0<e_j(\mathbf{x})<1\), for all \(j\).

With multiple treatments, the tilting function \(h(\mathbf{x})\) specifies the target population, estimand, and balancing weights. For a given \(h(\mathbf{x})\), the balancing weights for the \(j\)th treatment group \(w_j(\mathbf{x})\propto {h(\mathbf{x})}/{e_j(\mathbf{x})}\). Then the Hájek estimator for \(\mu_j^h\) is \[\begin{aligned} \label{eq:simplemean} \hat{\mu}^h_j=\frac{\sum_{i=1}^N w_j(\mathbf{x}_i)D_{ij}Y_i }{\sum_{i=1}^N w_j(\mathbf{x}_i)D_{ij}}. \end{aligned} \tag{7}\] Contrasts based on \(\hat{\mu}^h_j\) can be obtained for any \(\mathbf{a}\) to estimate the additive causal estimand \(\tau^h(\mathbf{a})\). Of note, we only consider types of estimands that are transitive, and therefore the ATT estimands introduced in Lechner (2001) is not implemented. In parallel to binary treatments PSweight implements five types of balancing weights with multiple treatments: IPW, treated weights, OW, MW, and EW, and the corresponding target estimand of each weighting scheme is its pairwise (between each pair of treatments) counterpart in binary treatments.

Among all the weights, OW minimizes the total asymptotic variances of all pairwise comparisons, and has been shown to have the best finite-sample efficiency in estimating pairwise WATEs (Li and Li 2019). Table 3 summarizes the target population, tilting function and balancing weight for multiple treatments that are available in PSweight.

Table 3: Target populations, tilting functions, and the corresponding balancing weights for multiple treatments in PSweight.
Target population Tilting function \(h(\mathbf{x})\) Balancing weights \(\left\{w_j(\mathbf{x}),~j=1,\ldots,J\right\}\)
Combined 1 \(\left\{1/e_j(\mathbf{x})\right\}\)
Treated (\(j'\)th group) \(e_{j'}(\mathbf{x})\) \(\left\{e_{j'}(\mathbf{x})/e_j(\mathbf{x})\right\}\)
Overlap \(\{\sum_{k=1}^J1/e_k(\mathbf{x})\}^{-1}\) \(\left\{\{\sum_{k=1}^J1/e_k(\mathbf{x})\}^{-1}/e_j(\mathbf{x})\right\}\)
Matching \(\min_k \{e_k(\mathbf{x})\}\) \(\left\{ \min_k \{e_k(\mathbf{x})\}/e_j(\mathbf{x}) \right\}\)
Entropy \(-\sum_{k=1}^J e_k(\mathbf{x})\log\{e_k(\mathbf{x})\}\) \(\left\{-\sum_{k=1}^J e_k(\mathbf{x})\log\{e_k(\mathbf{x})\}/ e_j(\mathbf{x}) \right\}\)

To estimate the generalized propensity scores for multiple treatments, the default model in PSweight is a multinomial logistic model. PSweight also allows for externally estimated generalized propensity scores. Goodness-of-fit of the generalized propensity score model is assessed by the resulting covariate balance, which is measured by the pairwise versions of the ASD and PSD. The detailed formula of these metrics can be found in (Li and Li 2019). A common threshold for balance is that the maximum pairwise ASD or maximum PSD is below 0.1.

Propensity score trimming

Propensity score trimming excludes units with estimated (generalized) propensity scores close to zero (or one). It is a popular approach to address the extreme weights problem of IPW. PSweight implements the symmetric trimming rules in Crump et al. (2009) and Yoshida et al. (2018). Operationally, we allow users to specify a single cutoff \(\delta\) on the estimated generalized propensity scores, and only includes units for analysis if \(\min_j\{e_{j}(\mathbf{x})\}\in[\delta,1]\). With binary treatments, the symmetric trimming rule reduces to \(e(\mathbf{x})\in[\delta,1-\delta]\). The natural restriction \(\delta<1/J\) must be satisfied due to the constraint \(\sum_{j=1}^J e_j(\mathbf{x})=1\). To avoid specifying an arbitrary trimming threshold \(\delta\), PSweight also implements the optimal trimming rules of Crump et al. (2009) and Yang et al. (2016), which minimizes the (total) asymptotic variance(s) for estimating the (pairwise) ATE among the class of all trimming rules. OW can be viewed as a continuous version of trimming because it smoothly down-weigh the units with propensity scores close to 0 or 1, and thus avoids specifying a threshold.

Augmented weighting estimators

PSweight also implements augmented weighting estimators, which augment a weighting estimator by an outcome regression and improves the efficiency. With IPW, the augmented weighting estimator is known as the doubly-robust estimator (Lunceford and Davidian 2004; Bang and Robins 2005; Funk et al. 2011). With binary treatments, the augmented estimator with general balancing weights are discussed Hirano et al. (2003) and Mao et al. (2018). Below, we briefly outline the form of this estimator with multiple treatments. Recall the conditional mean of \(Y_i(j)\) given \(\mathbf{X}_i\) and treatment \(Z_i=j\) as \(m_{j}(\mathbf{x}_i)=\mathop{\mathrm{\mathbb E}}[Y_i(j)|\mathbf{X}_i=\mathbf{x}_i]=\mathop{\mathrm{\mathbb E}}[Y_i|\mathbf{X}_i=\mathbf{x}_i,Z_{i}=j]\). This conditional mean can be estimated by generalized linear models, kernel estimators, or machine learning models. PSweight by default employs the generalized linear models, but also allows estimated values from other routines. When \(m_{j}(\mathbf{x}_i)\) is estimated by generalized linear models, PSweight currently accommodates three types of outcomes: continuous, binary and count outcomes (with or without an offset), using the canoncal link function.

With a pre-specified tilting function, the augmented weighting estimator for group \(j\) is \[\begin{aligned} \label{eq:augest} \hat{\mu}^{h,\tiny{aug}}_j=\frac{\sum_{i=1}^N w_j(\mathbf{x}_i)D_{ij}\{Y_i-m_{j}(\mathbf{x}_i)\} }{\sum_{i=1}^N w_j(\mathbf{x}_i)D_{ij}}+\frac{\sum_{i=1}^N h(\mathbf{x}_i)m_{j}(\mathbf{x}_i) }{\sum_{i=1}^N h(\mathbf{x}_i)}. \end{aligned} \tag{8}\] The first term of ((8)) is the Hájek estimator of the regression residuals, and the second term is the standardized average potential outcome (a \(g\)-formula estimator). With IPW, ((8)) is consistent to \(\mathop{\mathrm{\mathbb E}}[Y(j)]\) when either the propensity score model or the outcome model is correctly specified, but not necessarily both. For other balancing weights, ((8)) is consistent to the WATE when the propensity model is correctly specified, regardless of outcome model specification. When both models are correctly specified, ((8)) achieves the lower bound of the variance for regular and asymptotic linear estimators (Robins et al. 1994; Hirano et al. 2003; Mao et al. 2018).

Ratio causal estimands

With binary and count outcomes, ratio causal estimands are often of interest. Using notation from the multiple treatments as an example, once we use weighting to obtain estimates for the set of average potential outcomes \(\{\mu_j^h,j=1,\ldots,J\}\), we can directly estimate the causal relative risk (RR) and causal odds ratio (OR), defined as \[\label{eq:RROR} \tau^{h,\tiny{RR}}_{j,j^{\prime}}=\frac{\mu_{j}^h}{\mu_{j^{\prime}}^h},~~~~~~\tau^{h,\tiny{OR}}_{j,j^{\prime}}=\frac{\mu_{j}^h/(1-\mu_{j}^h)}{\mu_{j^{\prime}}^h/(1-\mu_{j^{\prime}}^h)}. \tag{9}\] Here the additive estimand \(\tau^{h,\tiny{RD}}_{j,j^{\prime}}=\mu_{j}^h-\mu_{j^{\prime}}^h\) is the causal risk difference (RD). PSweight supports a class of ratio estimands for any given contrasts \(\mathbf{a}\). Specifically, we define the log-RR type parameters by \[\label{eq:meanRR} \lambda^{h,\tiny{RR}}(\pmb{a})=\sum\limits_{j=1}^{J}a_j\log\left(\mu_j^h\right), \tag{10}\] and the log-OR type parameters by \[\label{eq:meanOR} \lambda^{h,\tiny{OR}}(\pmb{a})=\sum\limits_{j=1}^{J}a_j\left\{\log\left(\mu_j^h\right)-\log\left(1-\mu_j^h\right)\right\}. \tag{11}\] With nominal treatments, the contrast vector \(\mathbf{a}\) can be specified to encode pairwise comparisons in the log scale (as in (10)) or in the log odds scale (as in (11)), in which case \(\exp\{\lambda^{h,\tiny{RR}}(\pmb{a})\}\) and \(\exp\{\lambda^{h,\tiny{OR}}(\pmb{a})\}\) become the causal RR and causal OR in (9). User-specified contrasts \(\mathbf{a}\) can provide a variety of nonlinear estimands. For example, when \(J=3\), with \(\mathbf{a}=(1,-2,1)^T\) one can use PSweight to assess the equality of two consecutive causal RR: \(H_0: \mu_3^h/\mu_2^h=\mu_2^h/\mu_1^h\).

Variance and interval estimation

PSweight by default implements the empirical sandwich variance for propensity score weighting estimators (Lunceford and Davidian 2004; Mao et al. 2018; Li et al. 2019) based on the M-estimation theory (Stefanski and Boos 2002). The variance adjusted for the uncertainty in estimating the propensity score and outcome models, and are sometime referred to as the nuisance-adjusted sandwich variance. Below we illustrate the main steps with multiple treatments and general balancing weights. Write \(\mathbf{\theta}=\left(\nu_1,\ldots,\nu_J,\eta_1,\ldots,\eta_J,\mathbf{\beta}^T,\mathbf{\alpha}^T\right)^T\) as the collection of parameters to be estimated. Then \(\left\{\hat{\mu}^{h,\tiny{aug}}_j=\hat{\nu}_j+\hat{\eta}_j:j=1,\ldots,J\right\}\) jointly solve \[\begin{aligned} \sum_{i=1}^{N}\Psi_{i}(\mathbf{\theta})=\sum_{i=1}^{N} \left( \begin{array}{c} w_1(\mathbf{x}_i)D_{i1}\{Y_{i}-m_{1}(\mathbf{x}_{i};\mathbf{\alpha})-\nu_1\}\\ \vdots\\ w_J(\mathbf{x}_i)D_{iJ}\{Y_{i}-m_{J}(\mathbf{x}_{i};\mathbf{\alpha})-\nu_J\}\\ h(\mathbf{x}_i)\{m_{1}(\mathbf{x}_{i};\mathbf{\alpha})-\eta_1\}\\ \vdots\\ h(\mathbf{x}_i)\{m_{J}(\mathbf{x}_{i};\mathbf{\alpha})-\eta_J\}\\ S_{\mathbf{\beta}}(Z_i,\mathbf{x}_{i};\mathbf{\beta})\\ S_{\mathbf{\alpha}}(Y_i,Z_i,\mathbf{x}_{i};\mathbf{\alpha}) \end{array} \right)=\mathbf{0}, \end{aligned}\] where \(S_{\mathbf{\beta}}(Z_i,\mathbf{x}_{i};\mathbf{\beta})\) and \(S_{\mathbf{\alpha}}(Y_i,Z_i,\mathbf{x}_{i};\mathbf{\alpha})\) are the score functions of the propensity score model and the outcome model. The empirical sandwich variance estimator is \[\begin{aligned} \widehat{\mathop{\mathrm{\mathbb V}}}(\hat{\mathbf{\theta}})=\left\{\sum_{i=1}^{N} \frac{\partial}{\partial\mathbf{\theta}^T}\Psi_{i}(\hat{\mathbf{\theta}})\right\}^{-1} \left\{\sum_{i=1}^{N}\Psi_{i}(\hat{\mathbf{\theta}}) \Psi_{i}^{T}(\hat{\mathbf{\theta}})\right\} \left\{\sum_{i=1}^{N} \frac{\partial}{\partial\mathbf{\theta}}\Psi_{i}^T(\hat{\mathbf{\theta}})\right\}^{-1}. \end{aligned}\] Because \(\hat{\mu}^{h,\tiny{aug}}_j=\hat{\nu}_j+\hat{\eta}_j\), the variance of arbitrary linear contrasts based on the average potential outcomes can be easily computed by applying the Delta method to the joint variance \(\widehat{\mathop{\mathrm{\mathbb V}}}(\hat{\mathbf{\theta}})\). For the Hájek weighting estimators, variance is estimated by removing \(S_{\mathbf{\alpha}}(Y_i,Z_i,\mathbf{x}_{i};\mathbf{\alpha})\) as well as the components involving \(m_j(\mathbf{x}_i;\mathbf{\alpha})\) in \(\Psi_{i}(\mathbf{\theta})\). Finally, when propensity scores and potential outcomes are not estimated through the generalized linear model or are supplied externally, or MW are used (since the tilting function is not everywhere differentiable), PSweight ignores the uncertainty in estimating \(\mathbf{\beta}\) and \(\mathbf{\alpha}\) and removes \(S_{\mathbf{\beta}}(Z_i,\mathbf{x}_{i};\mathbf{\beta})\) and \(S_{\mathbf{\alpha}}(Y_i,Z_i,\mathbf{x}_{i};\mathbf{\alpha})\) in \(\Psi_{i}(\mathbf{\theta})\) in the calculation of the empirical sandwich variance. Based on the estimated variance, PSweight computes the associated symmetric confidence intervals and p-values via the normal approximation.

For ratio causal estimands, PSweight applies the logarithm transformation to improve the accuracy of the normal approximation (Agresti 2003). For estimating the variance of causal RR, we first obtain the joint variance of \(\left(\log\left(\hat{\mu}^{h,\tiny{aug}}_1\right),\ldots \log\left(\hat{\mu}^{h,\tiny{aug}}_J\right)\right)^T\) using the Delta method, and then estimate the variance of \(\lambda^{h,\tiny{RR}}(\pmb{a})\). Once the symmetric confidence intervals are obtained for \(\lambda^{h,\tiny{RR}}(\pmb{a})\) using the normal approximation, we can exponentiate the upper and lower confidence limits to derive the asymmetric confidence intervals for the causal RR. Confidence intervals for the causal OR are computed similarly.

PSweight also allows using bootstrap to estimate variances, which can be much more computationally intensive than the closed-form sandwich estimator but sometimes give better finite-sample performance in small samples. By default, PSweight resamples \(R=50\) bootstrap replicates with replacement. For each replicate, the weighting estimator ((7)) or the augmented weighting estimtor ((8)) is implemented, providing \(R\) estimates of the \(J\) average potential outcomes (an \(R\times J\) matrix). Then for any contrast vector \(\pmb{a}=(a_{1},\cdots, a_{J})^T\), PSweight obtains \(R\) bootstrap estimates: \[\hat{\mathbb{T}}^h(\pmb{a})_{bootstrap}= \left\{\hat{\tau}^h(\pmb{a})_1=\sum_{j=1}^{J} a_{j} \hat{\mu}^h_{j,1},~\ldots~, \hat{\tau}^h(\pmb{a})_R=\sum_{j=1}^{J} a_{j} \hat{\mu}^h_{j,R} \right\}.\] The sample variance of \(\hat{\mathbb{T}}^h(\pmb{a})_{bootstrap}\) is reported by PSweight as the bootstrap variance; the lower and upper \(2.5\%\) quantiles of \(\hat{\mathbb{T}}^h(\pmb{a})_{bootstrap}\) form the \(95\%\) bootstrap interval estimate

3 Overview of package

The PSweight package includes two modules tailored for design and analysis of observational studies. The design module provides diagnostics to assess the adequacy of the propensity score model and the weighted target population, prior to the use of outcome data. The analysis module provides functions to estimate the causal estimands. We briefly describe the two modules below.

Design module

PSweight offers the SumStat() function to visualize the distribution of the estimated propensity scores, to assess the balance of covariates under different weighting schemes, and to characterize the weighted target population. It uses the following code snippet:

SumStat(ps.formula, ps.estimate = NULL, trtgrp = NULL, Z = NULL, covM = NULL, 
    zname = NULL,  xname = NULL, data = NULL,weight = "overlap", delta = 0,
    method = "glm", ps.control = list())

By default, the (generalized) propensity scores are estimated by the (multinomial) logistic regression, through the argument ps.formula. Alternatively, gbm() functions in the gbm package (Greenwell et al. 2019) or the SuperLearner() function in the SuperLearner package (Polley et al. 2019) can also be called by using method = "gbm" or method = "SuperLearner". Additional parameters of those functions can be supplied through the ps.control argument. The argument ps.estimate supports estimated propensity scores from external routines. SumStat() produces a SumStat object, with estimated propensity scores, unweighted and weighted covariate means for each treatment group, balance diagnostics, and effective sample sizes (defined in (Li and Li 2019)). We then provide a summary.SumStat() function, which takes the SumStat object and summarizes weighted covariate means by treatment groups and the between-group differences in either ASD or PSD. The default options in weighted.var = TRUE and metric = "ASD" yield ASD based on weighted standard deviations in (Austin and Stuart 2015). The weighted covariate means can be used to build a baseline characteristics “Table 1” to illustrate the target population where trimming or balancing weights are applied.

summary(object, weighted.var = TRUE, metric = "ASD")
Table 4: Functions in the design module of PSweight.
Function Description
SumStat() Generate a SumStat object with information of propensity scores and weighted covariate balance
summary.SumStat() Summarize the SumStat object and return weighted covariate means by treatment groups and weighted or unweighted between-group differences in ASD or PSD
plot.SumStat() Plot the distribution of propensity scores or weighted covariate balance metrics from the SumStat object
PStrim() Trim the data set based on estimated propensity scores

Diagnostics of propensity score models can be visualized with the plot.SumStat() function. It takes the SumStat object and produces a balance plot (type = "balance") based on the ASD and PSD. A vertical dashed line can be set by the threshold argument, with a default value equal to \(0.1\). The plot.SumStat() function can also supply density plot (type = "density"), or histogram (type = "hist") of the estimated propensity scores. The histogram, however, is only available for the binary treatment case. The plot function is implemented as follows:

plot(x, type = "balance", weighted.var = TRUE, threshold = 0.1, metric = "ASD")

In the design stage, propensity score trimming can be carried out with the PStrim() function. The trimming threshold delta is set to 0 by default. PStrim() also enables optimal trimming rules (optimal = TRUE) that give the most statistically efficient (pairwise) subpopulation ATE, among all possible trimming rules. A trimmed data set along with a summary of trimmed cases will be returned by PStrim(). This function is given below:

PStrim(data, ps.formula = NULL, zname = NULL, ps.estimate = NULL, delta = 0, 
    optimal = FALSE, method = "glm", ps.control = list())

Alternatively, trimming is also anchored in the SumStat() function with the delta argument. All functions in the design module are summarized in Table 4.

Analysis module

The analysis module of PSweight includes two functions: PSweight() and summary.PSweight(). The PSweight() function estimates the average potential outcomes in the target population, \(\{\mu_{j}^{h},j=1,\ldots,J\}\), and the associated variance-covariance matrix. By default, the empirical sandwich variance is implemented, but bootstrap variance can be obtained with the argument bootstrap = TRUE). The weight argument can take "IPW", "treated", "overlap", "matching" or "entropy", corresponding to the weights introduced in Tables 2 and 3. More detailed descriptions of each input argument in the PSweight() function can be found in Table 5. A typical PSweight() code snippet looks like

PSweight(ps.formula, ps.estimate, trtgrp, zname, yname, data, weight = "overlap", 
    delta = 0, augmentation = FALSE, bootstrap = FALSE, R = 50, out.formula = NULL,
    out.estimate = NULL, family = "gaussian", ps.method = "glm", ps.control = list(),
    out.method = "glm", out.control = list())
Table 5: Arguments for function PSweight() in the analysis module of PSweight.
Argument Description Default
ps.formula A symbolic description of the propensity score model.
ps.estimate An optional matrix or data frame with externally estimated (generalized) propensity scores for each observation; can also be a vector with binary treatments. NULL
trtgrp An optional character defining the treated population for estimating (pairwise) ATT. It can also be used to specify the treatment level when only a vector of values are supplied for ps.estimate in the binary treatment setting. Last value in alphabetic order
zname An optional character specifying the name of the treatment variable when ps.formula is not provided. NULL
yname A character specifying name of the outcome variable in data.
weight A character specifying the type of weights to be used. "overlap"
delta Trimming threshold for (generalized) propensity scores. 0
augmentation Logical value of whether augmented weighting estimators should be used. FALSE
bootstrap Logical value of whether bootstrap is used to estimate the standard error FALSE
R Number of bootstrap replicates if bootstrap = TRUE 50
out.formula A symbolic description of the outcome model to be estimated when augmentation = TRUE
out.estimate An optional matrix or data frame containing externally estimated potential outcomes for each observation under each treatment level. NULL
family A description of the error distribution and canonical link function to be used in the outcome model if out.formula is provided "gaussian"
ps.method a character to specify the method for propensity model. "glm"
ps.control A list to specify additional options when method is set to "gbm" or "SuperLearner". list()
out.method A character to specify the method for outcome model. "glm"
out.control A list to specify additional options when methodout is set to "gbm" or "SuperLearner". list()

Similar to the design module, the summary.PSweight() function synthesizes information from the PSweight object for statistical inference. A typical code snippet looks like

summary(object, contrast, type = "DIF", CI = TRUE)

In the summary.PSweight() function, the argument type corresponds to the three types estimands: type = "DIF" is the default argument that specifies the additive causal contrasts; type = "RR" specifies the contrast on the log scale as in equation (10); type = "OR" specifies the contrast on the log odds scale as in equation (11). Confidence intervals and p-values are obtained using normal approximation and reported by the summary.PSweight() function. The argument contrast represents a contrast vector \(\pmb{a}\) or matrix with multiple contrast row vectors. If contrast is not specified, summary.PSweight() provides all pairwise comparisons of the average potential outcomes. By default, confidence interval is printed (CI = TRUE); alternatively, one can print the test statistics and p-values by CI = FALSE.

4 Case study with the NCDS data

We demonstrate PSweight in a case study that estimates the causal effect of educational attainment on hourly wage, based on the National Child Development Survey (NCDS) data. The National Child Development Survey (NCDS) is a longitudinal study on children born in the United Kingdom (UK) in \(1958\) 1. NCDS collected information such as educational attainment, familial backgrounds, and socioeconomic and health well being on \(17,415\) individuals. We followed Battistin and Sianesi (2011) to pre-process the data and obtain a subset of 3,642 males employed in 1991 with complete educational attainment and wage information for analysis. For illustration, we use the Multiple Imputation by Chained Equations in (Buuren and Groothuis-Oudshoorn 2010) to impute missing covariates and obtain a single imputed data set for all subsequent analysis.2 The outcome variable wage is log of the gross hourly wage in Pound. The treatment variable is educational attainment. For the multiple treatment case, To start with, we created Dmult as a treatment variable with three levels: ">=A/eq", "O/eq" and "None", representing advanced qualification (\(1,806\) individuals), intermediate qualification (\(941\) individuals) and no qualification (\(895\) individuals). We consider twelve pre-treatment covariates or potential confounders. The variable white indicates whether an individual identified himself as white race; scht indicates the school type they attended at age \(16\); qmab and qmab2 are math test scores at age \(7\) and \(11\); qvab and qvab2 are two reading test scores at age \(7\) and \(11\); sib\(\_\)u stands for the number of siblings; agepa and agema are the ages of parents in year 1,974; in the same year, the employment status of mother maemp was also collected; paed\(\_\)u and maed\(\_\)u are the years of education for parents. For simplicity, we will focus on IPW and the three types of weights that improve covariate overlap: OW, MW and EW.

Estimating generalized propensity scores and balance assessment

We use Dmult, the three-level variable, as the treatment of interest. About one half of the population attained advanced academic qualification, there are approximately equal numbers of individuals with intermediate academic qualification or no academic qualification. To illustrate the estimation and inference for ratio estimands, we also introduce a binary outcome of wage, wagebin. The dichotomized wage was obtained with the cutoff of the average hourly wage of actively employed British male aged \(30\)-\(39\) in \(1991\)3. The averaged hourly wage is \(8.23\), and we take \(\log(8.23)\approx 2.10\) as the cutoff. Among the study participants, we observe \(1610\) and \(2032\) individuals above and below the average, and we are interested in estimating the pairwise (weighted) average treatment effect of the academic qualification for obtaining above-average hourly wage.

We specify a multinominal regression model, ps.mult, to estimate the generalized propensity scores.

ps.mult <- Dmult ~ white + maemp + as.factor(scht) + as.factor(qmab)
    as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u + maed_u + 
    agepa + agema + sib_u + paed_u * agepa + maed_u * agema

Then we obtain the propensity score estimates and assess weighted covariate balance with the SumStat() function.

bal.mult <-  SumStat(ps.formula = ps.mult,  
    weight = c("IPW", "overlap", "matching", "entropy"), data = NCDS)
plot(bal.mult, type = "density")

graphic without alt text graphic without alt text graphic without alt text

Figure 1: Density plots of estimated generalized propensity scores with respect to the three-level treatment variable Dmult generated by plot.SumStat() function in the PSweight package.

The distributions of generalized propensity scores are given in Figure 1 (in alphabetic order of the names of treatment groups). For the generalized propensity score to receive the advanced qualification (">=A/eq") or no qualification ("None"), there is a mild lack of overlap due to separation of the group-specific distribution. Since bal.mult includes four weighting schemes, we plot the maximum pairwise ASD and assess the (weighted) covariate balance in a single Love plot.

plot(bal.mult, metric = "ASD")

The covariates are imbalanced across the three groups prior to any weighting. Although IPW can generally improve covariate balance, the maximum pairwise ASD still ocassionally exceeds the threshold \(0.1\) due to lack of overlap. In contrast, OW, MW and EW all emphasize the subpopulation with improved overlap and provide better balance across all covariates.