Rizopoulos 2016 - Joint Model Specification and Posterior

Summary

Observed data are $(T_i,{\delta }_i,{\mathbf{y}}_i)$. The longitudinal outcome follows a GLMM with linear predictor ${\eta }_i(t)={\mathbf{x}}_i^{\top }\mathbf{\beta }+{\mathbf{z}}_i^{\top }{\mathbf{b}}_i$; the hazard is $h_i(t)=h_0(t)\exp \{{\mathbf{\gamma }}^{\top }{\mathbf{w}}_i+f[H_i(t),{\mathbf{b}}_i,\mathbf{\alpha }]\}$. Conditional on ${\mathbf{b}}_i$ the processes are independent, which gives a posterior factorisation. The survival integral has no closed form and is computed by Gauss-Kronrod or Gauss-Legendre quadrature.

Main Content

Joint model (Eqs. 1-3)

Notation: $T_i^{\ast }$ true event time, $C_i$ censoring, $T_i=\min (T_i^{\ast },C_i)$, ${\delta }_i=I(T_i^{\ast }\le C_i)$, ${\mathbf{y}}_i$ the $n_i$ longitudinal responses at times $t_{il}$. Longitudinal (GLMM, Eq. 1): $g(\mathbb{E}\{y_i(t)\mid {\mathbf{b}}_i\})={\eta }_i(t)={\mathbf{x}}_i(t{)}^{\top }\mathbf{\beta }+{\mathbf{z}}_i(t{)}^{\top }{\mathbf{b}}_i$, $g$ monotone link, $y$ from the exponential family, ${\mathbf{b}}_i\sim N(0,D)$ or multivariate $t_{df}(0,D)$. Survival (Eq. 2): $h_i(t\mid H_i(t),{\mathbf{w}}_i(t))=h_0(t)\exp \{{\mathbf{\gamma }}^{\top }{\mathbf{w}}_i(t)+f[H_i(t),{\mathbf{b}}_i,\mathbf{\alpha }]\}$, $H_i(t)=\{{\eta }_i(s),0\le s<t\}$ the marker history, ${\mathbf{w}}_i$ exogenous covariates, $\mathbf{\alpha }$ the association parameters. Baseline hazard (Eq. 3): $\log h_0(t)={\gamma }_{h_0,0}+{\sum }_{q=1}^Q{\gamma }_{h_0,q}{B}_q(t,\mathbf{v})$, B-splines with knots $\mathbf{v}$. Rule of thumb: total parameters between $1/10$ and $1/20$ of the number of events (Harrell 2001); or use 15-20 knots with a penalty (Eilers and Marx 1996).

Conditional independence and posterior (Eqs. 4-7)

Given ${\mathbf{b}}_i$ and $\theta$, longitudinal and event processes are independent and longitudinal measurements are mutually independent: $p({\mathbf{y}}_i,T_i,{\delta }_i\mid {\mathbf{b}}_i,\theta )=p({\mathbf{y}}_i\mid {\mathbf{b}}_i,\theta )p(T_i,{\delta }_i\mid {\mathbf{b}}_i,\theta )$ and $p({\mathbf{y}}_i\mid {\mathbf{b}}_i,\theta )={\prod }_{l}p(y_{il}\mid {\mathbf{b}}_i,\theta )$. Hence

$$p(\theta ,\mathbf{b})\propto \prod _{i=1}^n\{\prod _{l=1}^{n_i}p(y_{il}\mid {\mathbf{b}}_i,\theta )\}p(T_i,{\delta }_i\mid {\mathbf{b}}_i,\theta )p({\mathbf{b}}_i\mid \theta )p(\theta ),$$

with the exponential-family density $p(y_{il}\mid {\mathbf{b}}_i,\theta )=\exp \{[y_{il}{\psi }_{il}({\mathbf{b}}_i)-c({\psi }_{il}({\mathbf{b}}_i))]/a(\phi )-d(y_{il},\phi )\}$ and the survival part

$$p(T_i,{\delta }_i\mid {\mathbf{b}}_i,\theta )=h_i(T_i\mid H_i(T_i){)}^{{\delta }_i}\exp \{-{\int }_0^{T_i}{h}_i(s\mid H_i(s))ds\},S_i(t)=\exp \{-{\int }_0^t{h}_0(s)\exp [{\mathbf{\gamma }}^{\top }{\mathbf{w}}_i(t)+f\{H_i(s),\mathbf{\alpha }\}]ds\}\text{(Eq. 7)}.$$

The integral requires numerical quadrature (Gauss-Kronrod, Gauss-Legendre).

Priors. Independent diffuse normals for $\mathbf{\beta },\mathbf{\gamma },{\mathbf{\gamma }}_{h_0},\mathbf{\alpha }$. Penalised spline: $p({\gamma }_{h_0}\mid {\tau }_h)\propto {\tau }_h^{\rho (K)/2}\exp \{-\frac{{\tau }_h}{2}{\gamma }_{h_0}^{\top }K{\gamma }_{h_0}\}$, $K={\Delta }_r^{\top }{\Delta }_r+10^{-6}I$ (${\Delta }_r$ the $r$-th difference penalty, $\rho (K)$ its rank), ${\tau }_h\sim \text{Gamma}(1,0.005)$ (Lang and Brezger 2004). Inverse-Wishart for $D$, inverse-Gamma for the error variance ${\sigma }^2$.

Examples

Cumulative hazard by Gauss-Legendre

Setup: $\eta (t)=b_0+b_{1}t$, baseline $h_0=\lambda$, $f=\alpha \eta$: $h(t)=\lambda e^{\alpha (b_0+b_{1}t)}$. With $\lambda =0.05,b_0=0.5,b_1=0.3,\alpha =1.2,T=6$ the closed-form cumulative hazard is $\lambda e^{\alpha b_0}(e^{\alpha b_{1}T}-1)/(\alpha b_1)=1.941351$; a 15-point Gauss-Legendre rule reproduces it to $10^{-15}$. Interpretation: for smooth hazards a modest number of quadrature nodes per subject suffices, which is how E25 evaluates $H_i(T_i)$ inside PyMC.

import numpy as np
def cum_hazard(T, b0, b1, alpha, lam=0.05, k=15):
    x, w = np.polynomial.legendre.leggauss(k)
    s = 0.5*T*(x+1); ww = 0.5*T*w
    return (ww * lam * np.exp(alpha*(b0 + b1*s))).sum()

PyMC-style sketch (untested).

with pm.Model():
    chol, corr, sd = pm.LKJCholeskyCov("chol", n=2, eta=2, sd_dist=pm.HalfNormal.dist(1))
    b = pm.MvNormal("b", 0, chol=chol, shape=(N, 2))             # random intercept, slope
    eta = beta0 + beta1*t_obs + b[sub,0] + b[sub,1]*t_obs
    pm.Normal("y", eta, sigma, observed=y)
    alpha = pm.Normal("alpha", 0, 1)
    H = cum_hazard_gl(T, beta0 + b[:,0], beta1 + b[:,1], alpha, shape, scale)  # nodes via leggauss
    pm.Potential("surv", (delta*(log_h(T)) - H).sum())

Connections

Survival Analysis, Hierarchical Models, Rizopoulos 2016 - Implementation MCMC and Baseline Hazard, Rizopoulos 2016 - Association Structures and Transformation Functions.

See Also

Rizopoulos 2016 - JMbayes Joint Models Overview