# Plotting smoothed hazard ratio intervals for interaction terms

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?

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) • This is done easily with the R rms package cph and contrast.rms functions. – Frank Harrell May 19 at 3:25