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

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.

$\endgroup$
1
  • $\begingroup$ Welcome to the site! Please do not use the "Your Answer" field to ask new or follow-up questions. Instead, click the blue ASK QUESTION link at the top of the page & ask there. If you want to address the author of an answer, you will need 50 reputation first to leave comments. Since you're new here, you may want to take our tour, which has information for new users. $\endgroup$ – vkehayas Sep 21 '20 at 7:39

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.