3
$\begingroup$

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 enter image description here

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?

$\endgroup$
1
$\begingroup$

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 enter image description here


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)

enter image description here

$\endgroup$
  • 2
    $\begingroup$ This is done easily with the R rms package cph and contrast.rms functions. $\endgroup$ – Frank Harrell May 19 at 3:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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