Propensity score weighting

Description

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.

Usage

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)

Arguments

Details

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.

Value

psw returns a list of elements depending on the supplied arguments.

References

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.

See Also

psw.balance, psw.spec.test, psw.wt, psw.aug, psw.mirror.hist.

Examples

# 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 );

Propensity score weighting with augmented estimation

Description

psw.aug is the function to estimate the augmented estimator for mean difference (mean outcome difference for "gaussian" family and risk difference for "binomial" family). The augmented estimator is consistent for the estimand defined by the corresponding propensity score model.

Usage

psw.aug(data, form.ps, weight, form.outcome, family = "gaussian", K = 4)

Arguments

Details

psw.aug is used to estimate the augmented estimator, \hat{\Delta}_{aug}, and make inference using the sandwich variance that adjusts for the sampling variability in the estimated propensity score.

Value

A list of weighting method, fitted propensity score model, estimated propenstity scores, estimated propensity score weights, augmented estimator and associated standard error.

Examples

# Load the test data set
data(test_data);
# Propensity score model
form.ps <- "Z ~ X1 + X2 + X3 + X4";
# Outcome model
form.out <- "Y ~ X1 + X2 + X3 + X4";
tmp <- psw.aug( data = test_data, form.ps = form.ps, weight = "ATE",
form.outcome = form.out, family="gaussian" );

Balance checking using standardized mean difference

Description

psw.balance is used to compute the standardized mean difference (in percentage) for balance diagnosis.

Usage

psw.balance(data, form.ps, weight, V.name = NULL, K = 4)

Arguments

Value

A list of weighting method, fitted propensity score model, estimated propenstity scores, estimated propensity score weights, standardized mean difference before and after weighting adjustment.

See Also

psw, psw.spec.test

Examples

# 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" );
tmp <- psw.balance( data = test_data, weight = "MW", form.ps = form.ps,
V.name = V.name );

Mirror histogram

Description

psw.mirror.hist is used to plot the mirror histogram that visualizes the propensity score distributions in both treatment groups.

Usage

psw.mirror.hist(data, form.ps, weight, add.weight = FALSE, nclass = 50,
  K = 4)

Arguments

Details

See psw.

Value

NULL.

See Also

psw

Examples

# Load the test data set
data(test_data);
# Propensity score model
form.ps <- "Z ~ X1 + X2 + X3 + X4";
tmp <- psw.mirror.hist( data = test_data, weight = "MW", form.ps = form.ps,
add.weight = TRUE );

Propensity score model specification test

Description

psw.spec.test is used to test the sufficiency of propensity score model in balancing covariates between groups.

Usage

psw.spec.test(data, form.ps, weight, V.name, trans.type, K = 4)

Arguments

Details

In the data set, treatment indicator should be numerically specified such that a value of 1 indicates the treated and a value of 0 indicates the control. The null hypothesis is that the propensity score model is correctly specified; the alternative is that the propensity score model is misspecified. Therefore, this test is a goodness-of-fit test of propensity score model, with the test statistic being a metric of covariate balance. #' Rejection of the specification test implies current propensity score model is inadquate for balancing covariates between groups.

Value

A list of model specification test results.

See Also

psw, psw.balance

Examples

# 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" );
# A vector of transformation types for covariates in V.name.
trans.type <- c( "identity", "identity", "logit", "logit" );
tmp <- psw.spec.test( data = test_data, form.ps = form.ps,
weight = "MW", V.name = V.name, trans.type = trans.type );

Propensity score weighting estimator

Description

psw.wt is used to estimate the weighted treatment effect estimator (without double robustness).

Usage

psw.wt(data, form.ps, weight, out.var, family = "gaussian", K = 4)

Arguments

Details

psw.wt is used to estimate the weighted estimator, \hat{\Delta}, and make inference using the sandwich variance estimator that takes into account the sampling variability in the estimated propensity score.

Value

A list of weighting method, fitted propensity score model, estimated propenstity scores, estimated propensity score weights, weighted estimator and standard error estimator

See Also

psw

Examples

# Load the test data set
data(test_data);
# Propensity score model
form.ps <- "Z ~ X1 + X2 + X3 + X4";
tmp <- psw.wt( data = test_data, weight = "ATE", form.ps = form.ps,
out.var = "Y", family = "gaussian" );

Test data

Description

A simulated data frame for illustration. In the test data, X_1 and X_2 are continuous variables, X_3 and X_4 are binary variables, Y is the continuous outcome, and Z is the dichotomous treatment indicator.

Usage

test_data