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) nuisanceadjusted 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.
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 onestop 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.
Multiple  Balance  IPW/ATT  OW/other  Ratio  Augmented  Nuisanceadj  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.13): 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 weightingbased estimation with multiple treatments, but focus on IPW. The PSW R package (Mao and Li 2018) implements both OW and MW and allows for nuisanceadjusted 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.
Before diving into the implementation details of PSweight, we briefly introduce the basics of the propensity score weighting framework.
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)+(1Z_{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})/\{1e(\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)(1Z_i) Y_i}{\sum_{i=1}^Nw_0(\mathbf{x}_i)(1Z_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 wellknown 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 finitesample) 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.
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}{1e(\mathbf{x})}\right)\) 
Treated  \(e(\mathbf{x})\)  ATT  \(\left(1, \frac{e(\mathbf{x})}{1e(\mathbf{x})}\right)\) 
Overlap  \(e(\mathbf{x})(1e(\mathbf{x}))\)  ATO  \((1e(\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})}{1e(\mathbf{x})} \right)\) 
Entropy  \(\xi_2(\mathbf{x})\)  ATEN  \(\left( \frac{ \xi_2(\mathbf{x})}{e(\mathbf{x})}, \frac{ \xi_2(\mathbf{x})}{1e(\mathbf{x})} \right)\) 
Notes: \(\xi_1(\mathbf{x}) = \min \{e(\mathbf{x}),1e(\mathbf{x})\}\) and \(\xi_2(\mathbf{x})=\{e(\mathbf{x})\log(e(\mathbf{x}))+(1e(\mathbf{x}))\log(1e(\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.
Goodnessoffit 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)(1Z_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).
(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 finitesample 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.
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. Goodnessoffit 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 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 downweigh the units with propensity scores close to 0 or 1, and thus avoids specifying a threshold.
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 doublyrobust 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 prespecified 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_im_{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).
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 logRR 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 logOR 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). Userspecified 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\).
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 Mestimation 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 nuisanceadjusted 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 pvalues 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 closedform sandwich estimator but sometimes give better finitesample 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
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.
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
betweengroup 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")
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 betweengroup 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.
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 variancecovariance
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())
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 pvalues 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
pvalues by CI = FALSE
.
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 preprocess 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 GroothuisOudshoorn 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 pretreatment 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.
We use Dmult
, the threelevel 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 aboveaverage hourly wage.
We specify a multinominal regression model, ps.mult
, to estimate the
generalized propensity scores.
< Dmult ~ white + maemp + as.factor(scht) + as.factor(qmab)
ps.mult as.factor(qmab2) + as.factor(qvab) + as.factor(qvab2) + paed_u + maed_u +
+ agema + sib_u + paed_u * agepa + maed_u * agema agepa
Then we obtain the propensity score estimates and assess weighted
covariate balance with the SumStat()
function.
< SumStat(ps.formula = ps.mult,
bal.mult weight = c("IPW", "overlap", "matching", "entropy"), data = NCDS)
plot(bal.mult, type = "density")
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 groupspecific
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.