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 )
?coxph
or asking on stackoverflow.com is a good idea. $\endgroup$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$