Scalable proximal methods for cause-specific hazard modeling with time-varying coefficients

Abstract

Survival modeling with time-varying coefficients has proven useful in analyzing time-to-event data with one or more distinct failure types. When studying the cause-specific etiology of breast and prostate cancers using the large-scale data from the Surveillance, Epidemiology, and End Results (SEER) Program, we encountered two major challenges that existing methods for estimating time-varying coefficients cannot tackle. First, these methods, dependent on expanding the original data in a repeated measurement format, result in formidable time and memory consumption as the sample size escalates to over one million. In this case, even a well-configured workstation cannot accommodate their implementations. Second, when the large-scale data under analysis include binary predictors with near-zero variance (e.g., only 0.6% of patients in our SEER prostate cancer data had tumors regional to the lymph nodes), existing methods suffer from numerical instability due to ill-conditioned second-order information. The estimation accuracy deteriorates further with multiple competing risks. To address these issues, we propose a proximal Newton algorithm with a shared-memory parallelization scheme and tests of significance and nonproportionality for the time-varying effects. A simulation study shows that our scalable approach reduces the time and memory costs by orders of magnitude and enjoys improved estimation accuracy compared with alternative approaches. Applications to the SEER cancer data demonstrate the real-world performance of the proximal Newton algorithm.

Appendix

This appendix is devoted to the derivation of the gradient $\nabla \ell _j(\varvec{\gamma }_j)$ and Hessian matrix $\nabla ^2 \ell _j(\varvec{\gamma }_j)$ of $\ell _j(\varvec{\gamma }_j)$ as in (4). We define

$$\begin{aligned} S_{ij}^{(u)}(\varvec{\gamma }_j, X_i) :=\sum _{r \in R(X_i)} \exp \{[\mathbf {Z}_r \otimes \mathbf {B}(X_{i})]^\top \varvec{\gamma }_j\} \mathbf {Z}_r^{\odot u}, \quad u = 0, 1 ,2, \end{aligned}$$

where for a vector $\mathbf {v}\in {\mathbb {R}}^p$, $\mathbf {v}^{\odot 0} :=1$, $\mathbf {v}^{\odot 1} :=\mathbf {v}$, and $\mathbf {v}^{\odot 2} :=\mathbf {v}\mathbf {v}^\top$. The gradient $\nabla \ell _j(\varvec{\gamma }_j)$ and Hessian $\nabla ^2\ell _j(\varvec{\gamma }_j)$ of $\ell _j(\varvec{\gamma }_j)$ are hence given by

$$\begin{aligned} \nabla \ell _j(\varvec{\gamma }_j)&=\frac{1}{n}\sum _{i=1}^{n} \varDelta _{ij} \left\{ \mathbf {Z}_i - {\overline{\mathbf {Z}}}_{ij} (\varvec{\gamma }_j, X_i) \right\} \otimes \mathbf {B}(X_i), \end{aligned}$$

(10)

$$\begin{aligned} \nabla ^2\ell _j(\varvec{\gamma }_j)&=-\frac{1}{n}\sum _{i=1}^{n} \varDelta _{ij} \mathbf {V}_{ij}(\varvec{\gamma }_j, X_i) \otimes \left\{ \mathbf {B}(X_i) \mathbf {B}^\top (X_i) \right\} , \end{aligned}$$

(11)

in which

$$\begin{aligned} {\overline{\mathbf {Z}}}_{ij}(\varvec{\gamma }_j, X_i) :=\frac{S_{ij}^{(1)}(\varvec{\gamma }_j, X_i)}{S_{ij}^{(0)}(\varvec{\gamma }_j, X_i)}, \quad \mathbf {V}_{ij}(\varvec{\gamma }_j, X_i) :=\frac{S_{ij}^{(2)}(\varvec{\gamma }_j, X_i)}{ S_{ij}^{(0)}(\varvec{\gamma }_j, X_i)} - {\overline{\mathbf {Z}}}^{\odot 2}_{ij}(\varvec{\gamma }_j, X_i). \end{aligned}$$