I frequently do survival analysis on large data sets. One million samples or more is typical, and this seems to be much more than typical research usage. Many algorithms I've used are prohibitively slow at this scale.

Fortunately, the coxph routine in R's survival package is remarkably fast and scales much better with sample size than other packages. Why is this? Is there a reference for the algorithm implemented in coxph?

I'd like to understand so that I can investigate numerically efficient methods for time-dependent covariates on large data. Existing methods seem to be very slow (including the time transform functionality in survival) and I suspect that they may consume much more time and/or memory than necessary.

At the very least, it seems that some delicate algorithmic techniques have been used. Here's a relevant note from the time-dependent covariates vignette:

Although handy, the computational impact of the tt argument should be considered before using it. The Cox model requires computation of a weighted mean and variance of the covariates at each event time, a process that is inherently $O(ndp^2)$ where $n =$ the sample size, $d =$ the number of events and $p =$ the number of covariates. Much of the algorithmic effort in coxph() is to use updating methods for the mean and variance matrices, reducing the compute time to $O((n+d)p^2)$. When a tt term appears updating is not possible; for even moderate size data sets the impact of $nd$ versus $n + d$ can be surprising.

This is a tantalizing hint but I don't know how the overall algorithm works and thus cannot really understand it. I'd prefer a mathematical exposition over wading through the source, and I wasn't able to find one in the references for coxph in the survival package manual on CRAN.


1 Answer 1


Believe it or not, it's just Newton-Raphson. It's right here. The weighted mean and covariance matrices mentioned in the vignette passage are Equations (3.4) through (3.6).

  • $\begingroup$ Why are you surprised? Newton-Raphson is extremely fast; it has quadratic convergence. Given one can get a reasonable and not too expensive approximation to the Hessian N-R is probably his best bet to optimise a function. The whole battle is when one cannot get a Hessian and has to use quasi-Newton stuff (eg. BFGS) or derivative-free deep magic (aka. linear/quadric approximations - eg. NEWUOA). $\endgroup$
    – usεr11852
    Aug 7, 2015 at 2:54
  • $\begingroup$ I like Newton-Raphson, but it is not always fast after accounting for the cost of computing the Hessian, as in this case. Sometimes specialty approaches (or special variants of N-R) can be faster. In this case it's N-R with tricks to speed up the gradient/Hessian computation. $\endgroup$
    – Paul
    Aug 7, 2015 at 3:11
  • $\begingroup$ Agreed. I specifically said "not too expensive approximation" and then pointed to quasi-Newton methods. $\endgroup$
    – usεr11852
    Aug 7, 2015 at 3:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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