I would like to perform a joint subclassification of some data on the propensity and prognostic scores as described in this paper "On the joint use of propensity and prognostic scores in estimation of the average treatment effect on the treated: a simulation study":
3.4.3 Subclassification on a propensity and prognostic score grid with k × k subclasses
Subclassification on a propensity and prognostic score grid with k × k subclasses is proposed here as an alternative means of combining the propensity and prognostic scores. It is inspired by the propensity function of Imai and van Dyk 42. Imai and van Dyk use this method to estimate the effects of smoking on medical expenditure 42, where they treat the duration and frequency of smoking as a bivariate treatment, estimating a separate propensity score for each.
The propensity and prognostic scores are first estimated. The data are then divided into an n × m grid of subclasses based on the quantiles of the estimated propensity and prognostic scores for pre-specified n and m values. In practice, these values could be chosen in an iterative process where the number of subclasses is varied until an acceptable level of propensity and/or prognostic score balance on each covariate is attained (while maintaining sufficient sample sizes in each subclass). Propensity score balance can be checked by comparing the distribution of each covariate across treatment groups. Prognostic score balance can be checked by testing for association between the outcome and each covariate in the control group along subclasses of the prognostic score. The ATT is calculated as a weighted average of the within-subclass estimates. Previous simulation work examining the performance of subclassification on the propensity function has suggested that for sufficiently large sample sizes increasing the total number of subclasses is associated with greater reductions in bias 42. Preliminary simulation work, not reported here, suggests that this observation also holds for subclassification on a propensity and prognostic score grid with k × k subclasses.
I would like to use the R package
MatchIt for this using
method = "subclass", but it seems like
MatchIt can only perform subclassification on a single response variable, and it complains if I pass in a formula like
cbind(prop, prog) ~ X + Y + Z (which I'm not even sure is correct anyway!).
Can someone tell me if it is possible to do subclassification over a grid with
MatchIt, or another package? Otherwise, can someone tell me how I should go about computing the bounds for my different strata? I want to replicate the same 5x5 grid.
Thank you! :)