Propensity Score Weighting Using Generalized Linear Models

This page explains the details of estimating weights from generalized linear model-based propensity scores by setting method = "ps" in the call to weightit() or weightitMSM(). This method can be used with binary, multinomial, and continuous treatments.

In general, this method relies on estimating propensity scores with a parametric generalized linear model and then converting those propensity scores into weights using a formula that depends on the desired estimand. For binary and multinomial treatments, a binomial or multinomial regression model is used to estimate the propensity scores as the predicted probability of being in each treatment given the covariates. For ordinal treatments, an ordinal regression model is used to estimate generalized propensity scores. For continuous treatments, a generalized linear model is used to estimate generalized propensity scores as the conditional density of treatment given the covariates.

Binary Treatments

For binary treatments, this method estimates the propensity scores using glm(). An additional argument is link, which uses the same options as link in family(). The default link is "logit", but others, including "probit", are allowed. The following estimands are allowed: ATE, ATT, ATC, ATO, ATM, and ATOS. Weights can also be computed using marginal mean weighting through stratification for the ATE, ATT, and ATC. See get_w_from_ps() for details.

Multinomial Treatments

For multinomial treatments, the propensity scores are estimated using multinomial regression from one of a few functions depending on the requested link: for logit ("logit") and probit ("probit") links, mlogit::mlogit() from the mlogit package is used; for the Bayesian probit ("bayes.probit") link, MNP::mnp() from the MNP package is used; and for the biased-reduced multinomial logistic regression ("br.logit"), brglm2::brmultinom() from the brglm2 package is used. If the treatment variable is an ordered factor, MASS::polr() from the MASS package is used to fit ordinal regression unless link = "br.logit", in which case brglm2::bracl() from brglm2 is used. Any of the methods allowed in the method argument of polr() can be supplied to link. The following estimands are allowed: ATE, ATT, ATC, ATO, and ATM. The weights for each estimand are computed using the standard formulas or those mentioned above. Weights can also be computed using marginal mean weighting through stratification for the ATE, ATT, and ATC. See get_w_from_ps() for details.

Continuous Treatments

For continuous treatments, the generalized propensity score is estimated using linear regression. The conditional density can be specified as normal or another distribution. In addition, kernel density estimation can be used instead of assuming a specific density for the numerator and denominator of the generalized propensity score by setting use.kernel = TRUE. Other arguments to density() can be specified to refine the density estimation parameters. plot = TRUE can be specified to plot the density for the numerator and denominator, which can be helpful in diagnosing extreme weights.

Longitudinal Treatments

For longitudinal treatments, the weights are the product of the weights estimated at each time point.

Sampling Weights

Sampling weights are supported through s.weights in all scenarios except for multinomial treatments with link = "bayes.probit" and for binary and continuous treatments with missing = "saem" (see below). Warning messages may appear otherwise about non-integer successes, and these can be ignored.

Missing Data

In the presence of missing data, the following value(s) for missing are allowed:

"ind" (default): First, for each variable with missingness, a new missingness indicator variable is created which takes the value 1 if the original covariate is NA and 0 otherwise. The missingness indicators are added to the model formula as main effects. The missing values in the covariates are then replaced with 0s (this value is arbitrary and does not affect estimation). The weight estimation then proceeds with this new formula and set of covariates. The covariates output in the resulting weightit object will be the original covariates with the NAs.
"saem": For binary treatments with link = "logit" or continuous treatments, a stochastic approximation version of the EM algorithm (SAEM) is used via the misaem package. No additional covariates are created. See Jiang et al. (2019) for information on this method. In some cases, this is a suitable alternative to multiple imputation.

Additional Arguments

The following additional arguments can be specified:

link: The link used in the generalized linear model for the propensity scores. For binary treatments, link can be any of those allowed by binomial(). A br. prefix can be added (e.g., "br.logit"); this changes the fitting method to the bias-corrected generalized linear models implemented in the brglm2 package. For multicategory treatments, link can be "logit", "probit", "bayes.probit", or "br.logit". For ordered treatments, link can be any of those allowed by the method argument of MASS::polr() or "br.logit". For continuous treatments, link can be any of those allowed by gaussian().

For continuous treatments only, the following arguments may be supplied:

density: A function corresponding the conditional density of the treatment. The standardized residuals of the treatment model will be fed through this function to produce the numerator and denominator of the generalized propensity score weights. If blank, dnorm() is used as recommended by Robins et al. (2000). This can also be supplied as a string containing the name of the function to be called. If the string contains underscores, the call will be split by the underscores and the latter splits will be supplied as arguments to the second argument and beyond. For example, if density = "dt_2" is specified, the density used will be that of a t-distribution with 2 degrees of freedom. Using a t-distribution can be useful when extreme outcome values are observed (Naimi et al., 2014). Ignored if use.kernel = TRUE (described below).
use.kernel: If TRUE, uses kernel density estimation through the density() function to estimate the numerator and denominator densities for the weights. If FALSE, the argument to the density parameter is used instead.
bw, adjust, kernel, n: If use.kernel = TRUE, the arguments to the density() function. The defaults are the same as those in density except that n is 10 times the number of units in the sample.
plot: If use.kernel = TRUE with continuous treatments, whether to plot the estimated density.

For binary treatments, additional arguments to glm() can be specified as well. The method argument in glm() is renamed to glm.method. This can be used to supply alternative fitting functions, such as those implemented in the glm2 package. Other arguments to weightit() are passed to ... in glm(). In the presence of missing data with link = "logit" and missing = "saem", additional arguments are passed to miss.glm and predict.miss.glm, except the method argument in predict.miss.glm is replaced with saem.method.

For multi-category treatments with link = "logit" or "probit", the default is to use multinomial logistic or probit regression using the mlogit package. To request that separate binary logistic or probit regressions are run instead, set use.mlogit = FALSE. This can be helpful when mlogit is slow or fails to converge. With link = "logit", the option use.mclogit = TRUE can be specified to request that mclogit::mblogit() from the mclogit package is used instead, which can be faster and is recommended.

For continuous treatments in the presence of missing data with missing = "saem", additional arguments are passed to miss.lm and predict.miss.lm.

Additional Outputs

obj: When include.obj = TRUE, the (generalized) propensity score model fit. For binary treatments, the output of the call to glm(). For ordinal treatments, the output of the call to MASS::polr(). For multinomial treatments with link = "logit" or "probit" and use.mlogit = TRUE, the output of the call to mlogit::mlogit(). For multinomial treatments with use.mlogit = FALSE, a list of the glm() fits. For multinomial treatments with link = "br.logit", the output of the call to brglm2::brmultinom(). For multinomial treatments with link = "bayes.probit", the output of the call to MNP::mnp(). For continuous treatments, the output of the call to glm() for the predicted values in the denominator density.

References

Binary treatments

- estimand = "ATO"

Li, F., Morgan, K. L., & Zaslavsky, A. M. (2018). Balancing covariates via propensity score weighting. Journal of the American Statistical Association, 113(521), 390–400. doi:10.1080/01621459.2016.1260466

- estimand = "ATM"

Li, L., & Greene, T. (2013). A Weighting Analogue to Pair Matching in Propensity Score Analysis. The International Journal of Biostatistics, 9(2). doi:10.1515/ijb-2012-0030

- estimand = "ATOS"

Crump, R. K., Hotz, V. J., Imbens, G. W., & Mitnik, O. A. (2009). Dealing with limited overlap in estimation of average treatment effects. Biometrika, 96(1), 187–199. doi:10.1093/biomet/asn055

- Other estimands

Austin, P. C. (2011). An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies. Multivariate Behavioral Research, 46(3), 399–424. doi:10.1080/00273171.2011.568786

- Marginal mean weighting through stratification

Hong, G. (2010). Marginal mean weighting through stratification: Adjustment for selection bias in multilevel data. Journal of Educational and Behavioral Statistics, 35(5), 499–531. doi:10.3102/1076998609359785

- Bias-reduced logistic regression

See references for the brglm2 package.

- SAEM logistic regression for missing data

Jiang, W., Josse, J., & Lavielle, M. (2019). Logistic regression with missing covariates — Parameter estimation, model selection and prediction within a joint-modeling framework. Computational Statistics & Data Analysis, 106907. doi:10.1016/j.csda.2019.106907

Multinomial Treatments

- estimand = "ATO"

Li, F., & Li, F. (2019). Propensity score weighting for causal inference with multiple treatments. The Annals of Applied Statistics, 13(4), 2389–2415. doi:10.1214/19-AOAS1282

- estimand = "ATM"

Yoshida, K., Hernández-Díaz, S., Solomon, D. H., Jackson, J. W., Gagne, J. J., Glynn, R. J., & Franklin, J. M. (2017). Matching weights to simultaneously compare three treatment groups: Comparison to three-way matching. Epidemiology (Cambridge, Mass.), 28(3), 387–395. doi:10.1097/EDE.0000000000000627

- Other estimands

McCaffrey, D. F., Griffin, B. A., Almirall, D., Slaughter, M. E., Ramchand, R., & Burgette, L. F. (2013). A Tutorial on Propensity Score Estimation for Multiple Treatments Using Generalized Boosted Models. Statistics in Medicine, 32(19), 3388–3414. doi:10.1002/sim.5753

- Marginal mean weighting through stratification

Hong, G. (2012). Marginal mean weighting through stratification: A generalized method for evaluating multivalued and multiple treatments with nonexperimental data. Psychological Methods, 17(1), 44–60. doi:10.1037/a0024918

Continuous treatments

Robins, J. M., Hernán, M. Á., & Brumback, B. (2000). Marginal Structural Models and Causal Inference in Epidemiology. Epidemiology, 11(5), 550–560.

- Using non-normal conditional densities

Naimi, A. I., Moodie, E. E. M., Auger, N., & Kaufman, J. S. (2014). Constructing Inverse Probability Weights for Continuous Exposures: A Comparison of Methods. Epidemiology, 25(2), 292–299. doi:10.1097/EDE.0000000000000053

- SAEM linear regression for missing data

Examples

library("cobalt")
data("lalonde", package = "cobalt")

#Balancing covariates between treatment groups (binary)
(W1 <- weightit(treat ~ age + educ + married +
                  nodegree + re74, data = lalonde,
                method = "ps", estimand = "ATT",
                link = "probit"))
#> A weightit object
#>  - method: "ps" (propensity score weighting)
#>  - number of obs.: 614
#>  - sampling weights: none
#>  - treatment: 2-category
#>  - estimand: ATT (focal: 1)
#>  - covariates: age, educ, married, nodegree, re74
summary(W1)
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>            Min                                  Max
#> treated 1.0000                 ||            1.0000
#> control 0.0176 |---------------------------| 1.8338
#> 
#> - Units with the 5 most extreme weights by group:
#>                                         
#>              10     9     8     4      3
#>  treated      1     1     1     1      1
#>             612   595   269   409    296
#>  control 1.2781 1.351 1.412 1.518 1.8338
#> 
#> - Weight statistics:
#> 
#>         Coef of Var   MAD Entropy # Zeros
#> treated       0.000 0.000  -0.000       0
#> control       0.804 0.691   0.322       0
#> 
#> - Effective Sample Sizes:
#> 
#>            Control Treated
#> Unweighted  429.       185
#> Weighted    260.83     185
bal.tab(W1)
#> Call
#>  weightit(formula = treat ~ age + educ + married + nodegree + 
#>     re74, data = lalonde, method = "ps", estimand = "ATT", link = "probit")
#> 
#> Balance Measures
#>                Type Diff.Adj
#> prop.score Distance   0.0252
#> age         Contin.   0.0716
#> educ        Contin.  -0.0565
#> married      Binary   0.0058
#> nodegree     Binary   0.0128
#> re74        Contin.  -0.0507
#> 
#> Effective sample sizes
#>            Control Treated
#> Unadjusted  429.       185
#> Adjusted    260.83     185

#Balancing covariates with respect to race (multinomial)
(W2 <- weightit(race ~ age + educ + married +
                  nodegree + re74, data = lalonde,
                method = "ps", estimand = "ATE",
                use.mlogit = FALSE))
#> A weightit object
#>  - method: "ps" (propensity score weighting)
#>  - number of obs.: 614
#>  - sampling weights: none
#>  - treatment: 3-category (black, hispan, white)
#>  - estimand: ATE
#>  - covariates: age, educ, married, nodegree, re74
summary(W2)
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>           Min                                   Max
#> black  1.4486 |---------------------------| 29.1514
#> hispan 1.8148  |------------------------|   26.9274
#> white  1.1180 |-|                            4.1126
#> 
#> - Units with the 5 most extreme weights by group:
#>                                                
#>             226     231     485     181     182
#>   black  7.6815  7.8028   8.566 12.2152 29.1514
#>             346     392     345     269     371
#>  hispan 17.5087 18.3157 18.4663 21.5593 26.9274
#>             432      68     404     599     531
#>   white  3.8059  3.8309  3.8835  4.0487  4.1126
#> 
#> - Weight statistics:
#> 
#>        Coef of Var   MAD Entropy # Zeros
#> black        0.861 0.426   0.188       0
#> hispan       0.546 0.407   0.134       0
#> white        0.385 0.313   0.069       0
#> 
#> - Effective Sample Sizes:
#> 
#>             black hispan  white
#> Unweighted 243.    72.   299.  
#> Weighted   139.74  55.66 260.44
bal.tab(W2)
#> Call
#>  weightit(formula = race ~ age + educ + married + nodegree + re74, 
#>     data = lalonde, method = "ps", estimand = "ATE", use.mlogit = FALSE)
#> 
#> Balance summary across all treatment pairs
#>             Type Max.Diff.Adj
#> age      Contin.       0.0374
#> educ     Contin.       0.1187
#> married   Binary       0.0540
#> nodegree  Binary       0.0510
#> re74     Contin.       0.1682
#> 
#> Effective sample sizes
#>             black hispan  white
#> Unadjusted 243.    72.   299.  
#> Adjusted   139.74  55.66 260.44

#Balancing covariates with respect to re75 (continuous)
(W3 <- weightit(re75 ~ age + educ + married +
                  nodegree + re74, data = lalonde,
                method = "ps", use.kernel = TRUE))
#> A weightit object
#>  - method: "ps" (propensity score weighting)
#>  - number of obs.: 614
#>  - sampling weights: none
#>  - treatment: continuous
#>  - covariates: age, educ, married, nodegree, re74
summary(W3)
#>                  Summary of weights
#> 
#> - Weight ranges:
#> 
#>       Min                                   Max
#> all 0.088 |---------------------------| 83.1017
#> 
#> - Units with the 5 most extreme weights:
#>                                           
#>         482     481    484     483     485
#>  all 56.933 56.9865 70.237 77.1173 83.1017
#> 
#> - Weight statistics:
#> 
#>     Coef of Var   MAD Entropy # Zeros
#> all       2.908 1.049   1.156       0
#> 
#> - Effective Sample Sizes:
#> 
#>             Total
#> Unweighted 614.  
#> Weighted    65.03
bal.tab(W3)
#> Call
#>  weightit(formula = re75 ~ age + educ + married + nodegree + re74, 
#>     data = lalonde, method = "ps", use.kernel = TRUE)
#> 
#> Balance Measures
#>             Type Corr.Adj
#> age      Contin.  -0.0501
#> educ     Contin.   0.0016
#> married   Binary  -0.0427
#> nodegree  Binary  -0.0196
#> re74     Contin.  -0.0773
#> 
#> Effective sample sizes
#>             Total
#> Unadjusted 614.  
#> Adjusted    65.03