Plotting smoothed hazard ratio intervals for interaction terms

Question

I don't know if it is possible.

I am following Terry Therneau's Spline terms in a Cox model vignette available for survival package.

In Section 3, Splines in an interaction, he shows how to visualise the interaction between sex and age where age is included in the model using splines.

The command he uses are

library(survival)
library(splines)
nfit3 <- coxph(Surv(futime, death) ~ sex * ns(age, df=3), flchain)
pdata <- expand.grid(age= 50:99, sex=c("F", "M"))
ypred <- predict(nfit3, newdata=pdata, se=TRUE)
yy <- ypred$fit + outer(ypred$se, c(0, -1.96, 1.96), '*')
matplot(50:99, exp(matrix(yy, ncol=6)), type='l', lty=c(1,1,2,2,2,2),
        lwd=2, col=1:2, log='y',
        xlab="Age", ylab="Relative risk")
legend(55, 20, c("Female", "Male"), lty=1, lwd=2, col=1:2, bty='n')
abline(h=1)

Which end up with the with plot with the relative risk with respect average population

I would like to have a plot where x axis shows age and y axis shows the effect of being male of certain age with respect being a female of the same age. Is that possible?

Answer 1

I've just seen this question made by myself one year ago that I managed to deal with recently, I hope correctly. I put the script in case that it can be of help to another person and to see if somebody detects something wrong.

Following the example given in the question we can calculate the hazard ratio and their 95% interval with

age.spline = ns(flchain$age, df=3)

X.p = cbind(pdata$sex == 'M',
            predict(age.spline, pdata$age),
            predict(age.spline, pdata$age) * (pdata$sex == 'M'))

lHR = apply(matrix(rowSums( t(t(X.p) * coef(nfit3)) ), ncol = 2), 1, diff)
AGE.pred = (X.p[X.p[,1]==1,] - X.p[X.p[,1]==0,])
s2.age = apply(AGE.pred, 1, function(x){
  t(x) %*% vcov(nfit3) %*% x
})

plot(50:99, exp(lHR), type = 'l', ylab = 'Hazard ratio', xlab = 'Age', lwd=2)
points(50:99, exp(lHR - 1.96 * sqrt(s2.age)), type= 'l', lty = 2, lwd=2)
points(50:99, exp(lHR + 1.96 * sqrt(s2.age)), type= 'l', lty = 2, lwd=2)
abline(h=1)

which ens up with the following plot

---

Extended after a question appeared in ResearchGate.

The approach can be further generalized to other packages. For example the package mgcv, which allows the user to fit an additive cox model with certain smooth terms (p.e. thin plate regression splines).

library(mgcv)
nfittp = gam(futime~sex+s(age,by=sex,bs='tp'), 
             data=flchain, weights = death, 
             family = 'cox.ph', method = 'REML')

df.male  = expand.grid(sex = 'M', age = 50:99)
df.female = expand.grid(sex = 'F', age = 50:99)

y0.male_ = predict(nfittp,newdata=df.male,type="lpmatrix")
y0.female_ = predict(nfittp,newdata=df.female,type="lpmatrix")

y0.male = predict(nfittp,newdata=df.male,type="link", se.fit=TRUE)
y0.female = predict(nfittp,newdata=df.female,type="link", se.fit=TRUE)

dplot = data.frame(
  age = 50:99,
  b.male = y0.male$fit - y0.female$fit,
  s.male = sapply(1:length(50:99), function(i){
    l = matrix(y0.male_[i,] - y0.female_[i,])
    t(l) %*% vcov(nfittp) %*% l
  })
)


par(mfrow=c(1,2))
plot(50:99, exp(y0.female$fit), type='l', log='y', lwd=2, xlab="Age", ylab="Relative risk", ylim=c(0.1,50))
points(50:99, exp(y0.female$fit + 1.96 * y0.female$se.fit), type='l', col=1, lwd=2, lty=2)
points(50:99, exp(y0.female$fit - 1.96 * y0.female$se.fit), type='l', col=1, lwd=2, lty=2)
points(50:99, exp(y0.male$fit), type='l', col='red', lwd=2)
points(50:99, exp(y0.male$fit + 1.96 * y0.male$se.fit), type='l', col=2, lwd=2, lty=2)
points(50:99, exp(y0.male$fit - 1.96 * y0.male$se.fit), type='l', col=2, lwd=2, lty=2)
legend(50, 40, c("Female", "Male"), lty=1, lwd=2, col=1:2, bty='n')
abline(h=1)

plot(50:99, exp(dplot$b.male), type = 'l', ylab = 'Hazard ratio', xlab = 'Age', lwd=2, ylim = c(0.8,2))
points(50:99, exp(dplot$b.male - 1.96 * sqrt(dplot$s.male)), type= 'l', lty = 2, lwd=2)
points(50:99, exp(dplot$b.male + 1.96 * sqrt(dplot$s.male)), type= 'l', lty = 2, lwd=2)
abline(h=1)

Answer 2

Wow it is really a surprising plot. But i am interested that the confidence interval in the second plot is equal to the real CI of hazard ratio in univariate K-M analysis for age. Besides, is there a method could i include several other variates in the model to get a multi variate COX analysis and the corresponding HR and CI. I am looking forward to the answer,Thanks very much.