3
$\begingroup$

Plenty of literature exists on how to interpret subgroup analyses in clinical trials. One example is in this thread. Unfortunately, I have not found any paper explaining how to perform a subgroup analysis on a low-technical level. I think that many people calculate the treatment effect and its confidence interval separately for each subgroup stratum. However, that isn’t correct, is it?

I wrote an example and hope that someone can complete it. I do not know how to calculate the confidence interval of the treatment effect for subgroup = 1. I guess that one have to account for the uncertainty introduced by the interaction term. This will probably not done as I did it in the following example. I am looking forward to your help. References to low-technical papers are also highly appreciated.

library(survival)
data(veteran)

veteran$trt <- veteran$trt-1
veteran$subgroup <- ifelse(veteran$age > 65, 1, 0)
table(treatment=veteran$trt, subgroup=veteran$subgroup)

fit <- coxph(Surv(time,status) ~ trt*subgroup, data=veteran)
s <- summary(fit)
s

# P value for interaction
round(s$coefficients['trt:subgroup','Pr(>|z|)'], 2)

### I think wrong...
fit2a <- coxph(Surv(time,status) ~ trt, data=veteran[which(veteran$subgroup==0),])
s2a <- summary(fit2a)
s2a
est <- s2a$coefficients['trt','coef']
se  <- s2a$coefficients['trt','se(coef)']
paste(round(exp(est),2), "(",
      round(exp(est - 1.96 * se), 2), "-",
      round(exp(est + 1.96 * se), 2), ")")
# age <= 65 years: HR 0.87, 95% CI 0.57 - 1.32

fit2b <- coxph(Surv(time,status) ~ trt, data=veteran[which(veteran$subgroup==1),])
s2b <- summary(fit2b)
est <- s2b$coefficients['trt','coef']
se  <- s2b$coefficients['trt','se(coef)']
paste(round(exp(est),2), "(",
      round(exp(est - 1.96 * se), 2), "-",
      round(exp(est + 1.96 * se), 2), ")")
# age > 65 years: HR 1.84, 95% CI 0.92 - 3.70


### I think correct...
# hazard ratio for treatment 0 versus 1 in subgroup = 0 (age <= 65 years)
est.0 <- s$coefficients['subgroup','coef']
se.0  <- s$coefficients['subgroup','se(coef)']
paste(round(exp(est.0),2), "(",
      round(exp(est.0 - 1.96 * se.0), 2), "-",
      round(exp(est.0 + 1.96 * se.0), 2), ")")
# HR 1.01, 95% CI 0.57 - 1.78

# hazard ratio for treatment 0 versus 1 in subgroup = 1 (age > 65 years)
s$coefficients
est.0 <- s$coefficients['subgroup','coef'] + s$coefficients['trt:subgroup','coef']
se.0 <- s$coefficients['subgroup','se(coef)'] + s$coefficients['trt:subgroup','se(coef)']
paste(round(exp(est.0),2), "(",
      round(exp(est.0 - 1.96 * se.0), 2), "-",
      round(exp(est.0 + 1.96 * se.0), 2), ")")
# HR 2.01, 95% CI 0.52 - 7.85

UPDATE #1 ... corrected after EdM's comment

@AdamO: Thanks for your reply!

a) You wrote: “[…] note that the stratum specific diabetic vs. non-diabetic pump-hazard ratios are non-overlapping (should be statistically significant), but the p-interaction > 0.05”. However, as far as I read the plot both confidence intervals overlap considerably.

Diabetes = yes : HR 0.80, 95% CI 0.63 to 1.01

Diabetes = no: HR 1.02, 95% CI 0.82 to 1.28

b) I came up with a different option (#3) to calculate the treatment effect and its confidence interval for subgroup = 1. it gives exact the same results as option #2. Statistics is fascinating, even when I do not understand these formulas ;-)

c) What you call "repurposed" is rather a standard procedure in randomized controlled trials. The essence of my questions is: How do I get these numbers presented in such plot. Personally, I prefer option 2/3 ... but this is only a gut feeling.

You wrote that you prefer option #1. Why should one calculate treatment effects and respective confidence interval from stratum-specific models? Shouldn’t we use the baseline hazard function of the total population? Does a consensus consists on this subject?

Option #1

Subgroup = 0: HR 0.87, 95% CI 0.57 - 1.32

Subgroup = 1: HR 1.84, 95% CI 0.92 - 3.70

Option #2

Subgroup = 0: HR 0.85, 95% CI 0.56 - 1.30

Subgroup = 1: HR 1.70, 95% CI 0.87 - 3.32

Option #3

Subgroup = 0: HR 0.85, 95% CI 0.56 - 1.30

Subgroup = 1: HR 1.70, 95% CI 0.87 - 3.32

The p value for interaction would be 0.09 in every case.

What is the correct option to calculate numbers for such a forest plot and why?

########################
# Loading data
########################
library(survival)
data(veteran)

########################
# Preparing data
########################
veteran$trt <- veteran$trt-1
veteran$subgroup <- ifelse(veteran$age > 65, 1, 0)
table(treatment=veteran$trt, subgroup=veteran$subgroup)

########################
# Fitting interaction model
########################
fit <- coxph(Surv(time,status) ~ trt*subgroup, data=veteran)
s <- summary(fit)
s

########################
# Extracting P value for interaction
########################
round(s$coefficients['trt:subgroup','Pr(>|z|)'], 2)

########################
# Option #1: Getting treatment estimates and respective confidence intervals for subgroup stata
# - the wrong way?
########################

fit2a <- coxph(Surv(time,status) ~ trt, data=veteran[which(veteran$subgroup==0),])
s2a <- summary(fit2a)
s2a
est <- s2a$coefficients['trt','coef']
se  <- s2a$coefficients['trt','se(coef)']
paste(round(exp(est),2), "(",
      round(exp(est - 1.96 * se), 2), "-",
      round(exp(est + 1.96 * se), 2), ")")
# age <= 65 years: HR 0.87, 95% CI 0.57 - 1.32

fit2b <- coxph(Surv(time,status) ~ trt, data=veteran[which(veteran$subgroup==1),])
s2b <- summary(fit2b)
est <- s2b$coefficients['trt','coef']
se  <- s2b$coefficients['trt','se(coef)']
paste(round(exp(est),2), "(",
      round(exp(est - 1.96 * se), 2), "-",
      round(exp(est + 1.96 * se), 2), ")")
# age > 65 years: HR 1.84, 95% CI 0.92 - 3.70


########################
# Option #2: Getting treatment estimates and respective confidence intervals for subgroup stata
# - correct or even worse?
########################

est.0 <- s$coefficients['trt','coef']
se.0  <- s$coefficients['trt','se(coef)']
paste(round(exp(est.0),2), "(",
      round(exp(est.0 - 1.96 * se.0), 2), "-",
      round(exp(est.0 + 1.96 * se.0), 2), ")")
# age <= 65 years: HR 0.85, 95% CI 0.56 - 1.30

# reversing the subgroup indicator
veteran$subgroup_reversed <- as.numeric(veteran$subgroup == 0)
table(subgroup=veteran$subgroup, subgroup_reversed=veteran$subgroup_reversed)
fit_reversed <- coxph(Surv(time,status) ~ trt*subgroup_reversed, data=veteran)
s_reversed <- summary(fit_reversed)
s_reversed
est.1 <- s_reversed$coefficients['trt','coef']
se.1  <- s_reversed$coefficients['trt','se(coef)']
paste(round(exp(est.1),2), "(",
      round(exp(est.1 - 1.96 * se.1), 2), "-",
      round(exp(est.1 + 1.96 * se.1), 2), ")")
# age > 65 years: HR 1.7 ( 0.87 - 3.32 )


########################
# Option #3: Calculating treatment effect and confidence interval for subgroup = 1 from the original model
# - correct or even worse?
########################
s$coefficients
vcov(fit)

est.r <- s$coefficients['trt', 'coef'] + s$coefficients['trt:subgroup', 'coef'] # trt = 1 and subgroup = 1
se.r <- sqrt(s$coefficients['trt', 'se(coef)']^2 + s$coefficients['trt:subgroup', 'se(coef)']^2 + 2 * vcov(fit)['trt','trt:subgroup'])

paste(round(exp(est.r),2), "(",
      round(exp(est.r - 1.96 * se.r), 2), "-",
      round(exp(est.r + 1.96 * se.r), 2), ")")

# age > 65 years: 1.70 ( 0.87 - 3.32 )
Improve this question
$\endgroup$
8
  • $\begingroup$ You may want have a look here: stats.stackexchange.com/questions/401341/… and here stats.stackexchange.com/questions/400739/… $\endgroup$
    – user213325
    Apr 6, 2019 at 16:13
  • $\begingroup$ As a note questions about programming are off-topic for this site. We can explain how the Cox model is specified, but if the specifications are unclear, looking at ?coxph or asking on stackoverflow.com is a good idea. $\endgroup$
    – AdamO
    Apr 9, 2019 at 14:20
  • 1
    $\begingroup$ With respect to overlapping confidence intervals (CI) in the figure to which you have linked, note that 95% CI can overlap even if the point estimates are statistically distinguishable at p < 0.05. Non-overlap of 95% confidence intervals is equivalent to requiring something like p < 0.005. See this page for details about CI versus statistical distinguishability. What's more appropriate to check is whether the point estimate for each of the classes lies outside the 95% CI for the other. In the case noted by @AdamO they are. $\endgroup$
    – EdM
    Apr 11, 2019 at 19:38
  • 1
    $\begingroup$ Once the code is correct, you will see that Options 2 and 3 give the same HR and CI, as expected. The remaining issue is the choice between stratum-specific versus interaction modeling, as the answer by @AdamO points out. I don't know whether there is a consensus on which is the better choice. There can be a tradeoff between robustness (stratified) and the power potentially provided by an interaction model (although then the model must be correctly specified and meet the necessary assumptions). What's most important is to be clear in your presentation about the choice that you made. $\endgroup$
    – EdM
    Apr 11, 2019 at 20:55
  • 1
    $\begingroup$ Such people probably used a standard interaction term as you did, although it is possible to evaluate an interaction of a treatment with a stratified covariate in a Cox model. See page 482 of Harrell's Regression Modeling Strategies, 2nd edition. I strongly recommend that book and the associated rms package in R if you are going to pursue regression modeling. There can be an initially steep learning curve with these, but the payoff is well worth it. $\endgroup$
    – EdM
    Apr 11, 2019 at 21:28

1 Answer 1

Reset to default
4
$\begingroup$

The plot shown in the stackoverflow post is called a forest plot, it is part of the meta package. A meta analysis is a way of comparing treatment effects estimated from multiple studies. But the authors of the plot in the original post have "repurposed" this forest plot function for another purpose.

On the X-axis is the hazard ratio for survival comparing on-pump CABG to off-pump CABG... whatever that means. For each item on the Y-axis, a subgroup specific hazard ratio is presented.

There are two ways to calculate subgroup specific hazard ratios. You've identified both of them. I prefer the first approach because it's more robust. This approach calculates stratum-specific hazard ratios. Each stratum has a unique baseline hazard function. However, there are no simple methods to obtain a p-value for interaction in this case (a bootstrap or $\delta$-method could be used).

It is expected that the p-value for interaction comes from a third interaction model, which you have fit with model statement Surv(time,status) ~ trt*subgroup, data=veteran. The interaction term here does not necessarily agree with the subgroup specific findings because the interaction model uses the whole sample to calculate the baseline hazard function whereas the stratum-specific analyses use a fraction of the sample. But in most cases (well powered analyses), it is sloppy but interpretable.

To see how they can be inconsistent, note that the stratum specific diabetic vs. non-diabetic pump-hazard ratios are non-overlapping (should be statistically significant), but the p-interaction > 0.05.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.