What is the reference point for a Cox proprtional hazard spline model?

Question

When trying to fit a Cox model to a continuous term we have a problem interpreting the result. When fitting a linear model the hazard ratio is referenced to the mean of the predictor values, and at that point the SE is zero. In comparison, when fitting a spline model this is not the case, and the SE is non-zero at all points and it is therefore not clear what's the reference point for the hazard ratio.

The code below shows the problem -- on the right the linear model shows the SE becoming 0 at the mean (x=0 in this case). However on the left, the spline model has no point where the SE is zero

We think it may be due to the fact that the location of the means of the spline matrix do not coincide, however we couldn't find how to correct the spline matrix to have the means coincide and identify the reference point.

require("Hmisc")
require("survival")

N <- 50
x <- seq(from = -1, to = 1, length.out = N)
par(mfrow = c(1, 2))

# fitting a spline model
xx <- rcspline.eval(x, nk = 4, inclx = TRUE)
y <- rexp(N, exp(2 + 2 * xx[, 1] + 3 * xx[, 2] + 4 * xx[, 3]))
l <- (coxph(Surv(y) ~ xx))
y2 <- predict(l, se.fit = TRUE)


plot(y2$fit ~ xx[, 1], type = "l")
lines(y2$fit + y2$se ~ xx[, 1], type = "l")
lines(y2$fit - y2$se ~ xx[, 1], type = "l")
grid()

# fitting a linear model
y <- rexp(N, exp(2 + 2 * x))
l <- (coxph(Surv(y) ~ x))
y2 <- predict(l, se.fit = TRUE)

plot(y2$fit ~ x, type = "l")
lines(y2$fit + y2$se ~ x, type = "l")
lines(y2$fit - y2$se ~ x, type = "l")
grid()

Answer 1

I'm not sure if this will give you the answer, but it's just a guess.

As reading from Therneau's "A package for survival analysis in S", what you actually get is the $\eta_i = \beta' x_i$ where $\beta'$ is a vector of estimated log-hazard ratios. At $x_i = 0$ this will always be 0 for the linear model. The variance of this is $x_i' V x_i$ where $V$ is the covariance matrix of the $\beta$. If $x_i$ is 0, then also the variance of the linear predictor will be 0 by this.

Replicating your code for the splines part, even when the first column of xx is 0, the other two are not, i.e. there is no 0 0 0 row in the xx matrix, so the calculated variance matrix does not have a 0 element.

The built-in pspline function plays well with coxph for time-dependent effects and is quite well documented.

Answer 2

Following the answer by @Theodor we wrote a simplified predict function for this case. Using it will predict the hazard ratio relative to the ref.point

predict.zero <- function(model, x, ref.point = mean(x, na.rm = TRUE)) {

  ## offset so that each column has a zero at the location of the reference point
  idx.mean <- which.min(abs(x[, 1] - ref.point))
  for (column in 1:ncol(x)) {
    x[, column] <- x[, column] - x[idx.mean, column]
  }

  fit <- data.frame(
    se = apply(x, 1, function(v) {
      return(sqrt(t(v) %*% model$var %*% v))
    })
  )

  if (length(model$coefficients) == 1 ) {
      fit$fit <- x * model$coefficients
  }
   else {
      fit$fit <- x %*% model$coefficients
   }
   return(fit)
}