This vignette describes the mathematical framework and computational implementation of the JointODE model, which jointly models longitudinal biomarker trajectories and survival outcomes through a coupled ordinary differential equation (ODE) system.
Model Framework
Longitudinal Sub-Model
The observed biomarker measurements are modeled as: V_{ij}=m_i(T_{ij})+b_i+\varepsilon_{ij},\quad i=1,\ldots,n,\quad j=1,\ldots,n_i
where:
- V_{ij}: Observed biomarker value for subject i at time T_{ij}
- m_i(t): True underlying biomarker trajectory
- b_i\sim\mathcal{N}(0,\sigma_{b}^{2}): Subject-specific random intercept
- \varepsilon_{ij}\sim\mathcal{N}(0,\sigma_{e}^{2}): Measurement error
The biomarker trajectory evolution is characterized by the following second-order differential equation:
\ddot{m}_i(t) = f\big(m_i(t), \dot{m}_i(t), \mathbf{X}_i(t), t\big)
where f: \mathbb{R} \times \mathbb{R} \times \mathbb{R}^p \times \mathbb{R}^+ \to \mathbb{R} is a smooth function modeling the biomarker acceleration as a function of its current value m_i(t), velocity \dot{m}_i(t), time-varying covariates \mathbf{X}_i(t) \in \mathbb{R}^p, and time t.
Survival Sub-Model
The hazard function incorporates biomarker dynamics:
\lambda_i(t) = \lambda_{0}(t)\exp\left[\mathbf{W}_i^{\top}\boldsymbol{\phi}+\mathbf{m}_i(t)^{\top}\boldsymbol{\alpha}+b_{i}\right]
where:
- \lambda_{0}(t): Baseline hazard (e.g., Weibull, piecewise constant)
- \mathbf{m}_i(t)=\left(m_i(t), \dot{m}_i(t), \ddot{m}_i(t)\right)^{\top}: Biomarker value and derivatives
- \boldsymbol{\alpha}=(\alpha_0, \alpha_1, \alpha_2)^{\top}: Association parameters for value, velocity, and acceleration
- \mathbf{W}_i: Baseline covariates with coefficients \boldsymbol{\phi}
- b_i: Subject-specific random intercept
ODE System
The state vector \mathbf{s}_i(t) = (\Lambda_i(t), m_i(t), \dot{m}_i(t))^{\top} evolves according to:
\frac{d\mathbf{s}_i}{dt} = \begin{pmatrix} \lambda_i(t|b_i) \\ \dot{m}_i(t) \\ g(\boldsymbol{\beta}^{\top}\mathbf{Z}_i(t)) \end{pmatrix}
with initial conditions \mathbf{s}_i(0) = (0, m_{i0}, \dot{m}_{i0})^{\top}.
Statistical Inference
Model Components
Hazard Function: \lambda_i(t|b_i) = \exp\left[\boldsymbol{\eta}^{\top} \mathbf{B}^{(\lambda)}(t) + \mathbf{W}_i^{\top}\boldsymbol{\phi} + \mathbf{m}_i(t)^{\top}\boldsymbol{\alpha} + b_{i}\right]
where \mathbf{B}^{(\lambda)}(t) is the B-spline basis for baseline hazard, and \boldsymbol{\eta} are the corresponding coefficients.
Acceleration Function: g(u) = \boldsymbol{\gamma}^{\top} \mathbf{B}^{(g)}(u), \quad u = \boldsymbol{\beta}^{\top}\mathbf{Z}_i(t)
where \mathbf{B}^{(g)}(u) is the B-spline basis for acceleration, \boldsymbol{\gamma} are the basis coefficients, and \boldsymbol{\beta} are the single-index coefficients (constrained: \|\boldsymbol{\beta}\| = 1).
Likelihood
The joint likelihood for subject i integrates over the random effect:
L_i(\boldsymbol{\psi}) = \int f(\mathbf{V}_i | b_i) \cdot f(T_i, \delta_i | b_i) \cdot f(b_i) \, db_i
where \boldsymbol{\psi} = (\boldsymbol{\theta}, \boldsymbol{\beta}), with \boldsymbol{\theta} = (\boldsymbol{\eta}, \boldsymbol{\phi}, \boldsymbol{\alpha}, \boldsymbol{\gamma}, \sigma_e^2, \sigma_b^2).
Likelihood Components:
- Longitudinal: f(\mathbf{V}_i | b_i) = \prod_{j=1}^{n_i} \mathcal{N}(V_{ij}; m_i(T_{ij}) + b_i, \sigma_e^2)
- Survival: f(T_i, \delta_i | b_i) = [\lambda_i(T_i|b_i)]^{\delta_i} \exp[-\Lambda_i(T_i|b_i)]
- Random Effect: f(b_i) \sim \mathcal{N}(0, \sigma_b^2)
EM Algorithm
We use an Expectation-Maximization (EM) algorithm for parameter estimation.
E-Step: Posterior Computation
For each subject i, compute the posterior distribution of b_i given observed data \mathcal{O}_i.
Key simplification: The hazard and cumulative hazard factor as:
- \lambda_i(t|b_i) = e^{b_i} \lambda_i(t|0)
- \Lambda_i(t|b_i) = e^{b_i} \Lambda_i(t|0)
Implementation:
- Solve baseline ODE with b_i = 0 to obtain m_i(t), \lambda_i(t|0), \Lambda_i(T_i|0)
- Find posterior mode \tilde{b}_i by maximizing: \ell_i(b) = b\left[\frac{S_i}{\sigma_e^2} + \delta_i\right] - \frac{b^2}{2}\left[\frac{n_i}{\sigma_e^2} + \frac{1}{\sigma_b^2}\right] - e^b\Lambda_i(T_i|0) where S_i = \sum_j(V_{ij} - m_i(T_{ij}))
-
Compute posterior moments via adaptive
Gauss-Hermite quadrature:
- Mean: \hat{b}_i = E[b_i|\mathcal{O}_i]
- Variance: \hat{v}_i = \text{Var}[b_i|\mathcal{O}_i]
- Transform: E[e^{b_i}|\mathcal{O}_i] for survival updates
M-Step: Parameter Updates
Maximize the expected complete-data log-likelihood:
Q(\boldsymbol{\psi}) = Q_{\text{long}} + Q_{\text{surv}} + Q_{\text{RE}}
where:
- Q_{\text{long}} = -\frac{1}{2\sigma_e^2}\sum_{i,j} [(V_{ij} - m_i(T_{ij}) - \hat{b}_i)^2 + \hat{v}_i] - \frac{N}{2}\log(2\pi\sigma_e^2)
- Q_{\text{surv}} = \sum_i [\delta_i(\log\lambda_i(T_i|0) + \hat{b}_i) - E[e^{b_i}|\mathcal{O}_i]\Lambda_i(T_i|0)]
- Q_{\text{RE}} = -\frac{1}{2\sigma_b^2}\sum_i (\hat{b}_i^2 + \hat{v}_i) - \frac{n}{2}\log(2\pi\sigma_b^2)
Optimization Strategy:
- Update \boldsymbol{\beta}:
\hat{\boldsymbol{\beta}} = \arg\max_{\boldsymbol{\beta}:\|\boldsymbol{\beta}\|=1} Q(\boldsymbol{\beta};\widehat{\boldsymbol{\theta}})
- Update \boldsymbol{\theta}:
\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta};\widehat{\boldsymbol{\psi}}, \hat{\boldsymbol{\gamma}})
- Update \boldsymbol{\sigma}: \sigma_e^2 = \frac{1}{N}\sum_{i,j}[(V_{ij} - m_i(T_{ij}) - \hat{b}_i)^2 + \hat{v}_i],\quad\sigma_b^2 = \frac{1}{n}\sum_i(\hat{b}_i^2 + \hat{v}_i)
Computational Details
Overview of Gradient Computation
The M-step in the EM algorithm requires maximization of the expected complete-data log-likelihood Q(\boldsymbol{\psi}) where \boldsymbol{\psi} = (\boldsymbol{\theta}, \boldsymbol{\beta}). The gradient is:
\nabla_{\boldsymbol{\psi}} Q = \sum_{j=1}^{n_i} \frac{r_{ij}}{\sigma_e^2} \frac{\partial m_i(T_{ij})}{\partial \boldsymbol{\psi}} - E[e^{b_i}|\mathcal{O}_i] \frac{\partial \Lambda_i(T_i)}{\partial \boldsymbol{\psi}}
where r_{ij} = V_{ij} - m_i(T_{ij}) - \hat{b}_i is the residual.
The key computational challenge is obtaining the state sensitivities: - \frac{\partial m_i(T_{ij})}{\partial \boldsymbol{\psi}} at each observation time T_{ij} - \frac{\partial \Lambda_i(T_i)}{\partial \boldsymbol{\psi}} at the event time T_i
Two approaches are available: the forward sensitivity method and the adjoint method.
Sensitivity Analysis via Forward Method
Gradient Computation
Parameter Vector \boldsymbol{\theta}
The gradient components with respect to \boldsymbol{\theta} = (\boldsymbol{\eta}, \boldsymbol{\phi}, \boldsymbol{\alpha}, \boldsymbol{\gamma}) decompose as follows:
Baseline hazard coefficients (\boldsymbol{\eta}): \nabla_{\boldsymbol{\eta}} Q = \sum_{i=1}^{n} \left[\delta_i \mathbf{B}^{(\lambda)}(T_i) - E[e^{b_i}|\mathcal{O}_i] \frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\eta}}\right]
Baseline covariate effects (\boldsymbol{\phi}): \nabla_{\boldsymbol{\phi}} Q = \sum_{i=1}^{n} \left[\delta_i - E[e^{b_i}|\mathcal{O}_i] \cdot \Lambda_i(T_i|0)\right] \mathbf{W}_i
Association parameters (\boldsymbol{\alpha}): \nabla_{\boldsymbol{\alpha}} Q = \sum_{i=1}^{n} \left[\delta_i \mathbf{m}_i(T_i) - E[e^{b_i}|\mathcal{O}_i] \frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\alpha}}\right]
Acceleration spline coefficients (\boldsymbol{\gamma}): \nabla_{\boldsymbol{\gamma}} Q = \sum_{i=1}^{n}\sum_{j=1}^{n_i} \frac{r_{ij}}{\sigma_e^2} \frac{\partial m_i(T_{ij})}{\partial \boldsymbol{\gamma}} + \sum_{i=1}^{n} \left[\delta_i \boldsymbol{\alpha}^{\top} \frac{\partial \mathbf{m}_i(T_i)}{\partial \boldsymbol{\gamma}} - E[e^{b_i}|\mathcal{O}_i] \frac{\partial \Lambda_i(T_i|0)}{\partial \boldsymbol{\gamma}}\right]
where r_{ij} = V_{ij} - m_i(T_{ij}) - \hat{b}_i denotes the residual from the longitudinal model.
Single-Index Coefficients \boldsymbol{\beta}
The gradient with respect to the constrained single-index coefficients (\|\boldsymbol{\beta}\|=1) is:
\nabla_{\boldsymbol{\beta}} Q = \sum_{i=1}^{n}\sum_{j=1}^{n_i} \frac{r_{ij}}{\sigma_e^2} \frac{\partial m_i(T_{ij})}{\partial \boldsymbol{\beta}} + \sum_{i=1}^{n} \left[\delta_i \boldsymbol{\alpha}^{\top} \frac{\partial \mathbf{m}_i(T_i)}{\partial \boldsymbol{\beta}} - E[e^{b_i}|\mathcal{O}_i] \frac{\partial \Lambda_i(T_i|0)}{\partial \boldsymbol{\beta}}\right]
Sensitivity Propagation
Cumulative Hazard Sensitivities
The cumulative hazard function \Lambda_i(T_i|0) depends on parameters through the instantaneous hazard. For parameters affecting the hazard directly, the sensitivities are computed as:
\frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\eta}}=\int_{0}^{T_i}\lambda_i(t|0)\mathbf{B}^{(\lambda)}(t)\,\mathrm{d}t,\quad \frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\alpha}}=\int_{0}^{T_i}\lambda_i(t|0)\mathbf{m}_i(t)\,\mathrm{d}t,\quad \frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\phi}}=\int_{0}^{T_i}\lambda_i(t|0)\mathbf{W}_i\,\mathrm{d}t
For parameters affecting the trajectory m_i(t), the chain rule yields:
\frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\gamma}}=\int_{0}^{T_i}\lambda_i(t|0)\boldsymbol{\alpha}^{\top}\frac{\partial\mathbf{m}_i(t)}{\partial\boldsymbol{\gamma}}\,\mathrm{d}t,\quad\frac{\partial\Lambda_i(T_i|0)}{\partial\boldsymbol{\beta}}=\int_{0}^{T_i}\lambda_i(t|0)\boldsymbol{\alpha}^{\top}\frac{\partial\mathbf{m}_i(t)}{\partial\boldsymbol{\beta}}\,\mathrm{d}t
Trajectory Sensitivities
The biomarker trajectory sensitivities satisfy recursive integral relationships that arise naturally from the hierarchical structure of the second-order differential equation system. These sensitivities must be computed through the forward propagation method due to the coupling between trajectory dynamics and parameter dependencies.
Single-index coefficients (\boldsymbol{\beta}):
The trajectory sensitivities evolve according to:
\frac{\partial m_i(t)}{\partial \boldsymbol{\beta}} = \int_{0}^{t} \frac{\partial\dot{m}_i(s)}{\partial \boldsymbol{\beta}} \,\mathrm{d}s,\quad \frac{\partial\dot{m}_i(t)}{\partial \boldsymbol{\beta}} = \int_{0}^{t} \frac{\partial \ddot{m}_i(s)}{\partial \boldsymbol{\beta}} \,\mathrm{d}s
where the acceleration sensitivity is given by:
\frac{\partial \ddot{m}_i(t)}{\partial \boldsymbol{\beta}} = \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) \frac{\partial u}{\partial\boldsymbol{\beta}}
Here, u = \boldsymbol{\beta}^{\top}\mathbf{Z}_i(t) denotes the single-index value with \mathbf{Z}_i(t) = [m_i(t), \dot{m}_i(t), \mathbf{X}_i^{\top}(t), t]^{\top} being the feature vector. The derivative \mathbf{B}'_g(u) = \frac{d\mathbf{B}_g(u)}{du} represents the derivative of the B-spline basis. The sensitivity of the index is:
\frac{\partial u}{\partial\boldsymbol{\beta}} = \mathbf{Z}_i(t) + \boldsymbol{\beta}^{\top} \frac{\partial \mathbf{Z}_i(t)}{\partial\boldsymbol{\beta}}
where \frac{\partial \mathbf{Z}_i(t)}{\partial\boldsymbol{\beta}} = \begin{bmatrix} \frac{\partial m_i(t)}{\partial \boldsymbol{\beta}} \\ \frac{\partial \dot{m}_i(t)}{\partial \boldsymbol{\beta}} \\ \mathbf{0} \\ 0 \end{bmatrix} since the covariates \mathbf{X}_i(t) and time t do not depend on \boldsymbol{\beta}.
Acceleration spline coefficients (\boldsymbol{\gamma}):
Similarly, the sensitivities with respect to \boldsymbol{\gamma} follow:
\frac{\partial m_i(t)}{\partial \boldsymbol{\gamma}} = \int_{0}^{t} \frac{\partial\dot{m}_i(s)}{\partial \boldsymbol{\gamma}} \,\mathrm{d}s,\quad \frac{\partial\dot{m}_i(t)}{\partial \boldsymbol{\gamma}} = \int_{0}^{t} \frac{\partial \ddot{m}_i(s)}{\partial \boldsymbol{\gamma}} \,\mathrm{d}s
The acceleration sensitivity is:
\frac{\partial \ddot{m}_i(t)}{\partial \boldsymbol{\gamma}} = \mathbf{B}_g(u) + \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) \frac{\partial u}{\partial \boldsymbol{\gamma}}
This expression reveals two distinct contributions: - Direct effect: \mathbf{B}_g(u) captures the explicit linear dependence of the acceleration function on \boldsymbol{\gamma} - Indirect effect: \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) \frac{\partial u}{\partial \boldsymbol{\gamma}} accounts for the implicit dependence through the feedback of trajectory states
The index sensitivity is:
\frac{\partial u}{\partial \boldsymbol{\gamma}} = \boldsymbol{\beta}^{\top} \frac{\partial \mathbf{Z}_i(t)}{\partial \boldsymbol{\gamma}} = \boldsymbol{\beta}^{\top} \begin{bmatrix} \frac{\partial m_i(t)}{\partial \boldsymbol{\gamma}} \\ \frac{\partial \dot{m}_i(t)}{\partial \boldsymbol{\gamma}} \\ \mathbf{0} \\ 0 \end{bmatrix}
This dual dependence structure—explicit through the linear combination and implicit through the state feedback—is a key feature of the single-index acceleration model.
Computational Implementation via Augmented ODE Systems
The practical computation of these gradients necessitates solving augmented ordinary differential equation systems that simultaneously evolve both the primary state variables and their parameter sensitivities. This approach, known as the forward sensitivity method, extends the original three-dimensional state space to include sensitivity trajectories.
Augmented system for \boldsymbol{\theta} components:
\frac{d}{dt}\begin{bmatrix} \Lambda_i \\ m_i \\ \dot{m}_i \\ \partial\Lambda_{\eta,i} \\ \partial\Lambda_{\alpha,i} \\ \partial m_{i,\gamma} \\ \partial\dot{m}_{i,\gamma} \\ \partial\Lambda_{\gamma,i} \end{bmatrix} = \begin{bmatrix} \lambda_i(t|0) \\ \dot{m}_i(t) \\ g(\boldsymbol{\beta}^{\top}\mathbf{Z}_i(t)) \\ \mathbf{B}^{(\lambda)}(t) \lambda_i(t|0) \\ \mathbf{m}_i(t) \lambda_i(t|0) \\ \partial\dot{m}_{i,\gamma} \\ \mathbf{B}_g(u) + \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) \cdot \frac{\partial u}{\partial \boldsymbol{\gamma}} \\ \boldsymbol{\alpha}^{\top}\frac{\partial\mathbf{m}_i(t)}{\partial\boldsymbol{\gamma}} \lambda_i(t|0) \end{bmatrix}
Augmented system for \boldsymbol{\beta} sensitivities:
\frac{d}{dt}\begin{bmatrix} \Lambda_i \\ m_i \\ \dot{m}_i \\ \partial m_{i,\beta} \\ \partial\dot{m}_{i,\beta} \\ \partial\Lambda_{\beta,i} \end{bmatrix} = \begin{bmatrix} \lambda_i(t|0) \\ \dot{m}_i(t) \\ g(\boldsymbol{\beta}^{\top}\mathbf{Z}_i(t)) \\ \partial \dot{m}_{i,\beta} \\ \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) \cdot \frac{\partial u}{\partial\boldsymbol{\beta}} \\ \boldsymbol{\alpha}^{\top}\frac{\partial\mathbf{m}_i(t)}{\partial\boldsymbol{\beta}} \cdot \lambda_i(t|0) \end{bmatrix}
Here, \partial\Lambda_{\eta,i}, \partial\Lambda_{\alpha,i}, \partial\Lambda_{\gamma,i}, and \partial\Lambda_{\beta,i} denote the cumulative hazard sensitivities with respect to the corresponding parameter groups, while \frac{\partial u}{\partial\boldsymbol{\beta}} follows the recursive relationship defined previously. These augmented systems must be integrated from t=0 to t=T_i for each subject to obtain the complete sensitivity information required for gradient computation.
Adjoint Sensitivity Analysis
The adjoint method provides an alternative, more memory-efficient approach to compute the same state sensitivities, especially when the number of parameters is large.
Mathematical Foundation
Given the ODE system \frac{d\mathbf{s}}{dt} = f(t, \mathbf{s}; \boldsymbol{\psi}) with \mathbf{s}(0) = \mathbf{s}_0, we derive the adjoint sensitivity formula.
Since \mathbf{s}(t) satisfies the ODE, for any function \boldsymbol{\kappa}(t): \mathbf{s}(T) = \mathbf{s}(T) - \int_0^T \boldsymbol{\kappa}^{\top} \left[\frac{d\mathbf{s}}{dt} - f\right] dt
Differentiating with respect to \boldsymbol{\psi} and using integration by parts: \frac{\partial \mathbf{s}(T)}{\partial \boldsymbol{\psi}} = \int_0^T \boldsymbol{\kappa}^{\top} \frac{\partial f}{\partial \boldsymbol{\psi}} dt - \boldsymbol{\kappa}(T)^{\top} \frac{\partial \mathbf{s}(T)}{\partial \boldsymbol{\psi}} + \int_0^T \left[\frac{d\boldsymbol{\kappa}}{dt} + \left(\frac{\partial f}{\partial \mathbf{s}}\right)^{\top} \boldsymbol{\kappa}\right]^{\top} \frac{\partial \mathbf{s}}{\partial \boldsymbol{\psi}} dt
Define \tilde{\boldsymbol{\kappa}} = \boldsymbol{\kappa} + \mathbf{e}_k and choose it to satisfy: \frac{d\tilde{\boldsymbol{\kappa}}}{dt} = -\left(\frac{\partial f}{\partial \mathbf{s}}\right)^{\top} \tilde{\boldsymbol{\kappa}}, \quad \tilde{\boldsymbol{\kappa}}(T) = \mathbf{e}_k
This choice eliminates the \frac{\partial \mathbf{s}}{\partial \boldsymbol{\psi}} terms, yielding: \frac{\partial \mathbf{s}_k(T)}{\partial \boldsymbol{\psi}} = \int_0^T \tilde{\boldsymbol{\kappa}}^{\top} \frac{\partial f}{\partial \boldsymbol{\psi}}\bigg|_{\mathbf{s}} dt
where \mathbf{e}_k is the k-th unit vector (i.e., a vector with 1 in the k-th position and 0 elsewhere).
Application to JointODE
For the state vector \mathbf{s} = [\Lambda_i, m_i, \dot{m}_i]^{\top}, we need: 1. \frac{\partial m_i(T_{ij})}{\partial \boldsymbol{\psi}} for each observation time T_{ij} 2. \frac{\partial \Lambda_i(T_i)}{\partial \boldsymbol{\psi}} for the event time T_i
Adjoint System
The transposed Jacobian for the JointODE dynamics is: \left(\frac{\partial f}{\partial \mathbf{s}}\right)^{\top} = \begin{bmatrix} 0 & 0 & 0 \\ \alpha_0 \lambda_i(t|0) & 0 & g'(u) \beta_1 \\ \alpha_1 \lambda_i(t|0) & 1 & g'(u) \beta_2 \end{bmatrix}
where g'(u) = \boldsymbol{\gamma}^{\top}\mathbf{B}'_g(u) denotes the derivative of the acceleration function at u = \boldsymbol{\beta}^{\top}\mathbf{Z}_i(t).
This yields the adjoint differential equations (solved backward in time): \begin{aligned} \frac{d\tilde{\kappa}_1}{dt} &= 0 \\ \frac{d\tilde{\kappa}_2}{dt} &= -\alpha_0 \lambda_i(t|0) \tilde{\kappa}_1 - \tilde{\kappa}_3 \\ \frac{d\tilde{\kappa}_3}{dt} &= -\alpha_1 \lambda_i(t|0) \tilde{\kappa}_1 - g'(u) \beta_1 \tilde{\kappa}_2 - g'(u) \beta_2 \tilde{\kappa}_3 \end{aligned}
with terminal conditions:
- For m_i(T_{ij}) sensitivity: \tilde{\boldsymbol{\kappa}}(T_{ij}) = [0, 1, 0]^{\top}
- For \Lambda_i(T_i) sensitivity: \tilde{\boldsymbol{\kappa}}(T_i) = [1, 0, 0]^{\top}
Parameter Derivatives
The partial derivatives \frac{\partial f}{\partial \boldsymbol{\psi}}\bigg|_{\mathbf{s}} needed for the sensitivity integrals:
Trajectory parameters:
- Single-index coefficients: \frac{\partial f}{\partial \boldsymbol{\beta}} = [0, 0, g'(u) \mathbf{Z}_i(t)]^{\top}
- Acceleration basis coefficients: \frac{\partial f}{\partial \boldsymbol{\gamma}} = [0, 0, \mathbf{B}_g(u)]^{\top}
Survival parameters:
- Baseline hazard coefficients: \frac{\partial f}{\partial \boldsymbol{\eta}} = [\mathbf{B}^{(\lambda)}(t) \lambda_i(t|0), 0, 0]^{\top}
- Association parameters: \frac{\partial f}{\partial \boldsymbol{\alpha}} = [\mathbf{m}_i(t) \lambda_i(t|0), 0, 0]^{\top}
- Baseline covariate effects: \frac{\partial f}{\partial \boldsymbol{\phi}} = [\mathbf{W}_i \lambda_i(t|0), 0, 0]^{\top}