psw is the main function to perfrom propensity score weighting analysis for (1) visualization of the propensity score distribution in both treatment groups,
(2) covariate balance diagnosis, (3) propensity score model specification test, (4) treatment effect estimation and inference, and (5) augmented estimation with outcome regression
when applicable.
psw(data, form.ps, weight, std.diff = FALSE, V.name = NULL,
mirror.hist = FALSE, add.weight = FALSE, nclass = 50, wt = FALSE,
out.var = NULL, family = "gaussian", aug = FALSE, form.outcome = NULL,
spec.test = F, trans.type = NULL, K = 4)data frame to be used.
propensity score model.
weighting method to be used. Available methods are "ATE", "ATT", "ATC", "MW", "OVERLAP", and "TRAPEZOIDAL".
calculate standardized mean difference as a percentage, std.diff=FALSE by default.
a vector of covariates on which standardized mean difference is computed or the specification test is performed. If V.name = NULL, the covariates in propensity score model are used.
mirror histogram showing the propensity score distributions in both treatment groups, mirror.hist=FALSE by default.
add propensity score weights to the mirror histogram, add.weight=FALSE by default and it is not available for weight="ATE", "ATT", or "ATC".
number of breaks in the mirror histogram.
estimate the weighted estimator, wt=FALSE by default.
outcome variable, needed when wt=TRUE.
outcome family, either "gaussian" or "binomial", family="gaussian" by default.
estimate the augmented estimator, aug=FALSE by default.
outcome model, needed when aug=TRUE.
propensity score model specification test, spec.test=FALSE by default.
Note that specification test is not available for weight="OVERLAP".
a vector of the same length as V.name that specifies the transformation type for each element in V.name.
Available transformations are "identity", "log", "logit", "sqrt", "Fisher".
Needed when spec.test=T, and no transformation is perfomred with "identity". See Details.
value of \(K\) in \(\omega(e_i) = min(1, K min(e_i, 1-e_i)) \) for "TRAPEZOIDAL" weight. The estimand is
closer to the average treatment effect (ATE) with larger value of K. K=4 by default.
psw returns a list of elements depending on the supplied arguments.
weighting method.
object returned by fitting the propensity score model using glm with "binomial" family.
estimated propensity score.
estimated propensity score weight.
A data frame of weighed mean, variance, and standardized mean difference for covariates in V.name (see below) by treatment groups before weighting.
V.name is the row name and "treated.mean", "treated.var", "control.mean", "control.var", "std.diff.pct" are column names.
A data frame of weighed mean, variance, and standardized mean difference for covariates in V.name by treatment groups after weighting.
weighted estimator for mean difference when wt=T and family = "gaussian".
standard error for est.wt.
augmented estimator for mean difference when aug=T and family = "gaussian".
standard error for est.aug.
weighted estimator for risk difference when wt=T and family = "binomial".
standard error for est.risk.wt.
augmented estimator for risk difference when aug=T and family = "binomial".
standard error for est.risk.aug.
weighted estimator for relative risk when wt=T and family = "binomial".
standard error for est.rr.wt.
weighted estimator for odds ratio when wt=T and family = "binomial".
standard error for est.or.wt.
weighted estimator for log odds ratio when wt=T and family = "binomial".
standard error for est.lor.wt.
covariates for balance diagnosis and specification test.
test statistic for the specification test, which follows the \(\chi^2_{df}\) distribution under the null, available when spec.test=T.
degree of freedom for the specification test, df=rank(var.B.hat), available when spec.test=T.
pvalue of the specification test when spec.test=T.
In package PSW, treatment indicator (left handside of form.ps) should be dummy coded
such that a value of 1 indicates the treated and a value of 0 indicates the control. All categorical covariates need to be dummy coded too.
If the outcome belongs to the "gaussian" family, causal estimation is based on the mean differnce between treatment groups. If the
outcome belongs to the "binomial" family, causal estimation is based on risk difference, risk ratio, odds ratio or log odds ratio.
The Delta method is used for variance estimation when applicable.
Let \(Z_i\) be the treatment indicator of subject \(i\), \(e_i\) be the corresponding propensity score. Then
propensity score weight, \(W_i\), is defined as $$W_i = \frac{\omega(e_i)}{Z_i e_i + (1-Z_i)(1-e_i)},$$ where \(\omega(e_i)\) is a function depending
on \(e_i\). For "ATE", \(\omega(e_i) = 1\), which leads to estimating the average treatment effect. For "ATT", \(\omega(e_i) = e_i\),
which leads to estimating average treatment effect for the treated. For "ATC", \(\omega(e_i) = 1 - e_i\), which leads to estimating average treatment effect
for the controls. For "MW", \(\omega(e_i) = min( e_i, 1 - e_i )\). For "OVERLAP", \(\omega(e_i) = e_i(1 - e_i)\). For "TRAPEZOIDAL",
\(\omega(e_i) = min( 1, K min( e_i, 1 - e_i ) )\). This type of weights are studied by Hirano, Imbens and Ridder (2003) and Li et al (2016).
The \(\omega(e_i)\) function is specified by the weight argument.
The matching weight ("MW") was proposed by Li and Greene (2013). The overlap weight ("OVERLAP") was propsed by Li et al (2016).
These methods down weight subjects with propensity score close to 0 or 1. and hence improve the stability of computation.
A mirror histogram is produced to visualize the propensity score distributions in both treatment groups. In the mirror histogram, above the horizontal line
is the histogram of the propensiy scores of the control group, below is that of the treated group. The vertical axis of the histogram is the frequency. When
add.weight=TRUE, the height of the green bar added to mirror histogram is the summation of the weights of subjects within the corresponding propensity
score stratum. For weighting methods of "ATE", "ATT", "ATC", add.weight is not recommended for visualization because weights may
be larger than 1.
Standardized mean difference for a covariate is defiend as
$$ \frac{100 (\bar{x}_1 - \bar{x}_0)}{\sqrt{\frac{s_1^2 + s_0^2}{2} } },$$
where \(\bar{x}_1\) and \(s_1^2\) are weighted mean and standard deviation for the treated group, and \( \bar{x_0}\) and \(s_0^2\)
are defined similarly for the control group. A plot showing the standardized mean difference before and after weighting adjustement will be generated to
facilitate covariate balance diagnosis. It is sometimes recommended that the absolute values of standardized mean differences of all covariates should be less
than 10% in order to claim covariate balance.
For the proensity score model specification test (Li and Greene, 2013), the quantity of interest is
$$ \hat{B} = \boldsymbol{g} \left\{ \frac{ \sum_{i=1}^n W_i Z_i \boldsymbol{V}_i}{\sum_{i=1}^n W_i Z_i}\right\} - \boldsymbol{g} \left\{ \frac{ \sum_{i=1}^n W_i (1-Z_i) \boldsymbol{V}_i}{\sum_{i=1}^n W_i (1-Z_i)}\right\}, $$
where \(\boldsymbol{V}_i\) is a vector of covariates whose balance are examined, and \(\boldsymbol{g}(.)\) is a vector of monotone smooth transformations for the input.
Transformation type is specified by argument trans.type, and available transformation types are "identity", "log", "logit", "sqrt", "Fisher".
These transformations are recommended to improve the finite sample performance of the specification test. Log transformation ("log") and square root transformation ("sqrt")
are recommended for skewed data, logit transformation ("logit") for binary data, and Fisher z-transformation ("Fisher") for bounded data between \((-1, 1)\).
The current version of model specification test is not available for weight="OVERLAP" because it results in zero standardized difference.
For estimation of mean difference ("gaussian" family) or risk difference ("binomial" family), the weighted estimator is
$$ \hat{\Delta} = \frac{\sum_{i=1}^n W_i Z_i Y_i}{\sum_{i=1}^n W_i Z_i} - \frac{\sum_{i=1}^n W_i (1-Z_i) Y_i}{\sum_{i=1}^n W_i (1-Z_i)}, $$
and the augmented estimator is
$$\hat{\Delta}_{aug} = \frac{ \sum_{i=1}^n \omega(e_i) \{ m_{1i} - m_{0i} \}}{ \sum_{i=1}^n \omega(e_i) } + \frac{ \sum_{i=1}^n W_i Z_i \{ Y_i - m_{1i} \}}{ \sum_{i=1}^n W_i Z_i } - \frac{ \sum_{i=1}^n W_i (1-Z_i) \{ Y_i - m_{0i} \}}{ \sum_{i=1}^n W_i (1-Z_i)},$$
where \(m_{1i} = E[Y_i | \boldsymbol{X_i}, Z_i=1]\) is the conditional expectation of outcome when treated given covariates \(\boldsymbol{X}_i\),
and \(m_{0i} = E[Y_i | \boldsymbol{X_i}, Z_i=0]\) is the conditional expectation of outcome when control given covariates \(\boldsymbol{X}_i\).
When the outcome belongs to the "binomial" family, the marginal probability is used to estimate risk ratio, odds ratio and log odds ratio.
Sandwich variance estimation is used to adjust for the sampling variability in the estimated propensity scores (Li and Greene, 2013).
The augmented estimator \( \hat{\Delta}_{aug} \) incorporates regression models for the outcome variable and has simliar properties as the doubly robust IPW estimator (Lunceford and Davidian, 2004), but with one difference. The estimand of IPW estimator does not depend on the propensity score because \(\omega(e_i) = 1\), while the estimands of other weighting methods depend on propensity score specification. Nonetheless, the proposed augmented estimator converges to the estimand defined by the corresponding propensity score model.
Hirano K, Imbens GW and Ridder G. "Efficient estimation of average treatment effects using the estimated propensity score." Econometrica 2003; 71(4): 1161-1189.
Li F, Morgan KL and Zaslavsky AM. "Balancing covariates via propensity score weighting." J Am Stat Assoc 2016; DOI:10.1080/01621459.2016.1260466.
Li L and Greene T. "A weighting analogue to pair matching in propensity score analysis." Int J Biostat 2013; 9(2):215-234.
Lunceford JK and Davidian M. Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study. Stat Med. 2004; 23(19):2937-2960.
psw.balance, psw.spec.test, psw.wt, psw.aug, psw.mirror.hist.
# NOT RUN {
# Load the test data set
data(test_data);
# Propensity score model
form.ps <- "Z ~ X1 + X2 + X3 + X4";
# A vector of covariates
V.name <- c( "X1", "X2", "X3", "X4" );
#1. Standardized differnce with "ATE"
tmp1 <- psw( data = test_data, form.ps = form.ps, weight = "ATE",
std.diff = TRUE, V.name = V.name );
#2. Mirror histogram and add estimated matching weight to it
tmp2 <- psw( data = test_data, form.ps = form.ps, weight = "MW",
mirror.hist = TRUE, add.weight = TRUE );
#3. Estimate average treatment effect with "ATE"
tmp3 <- psw( data = test_data, form.ps = form.ps, weight = "ATE", wt = TRUE,
out.var = "Y", family = "gaussian" );
#4. Augmented estimator with "OVERLAP"
# outcome model
form.out <- "Y ~ X1 + X2 + X3 + X4";
tmp4 <- psw( data = test_data, form.ps = form.ps, weight = "OVERLAP", aug = TRUE,
form.outcome = form.out, family = "gaussian" );
#5. Propensity score model specification test with "MW".
# A vector of transformation types for covariates in V.name.
trans.type <- c( "identity", "identity", "logit", "logit" );
tmp5 <- psw( data = test_data, form.ps = form.ps, weight = "MW", spec.test = TRUE,
V.name = V.name, trans.type = trans.type );
# }
Run the code above in your browser using DataLab