method_glm | R Documentation |
This page explains the details of estimating weights from generalized linear model-based propensity scores by setting method = "glm"
in the call to weightit()
or weightitMSM()
. This method can be used with binary, multi-category, 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 multi-category 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.
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.
For multi-category treatments, the propensity scores are estimated using multinomial regression from one of a few functions depending on the argument supplied to multi.method
(see Additional Arguments below). 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. Ordinal treatments are treated exactly the same as non-order multi-category treatments except that additional models are available to estimate the generalized propensity score (e.g., ordinal logistic regression).
For continuous treatments, weights are estimated as w_i = f_A(a_i) / f_{A|X}(a_i)
, where f_A(a_i)
(known as the stabilization factor) is the unconditional density of treatment evaluated the observed treatment value and f_{A|X}(a_i)
(known as the generalized propensity score) is the conditional density of treatment given the covariates evaluated at the observed value of treatment. The shape of f_A(.)
and f_{A|X}(.)
is controlled by the density
argument described below (normal distributions by default), and the predicted values used for the mean of the conditional density are estimated using linear regression. Kernel density estimation can be used instead of assuming a specific density for the numerator and denominator by setting density = "kernel"
. Other arguments to density()
can be specified to refine the density estimation parameters.
For longitudinal treatments, the weights are the product of the weights estimated at each time point.
Sampling weights are supported through s.weights
in all scenarios except for multi-category 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.
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 the covariate medians (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 NA
s.
"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.
For binary treatments, M-estimation is supported when link
is neither "flic"
nor "flac"
(see below). For multi-category treatments, M-estimation is supported when multi.method
is "weightit"
(the default) or "glm"
. For continuous treatments, M-estimation is supported when density
is not "kernel"
. The conditional treatment variance and unconditional treatment mean and variance are included as parameters to estimate, as these all go into calculation of the weights. For all treatment types, M-estimation is not supported when missing = "saem"
. See glm_weightit()
and vignette("estimating-effects")
for details. For longitudinal treatments, M-estimation is supported whenever the underlying methods are.
For binary treatments, the following additional argument can be specified:
link
the link used in the generalized linear model for the propensity scores. 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. link
can also be either "flic"
or "flac"
to fit the corresponding Firth corrected logistic regression models implemented in the logistf package.
For multi-category treatments, the following additional arguments can be specified:
multi.method
the method used to estimate the generalized propensity scores. Allowable options include "weightit"
(the default) to use multinomial logistic regression implemented in WeightIt, "glm"
to use a series of binomial models using glm()
, "mclogit"
to use multinomial logistic regression as implemented in \pkgfunmclogitmblogit, "mnp"
to use Bayesian multinomial probit regression as implemented in \pkgfunMNPMNP, and "brmultinom"
to use bias-reduced multinomial logistic regression as implemented in \pkgfunbrglm2brmultinom. "weightit"
and "mclogit"
should give near-identical results, the main difference being increased robustness and customizability when using "mclogit"
at the expense of not being able to use M-estimation to compute standard errors after weighting. For ordered treatments, allowable options include "weightit"
(the default) to use ordinal regression implemented in WeightIt or "polr"
to use ordinal regression implemented in \pkgfunMASSpolr, unless link
is "br.logit"
, in which case bias-reduce ordinal logistic regression as implemented in \pkgfunbrglm2bracl is used. Ignored when missing = "saem"
. Using the defaults allows for the use of M-estimation and requires no additional dependencies, but other packages may provide benefits such as speed and flexibility.
link
The link used in the multinomial, binomial, or ordered regression model for the generalized propensity scores depending on the argument supplied to multi.method
. When multi.method = "glm"
, link
can be any of those allowed by binomial()
. When treatment is ordered and multi.method
is "weightit"
or "polr"
, link
can be any of those allowed by MASS::polr()
or "br.logit"
. Otherwise, link
should be "logit"
or not specified.
For continuous treatments, the following additional 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).
Can also be "kernel"
to use kernel density estimation, which calls density()
to estimate the numerator and denominator densities for the weights. (This used to be requested by setting use.kernel = TRUE
, which is now deprecated.)
bw
, adjust
, kernel
, n
If density = "kernel"
, the arguments to density()
. The defaults are the same as those in density()
except that n
is 10 times the number of units in the sample.
plot
If density = "kernel"
, whether to plot the estimated densities.
link
The link used to fit the linear model for the generalized propensity score. Can be any allowed by gaussian()
.
Additional arguments to glm()
can be specified as well when it is used for fitting. 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 \pkgfunmisaemmiss.glm and \pkgfunmisaempredict.miss.glm, except the method
argument in \pkgfunmisaempredict.miss.glm is replaced with saem.method
.
For continuous treatments in the presence of missing data with missing = "saem"
, additional arguments are passed to \pkgfunmisaemmiss.lm and \pkgfunmisaempredict.miss.lm.
obj
When include.obj = TRUE
, the (generalized) propensity score model fit. For binary treatments, the output of the call to glm()
or the requested fitting function. For multi-category treatments, the output of the call to the fitting function (or a list thereof if multi.method = "glm"
). For continuous treatments, the output of the call to glm()
for the predicted values in the denominator density.
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. \Sexpr[results=rd]{tools:::Rd_expr_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). \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3102/1076998609359785")}
Bias-reduced logistic regression
See references for the brglm2 package.
Firth corrected logistic regression
Puhr, R., Heinze, G., Nold, M., Lusa, L., & Geroldinger, A. (2017). Firth’s logistic regression with rare events: Accurate effect estimates and predictions? Statistics in Medicine, 36(14), 2302–2317. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1002/sim.7273")}
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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.csda.2019.106907")}
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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1037/a0024918")}
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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1097/EDE.0000000000000053")}
SAEM linear 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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.csda.2019.106907")}
weightit()
, weightitMSM()
, get_w_from_ps()
library("cobalt")
data("lalonde", package = "cobalt")
#Balancing covariates between treatment groups (binary)
(W1 <- weightit(treat ~ age + educ + married +
nodegree + re74, data = lalonde,
method = "glm", estimand = "ATT",
link = "probit"))
summary(W1)
bal.tab(W1)
#Balancing covariates with respect to race (multi-category)
(W2 <- weightit(race ~ age + educ + married +
nodegree + re74, data = lalonde,
method = "glm", estimand = "ATE"))
summary(W2)
bal.tab(W2)
#Balancing covariates with respect to re75 (continuous)
#with kernel density estimate
(W3 <- weightit(re75 ~ age + educ + married +
nodegree + re74, data = lalonde,
method = "glm", density = "kernel"))
summary(W3)
bal.tab(W3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.