Subclassification on a propensity and prognostic score grid with k × k subclasses, using R MatchIt

Question

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! :)

Answer 1

Unfortunately, matchit() with method = "subclass" cannot be used for this purpose, but there is a workaround. You should note that the formula you supplied has a few problems; the response variable is supposed to be the treatment, not the propensity or prognostic score. Matching variables will always go on the right-hand side, as they will in the workaround below.

The workaround is to consider subclassification as a version of coarsened exact matching. That is, you coarsen the variables into categories and create strata based on the combinations of coarsened variables. Here, the two variables of interest are the propensity score and prognostic score. I'll provide a little demo below:

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

#Estimate propensity score (using probit regression here)
m0 <- matchit(treat ~ age + educ + race + married + re74 + re75,
              data = lalonde, method = NULL, distance = "glm",
              link = "probit")
ps <- m0$distance

#Estimate prognostic score (using OLS in control group)
fit0 <- lm(re78 ~ age + educ + race + married + re74 + re75,
           data = lalonde, subset = treat == 0)
prog <- predict(fit0, newdata = lalonde)

#Do CEM on ps and prog
m.out <- matchit(treat ~ ps + prog, data = lalonde,
                 method = "cem", cutpoints = list(ps = "q5", prog = "q5"),
                 estimand = "ATT")

#Assess balance (output omitted for space)
summary(m.out, data = lalonde,
        addlvariables = ~ age + educ + race + married + re74 + re75)

In the matchit() call with method = "cem", you can change the number of subclasses of each variable by changing the arguments supplied to cutpoints. For example, to request 10 bins for the propensity score, you can set ps = "q10". These use the quantiles of the covariate in the full sample, unlike with method = "subclass", which uses the quantiles of the propensity score in the target sample (e.g., the treated group for the ATT). Also unlike with method = "subclass", using method = "cem" will drop any units that are in a subclass containing only treated units or only control units.

To estimate the treatment effect, you need to use bootstrapping around the whole process to ensure that the uncertainty due to estimating the prognostic score and propensity score is taken into account. To estimate the effect in each bootstrap sample, you can use the estimated matching weights. You don't need to worry about a standard error because that is computed using the bootstrap.

#Use bootstrapping to estimate effect and CI

boot_fun <- function(d, i) {
  #Estimate propensity score (using probit regression here)
  m0 <- matchit(treat ~ age + educ + race + married + re74 + re75,
                data = d[i,], method = NULL, distance = "glm",
                link = "probit")
  ps <- m0$distance
  
  #Estimate prognostic score (using OLS in control group)
  fit0 <- lm(re78 ~ age + educ + race + married + re74 + re75,
             data = d[i,], subset = treat == 0)
  prog <- predict(fit0, newdata = d[i,])
  
  #Do CEM on ps and prog
  m.out <- matchit(treat ~ ps + prog, data = d[i,],
                   method = "cem", cutpoints = list(ps = "q5", prog = "q5"),
                   estimand = "ATT")
  
  #Estimate treatment effect
  md <- match.data(m.out, data = d[i,])
  
  coef(lm(re78 ~ treat, data = md, weights = weights))["treat"]
}

set.seed(2022)
(b <- boot::boot(lalonde, boot_fun, R = 999))
#> Bootstrap Statistics :
#>     original    bias    std. error
#> t1* 1904.958 -341.4798    928.2113
boot::boot.ci(b, type = "bca")
#> Intervals : 
#> Level       BCa          
#> 95%   ( 402, 4269 )  

(I hid some of the output for space.) In this example, our effect estimate was 1905 with a confidence interval of (402, 4269). It's always better to use more bootstrap replications (e.g., ~10000). Make sure you set a seed as I did above for purposes of replicability.