I am interested in plotting the HR from a coxph() model for a categorical-continuous interaction. Basically estimating the binary HR at multiple values of the continuous predictor and then plotting them. I can do the estimation bit, but I am interested to see if there's an easier way to identify and filter the relevant contrasts from the generation of all pairwise contrasts. The code below produces:
library(survival)
library(emmeans)
# Fit model with covariate-covariate interaction
fit <- coxph(Surv(time, status) ~ trt * karno, data = veteran)
summary(fit)
# Estimate trt HR at different values of karno
emm_df <- emmeans(fit, ~ trt + karno, at = list(karno = seq(0, 100, by = 5)), type = "unlink") |> pairs(rev = T) |> summary(infer = T)
head(emm_df)
contrast ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
trt2 karno0 / trt1 karno0 2.984 1.8095 Inf 0.2785 31.970 1 1.803 0.9996
trt1 karno5 / trt1 karno0 0.885 0.0346 Inf 0.7592 1.031 1 -3.139 0.4244
trt1 karno5 / trt2 karno0 0.296 0.1715 Inf 0.0309 2.848 1 -2.102 0.9922
trt2 karno5 / trt1 karno0 2.438 1.4338 Inf 0.2445 24.314 1 1.516 1.0000
trt2 karno5 / trt2 karno0 0.817 0.0263 Inf 0.7206 0.927 1 -6.285 <.0001
trt2 karno5 / trt1 karno5 2.756 1.5413 Inf 0.3095 24.550 1 1.813 0.9996
How would I filter only those contrasts of trt2/trt1 at the same values of karno? e.g. I only want to keep:
trt2 karno0 / trt1 karno0
trt2 karno5 / trt1 karno5
trt2 karno10 / trt1 karno10
etc
It would be even handier if I could convert this variable to numeric on the fly (ready for plotting) e.g. converting those categories above to 0,5,10,etc
Thanks
Use the | notation: emmeans(fit, ~ trt | karno, ...)