How to construct a design matrix for coxph with pspline term?

Question

I am wondering how to reconstruct the design matrix for a coxph() model with a pspline() term. For example, if I fit the following model,

fit <- coxph(Surv(t,delta) ~ pspline(x,df=0))

how can I calculate the linear predictor by hand for a given value of $x$?

Formally, I am modeling $h(t|x) = h_0(t)\cdot e^{(X^\intercal\beta)}$, and I want to know how to calculate $X$ or $X^\intercal\beta$ by hand. I can't find detailed documentation of how the design matrix is constructed when the pspline() term is used.

As for why I care, I am running a large set of new predictions from a coxph() object with a pspline() term. I have stored the cumulative hazard $H_0(t)$ and want to be able to calculate $H_0(t)\cdot e^{(X^\intercal\beta)}$ efficiently. I'm finding that using

predict(fit, newdata=data.frame(x=x.new), type='lp')

is very inefficient! Any tips would be greatly appreciated.

Answer 1

If you submit the following in R

> pspline

you will see the details of this function and how the design matrix is constructed via the spline.des() function from the 'splines' package.

From there, you can create new.pspline(newx, x) function that uses the actual x variable values to determine the arguments for the underlying spline.des() function, but then output the design matrix corresponding to newx values.

For example,

new.pspline <- function(newx,x){ 
    df = 4; nterm = 2.5 * df; degree = 3; eps = 0.1; Boundary.knots = range(x); nterm = round(nterm)
    dx <- (Boundary.knots[2] - Boundary.knots[1])/nterm
    knots <- c(Boundary.knots[1] + dx * ((-degree):(nterm - 1)), Boundary.knots[2] + dx * (0:degree)) 
    xx <- spline.des(knots, newx, degree + 1, outer.ok = TRUE)$design
    xx <- xx[, -1, drop = FALSE]
    xx
}