Calculation of calculated HR from coxph object with time-dependent coefficient and splines

Question

I want to calculate the estimated risk ratio as a function of time in the case of the coxph model with a time-dependent coefficient based on a spline term. I created a time-dependent coefficient using the tt function, similar to this example, which comes directly from ?coxph :

 # Fit a time transform model using current age cox = coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung, tt=function(x,t,...) pspline(x + t/365.25)) 

The call to survfit(cox) leads to an error that survfit does not understand the model with the term tt ( as described in 2011 by Terry Terno ).

You can extract the linear predictor using cox$linear.predictors , but I need to somehow extract the age and less trivially, once to go with each. Since tt splits the data set into the event time, I cannot just map the columns of the input frame to the output of coxph . In addition, I really would like to build an estimation function, not just forecasts for observed data points.

There is a related issue involving splines, but it does not include tt .

Edit (7/7)

I still stick to that. I studied this object in detail:

 spline.obj = pspline(lung$age) str(spline.obj) # something that looks very useful, but I am not sure what it is # cbase appears to be the cardinal knots attr(spline.obj, "printfun") function (coef, var, var2, df, history, cbase = c(43.3, 47.6, 51.9, 56.2, 60.5, 64.8, 69.1, 73.4, 77.7, 82, 86.3, 90.6)) { test1 <- coxph.wtest(var, coef)$test xmat <- cbind(1, cbase) xsig <- coxph.wtest(var, xmat)$solve cmat <- coxph.wtest(t(xmat) %*% xsig, t(xsig))$solve[2, ] linear <- sum(cmat * coef) lvar1 <- c(cmat %*% var %*% cmat) lvar2 <- c(cmat %*% var2 %*% cmat) test2 <- linear^2/lvar1 cmat <- rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2, 1, 1 - pchisq(test2, 1)), c(NA, NA, NA, test1 - test2, df - 1, 1 - pchisq(test1 - test2, max(0.5, df - 1)))) dimnames(cmat) <- list(c("linear", "nonlin"), NULL) nn <- nrow(history$thetas) if (length(nn)) theta <- history$thetas[nn, 1] else theta <- history$theta list(coef = cmat, history = paste("Theta=", format(theta))) } 

So, I have nodes, but I'm still not sure how to combine coxph coefficients with nodes to actually build the function. Any findings are greatly appreciated.

Answer 1

I think that what you need can be generated by creating an input matrix using pspline and multiplying the matrix by the corresponding coefficients from the output of coxph . To get HR, you need to take an exhibitor.

i.e.

 output <- data.frame(Age = seq(min(lung$age) + min(lung$time) / 365.25, max(lung$age + lung$time / 365.25), 0.01)) output$HR <- exp(pspline(output$Age) %*% cox$coefficients[-1] - sum(cox$means[-1] * cox$coefficients[-1])) library("ggplot2") ggplot(output, aes(x = Age, y = HR)) + geom_line() 

Please note that age is the age at the time of interest (i.e. the sum of the base age and elapsed time since entering the study). It should use the range specified to match the parameters of the original model. It can also be calculated using output x when using x = TRUE , as shown:

 cox <- coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung, tt=function(x,t,...) pspline(x + t/365.25), x = TRUE) index <- as.numeric(unlist(lapply(strsplit(rownames(cox$x), "\\."), "[", 1))) ages <- lung$age[index] output2 <- data.frame(Age = ages + cox$y[, 1] / 365.25, HR = exp(cox$x[, -1] %*% cox$coefficients[-1] - sum(cox$means[-1] * cox$coefficients[-1])))