A Framework for the Estimation of Quality Adjusted Life Years Using Joint Models of Longitudinal and Survival Data

This repository includes example data and code supporting the manuscript titled “A Framework for the Estimation of Quality Adjusted Life Years Using Joint Models of Longitudinal and Survival Data”.

This analysis is based on a simulated dataset to ensure that it can be shared freely and without any restriction. The data-generating mechanism is described below, and the Stata code implementing the simulation is included in this repository, in the file named 01-simulate-data.do. The dataset is included in the data folder, and can be downloaded here.

Note that a .pdf version of this README file can be downloaded here, and that all code included in this repository is available under the MIT license.

Data-Generating Mechanism

For this simulation, we generate a cohort of 1,000 patients, randomised 1:1 to a hypothetical treatment arm $X$ (denoted with trt in the data) or control.

In our settings, we assume two longitudinal outcomes (in addition to the time-to-event outcome of all-cause mortality): prostate-specific antigen (PSA), modelled on the log scale, and quality of life (QoL), simulated and then modelled on the logit scale. Once back-transformed to the original scale, QoL is assumed to be bounded between 0 (worst QoL) and 1 (best QoL).

Then, we generate data according to the following trivariate joint model:

$$ \begin{cases} \log \ \text{PSA}(t) &= m_{\text{PSA}}(t) + \varepsilon_{\text{PSA}}(t) \\ \text{logit} \ \text{QoL}(t) &= m_{\text{QoL}}(t) + \varepsilon_{\text{QoL}}(t) \\ h(t) &= h_0(t) \exp \left[ X \gamma + f(\log \text{PSA}(t) | \alpha_{\text{PSA}}) + f(\text{logit} \text{QoL}(t) | \alpha_{\text{QoL}}) \right] \end{cases} $$

where $X$ is the randomised treatment assignment. Each sub-model is described in more detail below.

Note that we also assume:

1. Every patient is observed at baseline;
2. New observations are collected every 6 to 18 months;
3. Administrative censoring at ten years from randomisation (i.e., baseline) is applied.

Model for log PSA

The longitudinal sub-model for log PSA is defined as

$$ \begin{aligned} \log \ \text{PSA}(t) &= m_{\text{PSA}}(t) + \varepsilon_{\text{PSA}}(t) \\ &= (0.0552 + b_{\text{PSA}, 0}) + (1.0204 + b_{\text{PSA}, 1}) \times t + 0.000 \times X - 0.1796 \times t^2 + 0.0101 \times t^3 \\ & \quad \quad - 0.7771 \times X \times t + 0.0937 \times X \times t^2 - 0.0029 \times X \times t^3 + \varepsilon_{\text{PSA}}(t) \end{aligned} $$

where $t$ is time, $X$ is the binary treatment assignment, $b_{\text{PSA}, 0}, b_{\text{PSA}, 1}$ are a random intercept and slope, respectively, and $\varepsilon_{\text{PSA}}(t) \sim N(0, 0.8^2)$ is a residual error term. We draw the random effects for every study subject from the following bivariate normal distribution:

$$ \begin{bmatrix} b_{\text{PSA}, 0} \\ b_{\text{PSA}, 1} \end{bmatrix} \sim N_2 \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 4.343 & -0.598 \\ -0.598 & 1.578 \end{bmatrix} \right) $$

The population-level trajectories, by treatment arm, are as follows:

When applying the inverse transformation, we obtain the following population-level trajectories for PSA values:

Model for logit QoL

The longitudinal sub-model for the logit of QoL is defined as

$$ \begin{aligned} \text{logit} \ \text{QoL}(t) &= m_{\text{QoL}}(t) + \varepsilon_{\text{QoL}}(t) \\ &= (0.07281 + b_{\text{QoL}, 0}) + (-0.29413 + b_{\text{QoL}, 1}) \times t + 0.000 \times X + 0.00734 \times t^2 - 0.00006 \times t^3 \\ & \quad \quad + 0.07979 \times X \times t - 0.00212 \times X \times t^2 + 0.00002 \times X \times t^3 + \varepsilon_{\text{QoL}}(t) \end{aligned} $$

assuming the residual error term $\varepsilon_{\text{QoL}}(t)$ follows a $N(0, 0.2^2)$ distribution and the following random effects distribution:

$$ \begin{bmatrix} b_{\text{QoL}, 0} \\ b_{\text{QoL}, 1} \end{bmatrix} \sim N_2 \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 0.114 & -0.002 \\ -0.002 & 0.013 \end{bmatrix} \right) $$

The population-level trajectories of logit QoL are as follows:

If we apply the inverse transformation, we obtain the longitudinal, population-level trajectories from the following plot:

Model for survival

The time-to-event model for all-cause survival is a Weibull proportional hazards model

$$ \begin{aligned} h(t) &= \lambda p t^{p - 1} \exp \big[ -0.5 \times X + 1.0 \times b_{\text{PSA}, 0} \\ & \qquad - 0.5 \times b_{\text{QoL}, 0} - 0.1 \times b_{\text{QoL}, 1} \big] \end{aligned} $$

with scale parameter $\lambda = 0.05$ and shape parameter $p = 1.5$. As in the longitudinal sub-models, $X$ denotes the randomised treatment assignment with associated log hazard ratio of $-0.5$. Then, we link the longitudinal sub-model with the survival sub-model via the shared random effects association structure with association parameters $\alpha_{\text{PSA}} = 1.0, \alpha_{\text{QoL}, 0} = -0.5, \alpha_{\text{QoL}, 1} = -0.1$ for the random intercept of the PSA submodel and the random intercept and slope of the QoL submodel, respectively.

Data Preview

The dataset includes variable labels describing the content of each column:

. use data/example.dta
(Simulated dataset from a trivariate joint model.)

. describe

Contains data from data/example.dta
 Observations:         6,262                  Simulated dataset from a trivariate joint model.
    Variables:            11                  25 Nov 2025 16:29
----------------------------------------------------------------------------------------------------------------------------
Variable      Storage   Display    Value
    name         type    format    label      Variable label
----------------------------------------------------------------------------------------------------------------------------
id              int     %9.0g                 Unique patient identifier
trt             byte    %9.0g                 Randomised binary treatment assignment (0 = control, 1 = treatment)
stime           double  %10.0g                Survival time
event           byte    %8.0g                 Binary event indicator variable (0 = censored, 1 = event)
n               byte    %9.0g                 Sequential indicator variable for the 1st, 2nd, etc. measurement
N               byte    %9.0g                 Total number of measurements per patient
time            float   %9.0g                 Longitudinal measurement times
log_psa         float   %9.0g                 Longitudinal log PSA values
psa             float   %9.0g                 Longitudinal PSA values
logit_qol       float   %9.0g                 Longitudinal logit quality of life values
qol             float   %9.0g                 Longitudinal quality of life values
----------------------------------------------------------------------------------------------------------------------------
Sorted by: id  time

A short descriptive table of the simulated data is also included below:

. codebook, compact

Variable    Obs Unique       Mean        Min       Max  Label
----------------------------------------------------------------------------------------------------------------------------
id         6262   1000   503.7918          1      1000  Unique patient identifier
trt        6262      2   .5148515          0         1  Randomised binary treatment assignment (0 = control, 1 = treatment)
stime      1000    636   5.699492   .0023133        10  Survival time
event      1000      2       .635          0         1  Binary event indicator variable (0 = censored, 1 = event)
n          6262     14   4.892367          1        14  Sequential indicator variable for the 1st, 2nd, etc. measurement
N          6262     14   8.784733          1        14  Total number of measurements per patient
time       6262   5262   3.805555          0  9.998419  Longitudinal measurement times
log_psa    6262   6262   .7587043  -33.02409  38.96137  Longitudinal log PSA values
psa        6262   6262   1.49e+13   4.55e-15  8.33e+16  Longitudinal PSA values
logit_qol  6262   6262  -.7783099  -5.417343  1.518691  Longitudinal logit quality of life values
qol        6262   6262   .3443611   .0044193  .8203457  Longitudinal quality of life values
----------------------------------------------------------------------------------------------------------------------------