1 Bio-SHiFT Data Analysis

1.1 Summary

This section details the application of the joint modeling framework to a clinical dataset of patients with chronic heart failure. The primary objective is to investigate how the temporal evolution of two biomarkers, eGFR and NGAL, relates to the risk of a composite clinical endpoint. This composite endpoint includes heart failure hospitalization, cardiac death, left ventricular assist device placement, and heart transplantation.

  • Model Specifications: The analysis explores four bivariate joint models built from different combinations of longitudinal and survival submodels. The longitudinal submodels test both linear and nonlinear (natural cubic spline) subject-specific evolutions. The survival submodels test two functional forms: one linking the current biomarker value to the event risk, and another linking the average biomarker value over time to the event risk.

  • Evaluation Methodology: The fit of these four models is evaluated using posterior-posterior predictive checks, posterior-prior checks, and 10-fold cross-validation.

  • Key Findings: Standard information criteria (DIC, WAIC, and LPML) do not uniformly agree on the best model. However, the predictive checks reveal that nonlinear trajectory models better describe the mean function, and that the NGAL model provides a better overall fit than the eGFR model. Furthermore, the area/integral functional form yields a stronger association fit for eGFR.

1.2 Descriptive Statistics

Characteristic N = 2631
AGE 66.7 (12.6; [28.4, 88.0])
SEX
    female 74 (28%)
    male 189 (72%)
NYHA123
    1 77 (29%)
    >1 186 (71%)
DM 81 (31%)
IHD 117 (44%)
Diuretics 237 (90%)
1 Mean (SD; [Min, Max]); n (%)

1.3 Descriptive Figures

  • Sample Means: per longitudinal outcome; the red superimposed line is the LOESS curve.
  • Sample Patients: longitudinal trajectories from 25 randomly selected patients; patient numbers with an asterisk denote patients who experienced the composite endpoint.

1.3.1 Sample Means

1.3.2 Sample Patients

1.3.3 Kaplan Meier Estimate Composite Endpoint

1.4 Parameter Estimates and Information Criteria

  • The four models fitted to the Bio-SHiFT data are
    • jFit_lin_val: linear fixed and random effects structure, current value functional form.
    • jFit_lin_area: linear fixed and random effects structure, area/integral functional form.
    • jFit_nonlin_val: nonlinear (splines) fixed and random effects structure, current value functional form.
    • jFit_nonlin_area: nonlinear (splines) fixed and random effects structure, area/integral functional form.
  • We show parameter estimates and 95% credible intervals for each fitted joint model. For the survival submodel, hazard ratios (the exponent of the corresponding regression coefficients) are presented.
  • Comparison: The Deviance Information Criterion (DIC), the Watanabe-Akaike information criterion (WAIC), and the log pseudo marginal likelihood (LPML) for the fitted joint models. For DIC and WAIC smaller values are better, whereas for LPML larger values are better.

1.4.1 Linear Mixed Models - Value functional Form

Hazard ratios (95% credible intervals) for the survival submodel from the joint models with linear longitudinal submodels and the value functional form.
HR (95% CI)
AGE 0.99 (0.96; 1.02)
SEX 0.97 (0.47; 2.00)
NYHA123 1.97 (1.41; 2.82)
DM 1.52 (0.93; 2.53)
IHD 1.30 (0.78; 2.21)
Diuretics 3.18 (0.82; 20.62)
value(eGFR_CKDEPI_per10) 1.14 (0.80; 1.66)
value(Log2_NGAL_P) 3.04 (1.66; 5.69)
Posterior Means (95% credible intervals) for the eGFR submodel, from the joint models with linear longitudinal submodels and the value functional form.
Estimate (95% CI)
(Intercept) 10.52 (9.50; 11.55)
OBSTIME_YEAR 0.05 (-0.08; 0.18)
AGE -0.07 (-0.08; -0.05)
SEX 0.92 (0.55; 1.29)
NYHA123 -0.01 (-0.24; 0.22)
DM -0.36 (-0.72; -0.01)
IHD 0.42 (0.07; 0.76)
Diuretics -0.90 (-1.44; -0.36)
sigma 1.61 (1.55; 1.67)
Posterior Means (95% credible intervals) for the NGAL submodel, from the joint models with linear longitudinal submodels and the value functional form
Estimate (95% CI)
(Intercept) 5.67 (5.21; 6.12)
OBSTIME_YEAR 0.11 (0.08; 0.14)
AGE 0.02 (0.01; 0.02)
SEX 0.15 (-0.02; 0.31)
NYHA123 0.14 (0.04; 0.24)
DM 0.16 (0.00; 0.32)
IHD -0.10 (-0.26; 0.05)
Diuretics 0.31 (0.06; 0.55)
sigma 0.39 (0.38; 0.41)

1.4.2 Linear Mixed Models - Area functional Form

Hazard ratios (95% credible intervals) for the survival submodel from the joint models with linear longitudinal submodels and the area functional form.
HR (95% CI)
AGE 0.99 (0.96; 1.03)
SEX 0.95 (0.47; 1.90)
NYHA123 1.97 (1.39; 2.83)
DM 1.53 (0.94; 2.53)
IHD 1.29 (0.76; 2.17)
Diuretics 3.08 (0.81; 18.14)
area(eGFR_CKDEPI_per10) 1.16 (0.82; 1.67)
area(Log2_NGAL_P) 3.02 (1.69; 5.62)
Posterior Means (95% credible intervals) for the eGFR submodel, from the joint models with linear longitudinal submodels and the area functional form.
Estimate (95% CI)
(Intercept) 10.51 (9.47; 11.57)
OBSTIME_YEAR 0.05 (-0.08; 0.18)
AGE -0.07 (-0.08; -0.05)
SEX 0.92 (0.53; 1.30)
NYHA123 -0.01 (-0.24; 0.22)
DM -0.36 (-0.73; -0.01)
IHD 0.42 (0.07; 0.77)
Diuretics -0.90 (-1.45; -0.36)
sigma 1.61 (1.55; 1.66)
Posterior Means (95% credible intervals) for the NGAL submodel, from the joint models with linear longitudinal submodels and the area functional form
Estimate (95% CI)
(Intercept) 5.67 (5.21; 6.12)
OBSTIME_YEAR 0.11 (0.08; 0.14)
AGE 0.02 (0.01; 0.02)
SEX 0.15 (-0.01; 0.31)
NYHA123 0.14 (0.04; 0.24)
DM 0.16 (0.01; 0.32)
IHD -0.10 (-0.26; 0.05)
Diuretics 0.31 (0.06; 0.55)
sigma 0.39 (0.38; 0.41)

1.4.3 Nonlinear Mixed Models - Value functional Form

Hazard ratios (95% credible intervals) for the survival submodel from the joint models with nonlinear longitudinal submodels and the value functional form.
HR (95% CI)
AGE 1.00 (0.97; 1.03)
SEX 0.90 (0.44; 1.85)
NYHA123 1.98 (1.39; 2.82)
DM 1.54 (0.94; 2.55)
IHD 1.29 (0.77; 2.17)
Diuretics 3.15 (0.83; 18.05)
value(eGFR_CKDEPI_per10) 1.21 (0.88; 1.71)
value(Log2_NGAL_P) 3.12 (1.74; 5.68)
Posterior Means (95% credible intervals) for the eGFR submodel, from the joint models with nonlinear longitudinal submodels and the value functional form.
Estimate (95% CI)
(Intercept) 10.40 (9.39; 11.44)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))1 0.59 (0.17; 1.02)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))2 -0.46 (-1.03; 0.11)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))3 -0.61 (-1.31; 0.09)
AGE -0.07 (-0.08; -0.05)
SEX 0.95 (0.56; 1.32)
NYHA123 0.02 (-0.22; 0.24)
DM -0.37 (-0.73; -0.01)
IHD 0.37 (0.01; 0.71)
Diuretics -0.85 (-1.39; -0.31)
sigma 1.58 (1.52; 1.64)
Posterior Means (95% credible intervals) for the NGAL submodel, from the joint models with nonlinear longitudinal submodels and the value functional form
Estimate (95% CI)
(Intercept) 5.72 (5.28; 6.16)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))1 0.06 (-0.05; 0.15)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))2 0.37 (0.14; 0.47)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))3 0.54 (0.34; 0.69)
AGE 0.02 (0.01; 0.02)
SEX 0.15 (-0.02; 0.31)
NYHA123 0.13 (0.04; 0.23)
DM 0.17 (0.01; 0.32)
IHD -0.10 (-0.25; 0.06)
Diuretics 0.30 (0.06; 0.55)
sigma 0.38 (0.37; 0.40)

1.4.4 Nonlinear Mixed Models - Area functional Form

Hazard ratios (95% credible intervals) for the survival submodel from the joint models with nonlinear longitudinal submodels and the area functional form.
HR (95% CI)
AGE 0.99 (0.96; 1.02)
SEX 1.00 (0.50; 2.00)
NYHA123 1.98 (1.39; 2.80)
DM 1.51 (0.92; 2.48)
IHD 1.32 (0.78; 2.19)
Diuretics 2.98 (0.79; 17.12)
area(eGFR_CKDEPI_per10) 1.11 (0.80; 1.54)
area(Log2_NGAL_P) 2.81 (1.60; 5.03)
Posterior Means (95% credible intervals) for the eGFR submodel, from the joint models with nonlinear longitudinal submodels and the area functional form.
Estimate (95% CI)
(Intercept) 10.42 (9.39; 11.48)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))1 0.58 (0.14; 1.00)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))2 -0.48 (-1.00; 0.08)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))3 -0.62 (-1.31; 0.06)
AGE -0.07 (-0.08; -0.05)
SEX 0.95 (0.57; 1.31)
NYHA123 0.01 (-0.22; 0.24)
DM -0.37 (-0.73; -0.01)
IHD 0.37 (0.02; 0.72)
Diuretics -0.85 (-1.39; -0.30)
sigma 1.58 (1.52; 1.64)
Posterior Means (95% credible intervals) for the NGAL submodel, from the joint models with nonlinear longitudinal submodels and the area functional form
Estimate (95% CI)
(Intercept) 5.72 (5.27; 6.16)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))1 0.05 (-0.05; 0.15)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))2 0.39 (0.20; 0.49)
ns(OBSTIME_YEAR, 3, B = c(0, 3.2))3 0.55 (0.37; 0.71)
AGE 0.02 (0.01; 0.02)
SEX 0.15 (-0.01; 0.31)
NYHA123 0.13 (0.04; 0.23)
DM 0.17 (0.01; 0.33)
IHD -0.10 (-0.25; 0.06)
Diuretics 0.30 (0.06; 0.55)
sigma 0.38 (0.37; 0.40)

1.4.5 Comparison

Deviance information criterion (DIC), Watanabe-Akaike information criterion (WAIC), and log pseudo marginal likelihood (LPML) for the four fitted joint models.
Model DIC WAIC LPML
4 Nonlinear-Area 10765.51 10995.24 -5526.914
3 Nonlinear-Value 10756.25 11013.66 -5550.674
2 Linear-Area 11120.21 11113.62 -5607.279
1 Linear-Value 11117.45 11116.76 -5601.697

1.5 Longitudinal Models: Posterior-Posterior Checks

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.5.1 ECDF

1.5.2 Mean Function

1.5.3 Variance Function

1.5.4 Correlation Structure

1.6 Longitudinal Models: Posterior-Posterior Checks - transformed eGFR

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions from the joint model with a mixed-effects model for the square root transformation of eGFR with Student’s t error terms with four degrees of freedom. For NGAL the linear mixed model was used. In both models, the subject-specific profiles are modeled with natural cubic splines. The value functional form is use the survival submodel.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.6.1 ECDF

1.6.2 Mean Function

1.6.3 Variance Function

1.6.4 Correlation Structure

1.7 Longitudinal Models: Posterior-Prior Checks

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.7.1 ECDF

1.7.2 Mean Function

1.7.3 Variance Function

1.7.4 Correlation Structure

1.8 Longitudinal Models: Cross-Validated Posterior-Prior Checks

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The simulated datasets have been generated using the cross-validation procedure described in the main paper using \(V = 10\) folds.

1.8.1 ECDF

1.8.2 Mean Function

1.8.3 Variance Function

1.8.4 Correlation Structure

1.9 Longitudinal Models: Cross-Validated Dynamic-Posterior-Posterior Checks at landmark time t_L = 1

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines) using longitudinal information up to the landmark time \(t_L = 1\) year.
  • The simulated datasets have been generated using the cross-validation procedure described in the main paper using \(V = 10\) folds.

1.9.1 ECDF

1.9.2 Mean Function

1.9.3 Variance Function

1.9.4 Correlation Structure

1.10 Longitudinal Models: Cross-Validated Dynamic-Posterior-Posterior Checks at landmark time t_L = 1.5

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines) using longitudinal information up to the landmark time \(t_L = 1.5\) years.
  • The simulated datasets have been generated using the cross-validation procedure described in the main paper using \(V = 10\) folds.

1.10.1 ECDF

1.10.2 Mean Function

1.10.3 Variance Function

1.10.4 Correlation Structure

1.11 Longitudinal Models: Cross-Validated Dynamic-Posterior-Posterior Checks at landmark time t_L = 2

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines) using longitudinal information up to the landmark time \(t_L = 2\) years.
  • The simulated datasets have been generated using the cross-validation procedure described in the main paper using \(V = 10\) folds.

1.11.1 ECDF

1.11.2 Mean Function

1.11.3 Variance Function

1.11.4 Correlation Structure

1.12 Longitudinal Models: Individualized Posterior-Posterior Checks

  • Mean function for the two longitudinal outcomes and Patients 180 and 124 based on the two fitted joint models using the current value functional form. The results from the area/integral functional form were very similar. The asterisks denote the observed data points.
  • Patient 124 experienced the composite endpoint, whereas Patient 180 did not.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.12.1 Patient 180

1.12.2 Patient 124

1.13 Event Time Model: Posterior-Posterior Checks

  • The two tabs show the empirical cumulative distribution function and the probability integral transform of the subject-specific cumulative distribution functions for the four fitted joint models.
  • In the title of each figure, ‘Linear/Nonlinear’ refers to the specification of the linear mixed models (i.e., linear time trends or splines in both the fixed- and random-effects), and ‘Value/Area’ refers to the functional form (i.e., current value or integral/area).
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.13.1 ECDF

1.13.2 Survival Transform to Uniform

1.14 Event Time Model: Posterior-Prior Checks

  • The two tabs show the empirical cumulative distribution function and the probability integral transform of the subject-specific cumulative distribution functions for the four fitted joint models.
  • In the title of each figure, ‘Linear/Nonlinear’ refers to the specification of the linear mixed models (i.e., linear time trends or splines in both the fixed- and random-effects), and ‘Value/Area’ refers to the functional form (i.e., current value or integral/area).
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.14.1 ECDF

1.14.2 Survival Transform to Uniform

1.15 Association: Posterior-Posterior Checks

  • The two tabs show the concordance statistic for the two longitudinal outcomes based on the four fitted joint models.
  • In the title of each figure, ‘Linear/Nonlinear’ refers to the specification of the linear mixed models (i.e., linear time trends or splines in both the fixed- and random-effects), and ‘Value/Area’ refers to the functional form (i.e., current value or integral/area).
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).

1.15.1 eGFR

1.15.2 NGAL

2 Simulation

2.1 Summary

This section serves as a controlled proof-of-concept to demonstrate how the proposed posterior predictive checks perform when the true data-generating mechanism is known.

  • Design: Data for 300 subjects was simulated using an explicitly defined true model. This included a linear mixed-effects model with natural cubic splines for time, and a hazard model utilizing the current value functional form.

  • Model Testing: The simulated data is analyzed using the true model alongside three intentionally misspecified models. The misspecified models include one assuming a strict linear time trend, one transforming the longitudinal outcome to \(y^* = \exp(y)\), and one using the current slope or velocity for the survival functional form.

  • Key Findings: Standard information criteria incorrectly selected the misspecified current-slope model as the best fit. In contrast, the proposed predictive checks successfully identified the exact sources of misspecification. For instance, the empirical cumulative distribution function (eCDF) caught the transformation error, while the mean, variance, and correlation functions identified the incorrect linear time trend.

2.2 Model Comparison

  • The four models fitted to the simulated datasets are
    • jmFit_true: the true model from which the dataset was simulated.
    • jmFit_linear: the joint model with a misspecified the longitudinal submodel by assuming linear subject-specific trends.
    • jmFit_exp: the joint model with a misspecified the longitudinal submodel by fitting the mixed model for the outcome variable \(y^* = \exp(y)\).
    • jmFit_slope: the joint model with a misspecified the survival submodel using the current slope/velocity as the functional form.
  • Comparison: The Deviance Information Criterion (DIC), the Watanabe-Akaike information criterion (WAIC), and the log pseudo marginal likelihood (LPML) for the fitted joint models. For DIC and WAIC smaller values are better, whereas for LPML larger values are better.
## 
##                      DIC      WAIC      LPML
##   jmFit1_slope  2818.188  5997.565 -6506.728
##    jmFit1_true  2815.217  6035.847 -5950.646
##  jmFit1_linear 10680.521 11154.493 -5679.796
##     jmFit1_exp 14727.852 15691.203 -7858.881
## 
## The criteria are calculated on the basis of the marginal log-likelihood.
## 
##                      DIC      WAIC      LPML
##    jmFit2_true  3010.890  6694.725 -5374.090
##   jmFit2_slope  3032.507  6772.072 -5739.441
##  jmFit2_linear 10422.728 10480.900 -5251.739
##     jmFit2_exp 17681.645 25554.386 -9979.662
## 
## The criteria are calculated on the basis of the marginal log-likelihood.
## 
##                      DIC      WAIC      LPML
##    jmFit3_true  2240.666  4844.627 -5961.789
##   jmFit3_slope  2268.857  5005.813 -5622.762
##  jmFit3_linear 10346.653 10382.342 -5191.319
##     jmFit3_exp 16763.732 26213.332 -9690.560
## 
## The criteria are calculated on the basis of the marginal log-likelihood.

2.3 Longitudinal Model - Simulated Data 1

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the four fitted joint models.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.3.1 ECDF

2.3.2 Mean Function

2.3.3 Variance Function

2.3.4 Correlation Structure

2.4 Longitudinal Model - Simulated Data 2

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the four fitted joint models.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.4.1 ECDF

2.4.2 Mean Function

2.4.3 Variance Function

2.4.4 Correlation Structure

2.5 Longitudinal Model - Simulated Data 3

  • The four tabs show the empirical cumulative distribution, mean, variance, and semi-variogram functions for the four fitted joint models.
  • In each figure and for the two longitudinal outcomes, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.5.1 ECDF

2.5.2 Mean Function

2.5.3 Variance Function

2.5.4 Correlation Structure

2.6 Event Time Model - Simulated Data 1

  • The two tabs show the empirical cumulative distribution function and the probability integral transform of the subject-specific cumulative distribution functions for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.6.1 ECDF

2.6.2 Survival Transform to Uniform

2.7 Event Time Model - Simulated Data 2

  • The two tabs show the empirical cumulative distribution function and the probability integral transform of the subject-specific cumulative distribution functions for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.7.1 ECDF

2.7.2 Survival Transform to Uniform

2.8 Event Time Model - Simulated Data 3

  • The two tabs show the empirical cumulative distribution function and the probability integral transform of the subject-specific cumulative distribution functions for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.8.1 ECDF

2.8.2 Survival Transform to Uniform

2.9 Association - Simulated Data 1

  • The concordance statistic for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.10 Association - Simulated Data 2

  • The concordance statistic for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.

2.11 Association - Simulated Data 3

  • The concordance statistic for the four fitted joint models.
  • In each figure, we compare the observed data metric (black line) with the same metric calculated in 50 simulated datasets from the corresponding model (grey lines).
  • The posterior-posterior checks are shown.