Example

How is higher than normal breast density and other variables related to breast cancer?

Can be estimated with a joint model for the density, the marker, and the risk of breast cancer, the survival time.

Goals

The goals is a fast estimation method for general joint marker and survival models.

Start by defining the model.

Markers

Individuals’ markers are observed at different points in time.

Formally

\[ \begin{align*} Y_{ij} &= \mu_i(s_{ij}, \vec U_i) + \epsilon_{ij} \\ \epsilon_{ij} &\sim N(0, \sigma^2) \\ \mu_{i}(s, \vec U_i) &= \vec x_i^\top\vec\gamma + \vec g(s)^\top\vec\beta + \vec m(s)^\top\vec U_i \\ \vec U_i &\sim N^{(R)}(\vec0, \Psi) \end{align*} \]

\(\vec g(s)^\top\vec\beta\): the population mean curve.

\(\vec m(s)^\top\vec U_i\): difference for individual \(i\) to the population curve.

Time-to-event Outcomes

We have time-to-event outcomes.

Markers’ deviation from the population may be linked to higher or lower risk of the event for each individual.

Different ways the two can be linked.

Survival Model

\[ \begin{align*} h_i(t\mid \vec U_i, \xi_i) &= \exp\left( \vec z_i^\top\vec \delta + \omega^\top \vec b(t) + \alpha \vec m(t)^\top\vec U_i + \xi_i \right) \\ \xi_i &\sim N(0, \Xi) \end{align*} \]

\(\vec\omega^\top \vec b(t)\): the log baseline hazard function.

Informative Observation Process

Such dependencies can be handled with joint models (Gasparini et al. 2020).

More Markers

We may \(L\) different markers rather than one

\[ \begin{align*} \vec Y_{ij}&= \vec\mu_i(s_{ij}, \vec U_i) + \vec\epsilon_{ij} \\ \epsilon_{ij} &\sim N^{(L)}(\vec 0, \Sigma) \\ \mu_{i1}(s, \vec U_i) &= \vec x_{i1}^\top\vec\gamma_1 + \vec g_1(s)^\top\vec\beta_1 + \vec m_1(s)^\top\vec U_{i1} \\ \vdots &\hphantom{=}\vdots\\\ \mu_{iL}(s, \vec U_i) &= \vec x_{iL}^\top\vec\gamma_L + \vec g_L(s)^\top\vec\beta_L + \vec m_L(s)^\top\vec U_{iL} \\ \vec U_i &= \begin{pmatrix} \vec U_{i1} \\ \vdots \\ \vec U_{iL} \end{pmatrix}\sim N^{(R)}(\vec0, \Psi). \end{align*} \]

More Survival Types

\[\begin{align*} G(s) &= \begin{pmatrix} \vec g_1(s)^\top & 0^\top & \cdots & \vec 0^\top \\ \vec 0^\top & \vec g_2(s)^\top & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vec 0^\top \\ \vec 0^\top & \cdots & \vec 0^\top & \vec g_L(s)^\top \end{pmatrix} \\ \vec\gamma &=(\vec \gamma_1^\top,\dots,\vec\gamma_L^\top)^\top, \end{align*}\]

and define \(M(s)\), \(X_i\), and \(\vec\beta\) similarly. Then

\[\vec\mu_i(s_{ij}, \vec U_i) = X_i\vec\gamma + G(s_{ij})\vec\beta + M(s_{ij})\vec U_i.\]

More Survival Types

We may have \(H\) different survival types rather than one

\[ \begin{align*} h_{i1}(t\mid \vec U_i, \vec\xi_i) &= \exp\left( \vec z_{i1}^\top\vec \delta_1 + \omega_1^\top \vec b_1(t) + \vec\alpha_1^\top M(t)\vec U_i + \xi_{i1} \right) \\ \vdots &\hphantom{=}\vdots \\ h_{iH}(t\mid \vec U_i,\vec \xi_i) &= \exp\left( \vec z_{iL}^\top\vec \delta_H + \omega_H^\top\vec b_H(t) + \vec\alpha_H^\top M(t)\vec U_i + \xi_{iH} \right) \\ \vec\xi_i = \begin{pmatrix}\xi_{i1} \\ \vdots \\ \xi_{iH}\end{pmatrix} &\sim N^{(H)}(\vec 0, \Xi) \end{align*} \]

Alternative Parameterization

Simplify the notation to

\[ \mu_{i}(s, \vec U_i) = \vec x_i(s)^\top\vec\gamma + \vec m(s)^\top\vec U_i \]

Could let the hazard depend on \(\mu_{i}(s, \vec U_i)\) rather than \(\vec m(s)^\top\vec U_i\).

Alternative Parameterization (Cont.)

Suppose \(\vec z_i(t) = (\widetilde{\vec z}_i(t)^\top, \vec x_i(t)^\top)^\top\). Then

\[ \begin{align*} h_i(t\mid \vec U_i, \xi_i) &= \exp\big( \vec z_i(t)^\top\vec \delta + \alpha \vec m(t)^\top\vec U_i + \xi_i \big) \\ &= \exp\big( \widetilde{\vec z}_i(t)^\top\vec \delta_1 + \vec x_i(t)^\top\vec (\vec\delta_2 - \alpha\vec\gamma) \\ &\hspace{40pt} + \alpha \mu_{i}(t, \vec U_i) + \xi_i \big) \end{align*} \]

\(\alpha\) has the same interpretation unless \(\vec x_i(t)\) is not a linear combination of \(\vec z_i(t)\).

Alternative Parameterization (Cont.)

The Likelihood

Let \(\vec\zeta\) be the model parameters and \(l_i(\vec\zeta)\) be the log likelihood term of individual \(i\).

Variational Approximations

Given a distribution with density \(q^{(i)}_{\vec\theta}\) for \(\theta\in \Theta_i\) and outcomes \(\vec y_i\), the particular variational approximation we use is the lower bound

\[ l_i(\vec\zeta) \geq \int q^{(i)}_{\vec\theta}(\vec U) \log \left(\frac{p_i(\vec y_i, \vec U)} {q^{(i)}_{\vec\theta}(\vec U)}\right) d \vec U = \tilde l_i(\vec\zeta, \vec\theta) \]

Typically, much faster to evaluate.

Variational Approximations (Cont.)

\[ \begin{align*} l_i(\vec\zeta) &= \log p_i(\vec y_i) = \int q^{(i)}_{\vec\theta}(\vec U) \log \left(\frac{p_i(\vec y_i, \vec U)\big/ q^{(i)}_{\vec\theta}(\vec U)} {p_i(\vec U\mid\vec y_i)\big/ q^{(i)}_{\vec\theta}(\vec U)}\right) d \vec U\\ &= \int q^{(i)}_{\vec\theta}(\vec U) \Bigg( \log \left(\frac{p_i(\vec y_i, \vec U)} {q^{(i)}_{\vec\theta}(\vec U)}\right) d \vec U + \underbrace{\log \left(\frac{q^{(i)}_{\vec\theta}(\vec U)} {p_i(\vec U\mid\vec y_i)}\right)}_{\text{KL divergence}} \Bigg)d \vec U \\ &\geq \int q^{(i)}_{\vec\theta}(\vec U) \log \left(\frac{p_i(\vec y_i, \vec U)} {q^{(i)}_{\vec\theta}(\vec U)}\right) d \vec U = \tilde l_i(\vec\zeta, \vec\theta) \end{align*} \]

Equality if \(p_i(\vec U\mid\vec y_i) = q^{(i)}_{\vec\theta}(\vec U)\).

Approximate Maximum Likelihood

\[ \text{argmax}_{\vec\zeta} \sum_{i = 1}^n l_i(\vec\zeta) \approx \text{argmax}_{\vec\zeta} \max_{\vec\theta_1,\dots,\vec\theta_n} \sum_{i = 1}^n\tilde l_i(\vec\zeta, \vec\theta_i) \]

The latter problem has many parameters but the problem is partially separable.

Implementation

Implemented with a multivariate normal distribution.

Written in C++ with support for computation in parallel.

Overview

Make a comparison with the Monte Carlo based joineRML package (Hickey, Philipson, and Jorgensen 2018).

Simulation study to check the bias of the VA.

Comparison with joineRML

One terminal time-to-event outcome for each individual.

We sample 300 individuals.

Results

Results

The estimated covariance matrix of the random effects, \(\Psi\), is shown above. The true values are shown below.

Details

Sobol sequences are used with joineRML.

Both markers have a standard normal distributed covariate and the time-to-event outcome has \(\text{Unif}(-1, 1)\) covariate.

Observation times are drawn with increments that are exponentially distributed with rate 0.5.

Second Simulation Study

A terminal time-to-event outcome.

An informative observation process.

We sample 1000 individuals.

Four random effects for the markers and two for the frailties, \(\vec\xi_i\), per individual.

Results

Average computation time was \(48.619 \pm 1.866\) seconds. 100 data sets are sampled.

Approximate Likelihood Ratio (cont.)

The average computation time was 127.147 seconds.

Computationally efficient approximation Wald intervals can be implemented.

Results

Bias estimates for the fixed effects of the survival outcomes, \(\vec\delta_1\) and \(\vec\delta_2\), are shown above. The SE columns are the standard errors. The true values are shown below.

Results

Bias estimates for the log baseline hazard of the survival outcomes, \(\vec\omega_1\) and \(\vec\omega_2\), are shown above. The SE columns are the standard errors. The true values are shown below.

Results

Bias estimates for the covariance matrix of the random effects, \(\Psi\), are shown above. The last columns show the standard errors. The true values are shown below.

Results

Bias estimates for the covariance matrix of the markers’ error term, \(\Sigma\), are shown above. The last columns show the standard errors. The true values are shown below.

Results

Bias estimates for the covariance matrix of the frailties, \(\Xi\), are shown above. The last columns show the standard errors. The rue values are shown below.

Details

\(\vec g_1\) and \(\vec m_1\) are natural cubic splines with boundary knots at \((0, 10)\). \(\vec g_1\) has interior knots at \(\{3.33,6.67\}\). The first marker has a standard normally distrusted covariate.

\(\vec g_2\) and \(\vec m_2\) are polynomials without and with an intercept and degree 2 and 1. The second marker has no covariates.

The terminal event has a B-spline basis with boundary knots at \((0,10)\) and an interior knot at \(5\). There is a \(\text{Unif}(-1, 1)\) covariate.

The observation process has a natural cubic spline with boundary knots at \((0,10)\) and an interior knot at \(5\). There is no covariates.

Conclusions

Introduced variational approximations.

Showed that they can yield fast estimates of joint survival and marker models.

Simulation studies suggest that the bias is low.

Models

Extension to competing events and multi-state models.

Time-varying covariate effects and interval censoring are simple extensions mathematically but not implemented.

Variational Approximations

The lower bound requires

• the entropy: \(-E_{q^{(i)}_{\vec\theta}}(\log q^{(i)}_{\vec\theta}(\vec U))\).
• the first two moments: \(E_{q^{(i)}_{\vec\theta}}(\vec U)\) and \(E_{q^{(i)}_{\vec\theta}}(\vec U\vec U^\top)\).
• the moment generating function: \(E_{q^{(i)}_{\vec\theta}}(\exp(\vec t^\top\vec U))\).