weightitMSM() allows for the easy generation of balancing
weights for marginal structural models for time-varying treatments using a
variety of available methods for binary, continuous, and multi-category
treatments, as well as censoring. Some of these methods exist in other packages, which weightit()
calls; these packages must be installed to use the desired method.
Usage
weightitMSM(
formula.list,
data = NULL,
method = "glm",
stabilize = FALSE,
by = NULL,
s.weights = NULL,
num.formula = NULL,
missing = NULL,
verbose = FALSE,
include.obj = FALSE,
keep.mparts = TRUE,
is.MSM.method,
weightit.force = FALSE,
...
)Arguments
- formula.list
a list of formulas corresponding to each time point with the time-specific treatment variable on the left hand side and pre-treatment covariates to be balanced on the right hand side. The formulas must be in temporal order, and must contain all covariates to be balanced at that time point (i.e., treatments and covariates featured in early formulas should appear in later ones). Interactions and functions of covariates are allowed. As in
weightit(), a formula may have an empty right hand side (e.g.,A_1 ~ 1), which requests a marginal model at that time point; see Empty model formulas in Details atweightit().- data
an optional data set in the form of a data frame that contains the variables in the formulas in
formula.list. This must be a wide data set with exactly one row per unit.- method
a string of length 1 containing the name of the method that will be used to estimate weights. See
weightit()for allowable options. The default is"glm", which estimates the weights using generalized linear models.- stabilize
logical; whether or not to stabilize the weights. Stabilizing the weights involves fitting a model predicting treatment at each time point from treatment status at prior time points. IfTRUE, a fully saturated model will be fit (i.e., all interactions between all treatments up to each time point), essentially using the observed treatment probabilities in the numerator (for binary and multi-category treatments). This may yield an error if some combinations are not observed. Default isFALSE. To manually specify stabilization model formulas, e.g., to specify non-saturated models, usenum.formula. With many time points, saturated models may be time-consuming or impossible to fit.- by
a string containing the name of the variable in
datafor which weighting is to be done within categories or a one-sided formula with the stratifying variable on the right-hand side. For example, ifby = "gender"orby = ~gender, a separate propensity score model or optimization will occur within each level of the variable"gender". Only onebyvariable is allowed; to stratify by multiply variables simultaneously, create a new variable that is a full cross of those variables usinginteraction().- s.weights
an optional vector of sampling weights or the name of a variable in
datathat contains sampling weights. See the individual pages for each method for information on whether sampling weights can be supplied.- num.formula
an optional one-sided formula with the stabilization factors (other than the previous treatments) on the right hand side, which adds, for each time point, the stabilization factors to a model saturated with previous treatments. See Cole & Hernán (2008) for a discussion of how to specify this model; including stabilization factors can change the estimand without proper adjustment, and should be done with caution. Can also be a list of one-sided formulas, one for each entry of
formula.list, including any censoring entries. Unless you know what you are doing, we recommend settingstabilize = TRUEand ignoringnum.formula.- missing
character; how missing data should be handled. The options and defaults depend on themethodused. Ignored if no missing data is present. It should be noted that multiple imputation outperforms all available missingness methods available inweightit()and should probably be used instead. Consider the MatchThem package for the use ofweightit()with multiply imputed data.- verbose
logical; whether to print additional information output by the fitting function. Default isFALSEto suppress output.- include.obj
logical; whether to include in the output a list of the fit objects created in the process of estimating the weights at each time point. For example, withmethod = "glm", a list of theglmobjects containing the propensity score models at each time point will be included. See the help pages for each method for information on what object will be included ifTRUE.- keep.mparts
logical; whether to include in the output components necessary to estimate standard errors that account for estimation of the weights inglm_weightit(). Default isTRUEif such parts are present. See the individual pages for each method for whether these components are produced. Set toFALSEto keep the output object smaller, e.g., if standard errors will not be computed usingglm_weightit().- is.MSM.method
logical; whether the method estimates weights for multiple time points all at once (TRUE) or by estimating weights at each time point and then multiplying them together (FALSE). This is only relevant for user-specified functions.- weightit.force
logical; several methods are not valid for estimating weights with longitudinal treatments, and will produce an error message if attempted. Set toTRUEto bypass this error message.- ...
other arguments that control aspects of fitting that are not covered by the above arguments. See Details at
weightit().
Value
A weightitMSM object with the following elements:
- weights
The estimated weights, one for each unit.
- treat.list
A list of the values of the time-varying treatment variables.
- covs.list
A list of the covariates used in the fitting at each time point. Only includes the raw covariates, which may have been altered in the fitting process.
- estimand
"ATE", currently the only estimand for MSMs with binary or multi-category treatments.
- method
The weight estimation method specified.
- s.weights
The provided sampling weights.
- by
A data.frame containing the
byvariable when specified.- stabilization
The stabilization factors, if any.
When censoring is modeled (i.e., when any entry of formula.list has its left
side wrapped in .cens()), treat.list and covs.list describe the treatment
models only, while formula.list is kept exactly as supplied, markers included,
so that update() round-trips. The following additional components describe the
censoring models:
- cens.list
A list of the values of the censoring indicators, one entry per censoring time point. Each is 0 for units still under observation and 1 for units censored at that time point, and
NAfor units censored earlier.- cens.covs.list
A list of the covariates used to fit each censoring model. As with
covs.list, only the raw covariates are included.- cens.formula.list
A list of the censoring model formulas, with the
.cens()marker retained on the left side.- cens.time
The positions of the censoring models within
formula.list, so thatformula.list[cens.time]recovers them and their timing relative to the treatment models can be determined.- at.risk
A logical matrix with one row per unit and one column per entry of
formula.list, named for the treatment or censoring variable modeled at that entry. Each column records which units were still under observation when that model was fit, i.e., which units contributed to it. Useful for assessing balance; see.cens().
Censored units have a final weight of exactly 0, so weights is 0 for any unit
censored at any time point.
When keep.mparts is TRUE (the default) and the chosen method is
compatible with M-estimation, the components related to M-estimation for use
in glm_weightit() are stored in the "Mparts.list" attribute. When by is
specified, the per-stratum components are combined into "Mparts.list" so
that the standard errors produced by glm_weightit() are asymptotically
equivalent to those from estimating the weights from models in which the by
variable is fully interacted with all the covariates at every time point. (For
method = "cbps", this requires is.MSM.method = FALSE, i.e., estimating a
separate model at each time point, as M-estimation is not supported for the
single-model MSM version of CBPS.)
Details
In general, weightitMSM() works by separating the estimation of weights
into separate procedures for each time period based on the formulas provided.
For each formula, weightitMSM() simply applies weightit() to that formula,
collects the weights for each time period, and multiplies them together to
arrive at longitudinal balancing weights.
Each formula should contain all the covariates to be balanced on. For
example, the formula corresponding to the second time period should contain
all the baseline covariates, the treatment variable at the first time period,
and the time-varying covariates that took on values after the first treatment
and before the second. Currently only "wide" data sets are supported, where each
unit is represented by exactly one row that contains its covariate and
treatment history encoded in separate variables. You can use reshape() or
other functions to transform your data into this format; see example below.
Censoring weights (IPCW)
Censoring can be modeled by including entries in formula.list whose left side
is wrapped in .cens(), placed in temporal order among the treatment
models. For example,
weightitMSM(list(A_1 ~ X1_0 + X2_0,
A_2 ~ X1_1 + X2_1 + A_1,
.cens(C_2) ~ X1_1 + X2_1 + A_1 + A_2,
A_3 ~ X1_2 + X2_2 + A_2),
data = d, method = "glm")models censoring occurring after the second treatment. Each censoring indicator
must be 0 for units still under observation and 1 for units censored at that time
point. See .cens() for details of what the resulting weights
estimate.
Every model, treatment or censoring, is fit only among the units still under
observation when it is reached, and the resulting weights are multiplied together
across time points as usual. A unit censored at any time point therefore has a
final weight of exactly 0. Because such units drop out, missing values are
permitted in the treatments and covariates that follow their censoring; missing
values among units still under observation remain an error. at.risk in the
output has one column per time point recording which units were under
observation when that model was fit.
Censoring time points are stabilized in exactly the same way as treatment time
points: with stabilize = TRUE, the numerator of a censoring weight is a model
for that censoring indicator given the preceding treatments (or a marginal model
when no treatment precedes it), and num.formula adds stabilization factors to
it as it does for a treatment. When num.formula is supplied as a list, it must
have one entry per entry of formula.list, censoring entries included. The
numerator of a censoring weight is itself a censoring model, so the stabilized
weight is \(P(C = 0 | \cdot) / P(C = 0 | X)\) for the units still under
observation and remains exactly 0 for those censored.
The right side of a censoring formula may be empty, as in .cens(C_2) ~ 1, which
requests a marginal censoring model that assumes censoring at that time point is
independent of the covariates; its contribution to the product is \(1/P(C = 0)\)
for the units still under observation and 0 for those censored there. Everything
else is unaffected: the risk sets, the missing values permitted after censoring,
stabilize, by, and M-estimation all work as they do for a
covariate-dependent censoring model, and empty and non-empty censoring formulas
can be mixed freely. Each time point is fit separately, so only the empty ones
take the intercept-only shortcut described in Empty model formulas in Details at
weightit(); when is.MSM.method = TRUE there is no shortcut to take, because a
single set of weights is estimated for all time points at once, and a time point
with no covariates instead contributes only its intercept balance condition.
References
Cole, S. R., & Hernán, M. A. (2008). Constructing Inverse Probability Weights for Marginal Structural Models. American Journal of Epidemiology, 168(6), 656–664. doi:10.1093/aje/kwn164
See also
weightit() for information on the allowable methods
summary.weightitMSM() for summarizing the weights
Examples
data("msmdata")
(W1 <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
A_2 ~ X1_1 + X2_1 +
A_1 + X1_0 + X2_0,
A_3 ~ X1_2 + X2_2 +
A_2 + X1_1 + X2_1 +
A_1 + X1_0 + X2_0),
data = msmdata,
method = "glm"))
#> A weightitMSM object
#> - method: "glm" (propensity score weighting with GLM)
#> - number of obs.: 7500
#> - sampling weights: none
#> - number of time points: 3 (A_1, A_2, A_3)
#> - treatment:
#> + time 1: 2-category
#> + time 2: 2-category
#> + time 3: 2-category
#> - covariates:
#> + baseline: X1_0, X2_0
#> + after time 1: X1_1, X2_1, A_1, X1_0, X2_0
#> + after time 2: X1_2, X2_2, A_2, X1_1, X2_1, A_1, X1_0, X2_0
summary(W1)
#> Time 1
#> - Weight ranges:
#>
#> Min Max
#> treated 1.079 |---------------------------| 403.483
#> control 1.276 |-------------------| 284.764
#>
#> - Units with the 5 most extreme weights by group:
#>
#> 5488 3440 3593 1286 5685
#> treated 166.992 170.555 196.414 213.193 403.483
#> 2594 2932 5226 1875 2533
#> control 155.625 168.964 172.42 245.882 284.764
#>
#> - Weight statistics:
#>
#> Coef of Var MAD Entropy # Zeros
#> treated 1.914 0.816 0.649 0
#> control 1.706 0.862 0.67 0
#>
#> - Effective Sample Sizes:
#>
#> Control Treated
#> Unweighted 3306. 4194.
#> Weighted 845.79 899.4
#>
#> Time 2
#> - Weight ranges:
#>
#> Min Max
#> treated 1.079 |---------------------------| 403.483
#> control 1.276 |----------------| 245.882
#>
#> - Units with the 5 most extreme weights by group:
#>
#> 2932 3440 3593 2533 5685
#> treated 168.964 170.555 196.414 284.764 403.483
#> 2594 5488 5226 1286 1875
#> control 155.625 166.992 172.42 213.193 245.882
#>
#> - Weight statistics:
#>
#> Coef of Var MAD Entropy # Zeros
#> treated 1.892 0.819 0.652 0
#> control 1.748 0.869 0.686 0
#>
#> - Effective Sample Sizes:
#>
#> Control Treated
#> Unweighted 3701. 3799.
#> Weighted 912.87 829.87
#>
#> Time 3
#> - Weight ranges:
#>
#> Min Max
#> treated 1.079 |---------------------------| 403.483
#> control 1.276 |---------| 148.155
#>
#> - Units with the 5 most extreme weights by group:
#>
#> 3593 1286 1875 2533 5685
#> treated 196.414 213.193 245.882 284.764 403.483
#> 6862 168 3729 6158 3774
#> control 88.072 97.827 104.623 121.845 148.155
#>
#> - Weight statistics:
#>
#> Coef of Var MAD Entropy # Zeros
#> treated 1.832 0.975 0.785 0
#> control 1.254 0.683 0.412 0
#>
#> - Effective Sample Sizes:
#>
#> Control Treated
#> Unweighted 4886. 2614.
#> Weighted 1900.26 600.12
#>
cobalt::bal.tab(W1)
#> Balance summary across all time points
#> Times Type Max.Diff.Adj
#> X1_0 1, 2, 3 Contin. 0.0342
#> X2_0 1, 2, 3 Binary 0.0299
#> X1_1 2, 3 Contin. 0.0657
#> X2_1 2, 3 Binary 0.0299
#> A_1 2, 3 Binary 0.0262
#> X1_2 3 Contin. 0.0643
#> X2_2 3 Binary 0.0096
#> A_2 3 Binary 0.0054
#>
#> Effective sample sizes
#> - Time 1
#> Control Treated
#> Unadjusted 3306. 4194.
#> Adjusted 845.79 899.4
#> - Time 2
#> Control Treated
#> Unadjusted 3701. 3799.
#> Adjusted 912.87 829.87
#> - Time 3
#> Control Treated
#> Unadjusted 4886. 2614.
#> Adjusted 1900.26 600.12
# Using stabilization factors
W2 <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
A_2 ~ X1_1 + X2_1 +
A_1 + X1_0 + X2_0,
A_3 ~ X1_2 + X2_2 +
A_2 + X1_1 + X2_1 +
A_1 + X1_0 + X2_0),
data = msmdata,
method = "glm",
stabilize = TRUE,
num.formula = list(~ 1,
~ A_1,
~ A_1 + A_2))
# Same as above but with fully saturated stabilization factors
# (i.e., making the last entry in 'num.formula' A_1*A_2)
W3 <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
A_2 ~ X1_1 + X2_1 +
A_1 + X1_0 + X2_0,
A_3 ~ X1_2 + X2_2 +
A_2 + X1_1 + X2_1 +
A_1 + X1_0 + X2_0),
data = msmdata,
method = "glm",
stabilize = TRUE)