ngreifer / WeightIt

WeightIt: an R package for propensity score weighting
https://ngreifer.github.io/WeightIt/
102 stars 12 forks source link

Extremely slow queries #27

Closed maellecoursonnais closed 2 years ago

maellecoursonnais commented 3 years ago

Hi,

First, thanks for the really helpful packages you created, they are very useful when doing IPTW!

I am currently doing IPTW manually and wanted to try out the WeightIt package. I have a dataset with ~50,000 observations over 15 years, and approximately 10 covariates (baseline + time-varying). Using a manually build function based on GLM estimation of the PS, I get the results in less than 3 minutes.

However, with the WeightIt package, I never actually got a chance to let the code running till the end (maximum was 3 days and 2 hours). I use the following function: weightitMSM and arguments: estimand = "ATE", method = "PS", stabilize = T. Have you already experienced this?

Thanks

ngreifer commented 3 years ago

Hi Maël,

I looked into this and was surprised at what I saw. This is indeed a problem with high numbers of time points. The cause is actually the stabilization. The stabilization factor for time point t when stabilize = TRUE is a fully saturated model for the treatment at time t given the previous treatments. A fully saturated model means that all the treatments are multiplied together in an interaction so there is no smoothing over possible combinations of treatments. With 15 treatments (i.e., 14 prior treatments), this yields a model with 2^14=16384 predictors! This model takes a long time to fit (in addition to all the models fit prior to it). The solution would be to fit a non-saturated model and make smoothing assumptions. Unfortunately, I don't allow this in WeightIt because the literature on IPW for MSMs seems to require saturated models for the stabilization factor. Even if you supply an argument to num.formula and try to use your own stabilization model formula, the fully saturated specification will be added on to it. With few treatment periods, this is desirable; with many, it is clearly not.

I will need to think about how to deal with this for WeightIt, but you also need to think about this for your analysis. Are you using a fully saturated model for your stabilization factors when you manually specify them? If not, why not, and why don't you think it's a problem? I honestly don't know what the literature recommends in this case.

To get around this in WeightIt, the solution is to fit one set of weights, call them w, using the function call as you specified except with stabilize = FALSE. Then, run the same function call again, again with stabilize = FALSE, but replacing the model formulas with the formulas for the stabilization factors (i.e., the same models but without the covariates), and call the weights s. You can then use w/s as your final weights for effect estimation. This requires you to specify your models for the stabilization factors, which, as I mentioned before, is not straightforward without a complete cross of all the treatments. One potential option would be to use machine learning methods to estimate both the propensity scores and stabilization weights, letting the algorithm decide which interactions are important. BART would be a good choice for this purpose.

I think I will edit weightitMSM() so that users can supply their own stabilization formulas. Currently, num.formula is quite limited, but it would be straightforward to expand its capabilities, which would require more manual specification by the user.

Please let me know what you decide. Thanks for bringing this up.

Noah

maellecoursonnais commented 3 years ago

Thanks for this insightful answer. To answer your first question, I have indeed a rather sparse model of the stabilization factor, with no interaction (so only the 15 treatments). I think I've never seen a paper using a fully saturated model for the stabilization factor when the number of treatments was too big; and by experience, the fully saturated model is a lot more cumbersome but does not really have added value. I have a follow-up question, are you implying in your answer that the propensity score (the denominator) in the weightitMSM() is not a fully saturated model?

Anyway, your proposition sounds quite nice actually, I will try it and let you know if it actually worked.

Maël

ngreifer commented 3 years ago

Propensity score models are typically not fully saturated. If there are any continuous covariates or if any interactions are omitted, the model will not be saturated (unless the number of predictors equals the number of units). In a fully saturated model, no parametric assumptions are made.

maellecoursonnais commented 3 years ago

Hey Noah,

I've implemented some of your comments. They are right on point. Using the two-step weight computation, the propensity scores are computed much faster. Even more, using BART gives better results in terms of covariate balancing. Thanks!