Example

How is higher than normal breast density related to breast cancer?

Is the risk related to

• long periods of higher than normal density
• the current density only, or
• large changes over a short period of time?

We need a model for the density, the marker, and the risk of breast cancer, the survival time.

Markers

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

The deviation from the markers from the population may be linked to higher or lower risk of the event for the individuals than the population.

Different ways the two can be linked.

Hazards

\(T_i\) is the event time of individual \(i\). The hazard is defined as

\[ h_i(t) = \lim_{dt\rightarrow0}\frac{P(t\leq T_i + dt\mid T_i \geq t)}{dt} \]

Given the hazard, \(h_i\), we can compute the density of \(T_i = t\) and probabilities for events like \(T_i > t\) and \(T_i\in (t, v)\).

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*} \]

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

Survival Model

Can replace \(\alpha\vec m(t)^\top\vec U_i\) with

\[\alpha\left(\frac{\partial\vec m(t)}{\partial t}\right)^\top\vec U_i\]

or

\[\alpha\left(\int_0^t \vec m(s) ds\right)^\top\vec U_i.\]

… we can also add them.

Informative Visiting 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

\[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}\]

and \(M(s)\) and \(X_i\) 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*} \]

The Likelihood

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

\(l_i(\vec\zeta)\) is of the form

\[ \exp(l_i(\vec\zeta)) = E\left( g_i(\vec U_i)h_{ij}(t\mid \vec U_i, \vec\xi_i) \exp\left(-\int_0^t h_{ij}(s\mid \vec U_i, \vec\xi_i)ds\right)\right) \]

Variational Approximations

Given distribution \(q^{(i)}_{\vec\theta}\) for \(\theta\in \Theta_i\) and outcomes \(\vec y_i\)

\[ \begin{align*} l_i(\vec\zeta) &= \log p_i(\vec y_i) = \log \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\\ &= \log \int q^{(i)}_{\vec\theta}(\vec U) \left( \log \left(\frac{p_i(\vec y_i, \vec U)} {q^{(i)}_{\vec\theta}(\vec U)}\right) d \vec U + \log \left(\frac{q^{(i)}_{\vec\theta}(\vec U)} {p_i(\vec U\mid\vec y_i)}\right) \right)d \vec U \\ &\geq \log \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

We assume that

\[ \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.

The Variational Distribution

For a given variational distribution \(q^{(i)}_{\vec\theta}\), we need to be able to easily evaluate

• the entropy: \(-E_{q^{(i)}_{\vec\theta}}(\log q^{(i)}_{\vec\theta}(\vec U))\).
• the first and second moment: \(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))\).

Not strong requirements!

Implementation

Implemented with a multivariate normal distribution.

Implementation

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

Supports left-truncation, partly observed markers, and much more.

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

There is one (terminal) time to event outcome for each individual.

We sample 300 individuals.

Results

The estimated associations parameters, \(\vec\alpha_1\), are shown below.

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

Two different types of markers observed possible on the same points in time.

One (terminal) time to event outcome for each individual and a model for the observation times.

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 \(111.299 \pm 2.943\) seconds. 100 data sets are sampled.

Results

Bias estimates for the fixed effects of the survival outcomes, \(\vec\delta_1\) and \(\vec\delta_2\), are shown above multiplied by 100. 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 multiplied by 100. 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 multiplied by 100. 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 multiplied by 100. 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 multiplied by 100. 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 will be interesting.

Implementation

The software will soon be public on Github.

Variational Approximations