Subgroup analyses after propensity score matching

Question

Assume an already propensity score matched cohort.
One of the variables included in the propensity score model is sex.
The main cox model investigates the comparative effect of two drugs for cancer remission.
This is the main effect I was looking for:

coxph(Surv(time, status) ~ drug, data = data)

Now I want to investigate a possible moderation (by sex).

OPTION 1:
I could use the original cohort from above investigating a possible moderation:

coxph(Surv(time, status) ~ drug:sex, data = data)

Thus, I have two models: One for the main effect estimate and one for the moderation.

OPTION 2:
According to Green & Stuart (2014), estimating separate propensity score models could show the best balance.
Thus, I created two datasets (for men and women) and calculated the propensity score for matching separately.
Next, I combined both matched datasets and calculated:

coxph(Surv(time, status) ~ drug:sex, data = data)

As I used 1:1 matching using a caliper, the original data and the combined stratified approach do have different sample sizes due to the fact that not everyone got a matched reference as well as different matched reference persons.

Questions:

1. For OPTION 1, is it legit to calculate one model for the main effect drug and one model for possible direct moderation using drug:sex? Or do I need to include both in one model as drug*sex?

2. For OPTION 2, do I need to calculate the main effect in the combined sample? Or do I use the combined data just to investigate the moderation and use my original dataset for the main effect?
   Assume I want to investigate another moderation by age categories (<= 40 and > 40 years). Thus, I would need to calculate separate propensity score models, match them and combine both separate datasets. But the main effect will of course be different between the combined data that has previously been stratified by sex and the one that has been stratified by age category.

3. In general, how do I derive separate hazard ratios and CIs for both, men and women?

Answer 1

For moderation analysis, you need to ensure you have achieved balance within each subgroup of the moderator. The most straightforward way to do so is to match separately within each level of the moderator. Another option is to use a matching procedure that ensures balance within each subgroup, such as by exact matching on subgroup. After matching, you need to extract the matched datasets and fit a Cox model with a treatment-by-moderator interaction. The coefficient on treatment gives you the log hazard ratio in the reference group and the coefficient on the interaction term gives you the difference between the log hazard ratio in the non-reference group and the log hazard ratio in the reference group. The p-value on the interaction tests whether the subgroup log hazard ratios differ. To get the log hazard ratio in the non-reference group, switch which group is the reference group and re-fit the model.

Below I'll show you how to use MatchIt to do matching for subgroup balance either by doing the match separately within subgroups or by targeting subgroup balance using exact matching on subgroup.

# Approach 1: match within subgroups
mA <- matchit(treat ~ x1 + x2 + x3, data = subset(data, mod == "A"))
mdA <- match.data(m1, data = subset(data, mod == "A"))

mB <- matchit(treat ~ x1 + x2 + x3, data = subset(data, mod == "B"))
mdB <- match.data(m1, data = subset(data, mod == "B"))

md <- rbind(mdA, mdB)

# Approach 2: exact match on subgroup
m <- matchit(treat ~ mod * (x1 + x2 + x3), data = data,
             exact = ~mod)
md <- match.data(m, data = data)

If you are using approach 1, you can assess balance using summary() or cobalt::bal.tab() for each matchit object. For approach 2, you should use cobalt::bal.tab(m, cluster = "mod"), which produces a balance table for each level of the moderator.

If balance is achieved, you can estimate the treatment effect using a Cox model:

fitA <- coxph(Surv(time, event) ~ treat * mod, data = md,
              weights = weights, cluster = subclass)
fitA

md$mod <- relevel(md$mod, "B")
fitB <- coxph(Surv(time, event) ~ treat * mod, data = md,
              weights = weights, cluster = subclass)
fitB

Again, the coefficient on treatment is the log hazard ratio for the reference group, which we change in the second model to the second group. The p-vale for the interaction coefficient should be the same. Remember only to interpret the tests on the coefficients.

Answer 2

If you do the analysis separately for men and women then you cannot get statistics about the interaction -- no p values or confidence intervals or standard errors. On the other hand, the output may be a bit easier to interpret (interpreting interactions is a frequent source of confusion for many and, even if it doesn't confuse the data analyst, it may confuse the audience).

But I don't see a problem with using a different matching method for men and women, matching within sex, and then analyzing the data as a whole. (I might be missing something, why do you think this might be a problem?)