Introduction to RetroJM

Overview

RetroJM fits a two-stage retrospective joint model for longitudinal and survival data. The main idea is to analyze the longitudinal outcome on a retrospective timescale, where time is measured backward from death rather than forward from baseline.

This perspective is useful when the scientific question concerns how a longitudinal outcome, such as quality of life, changes as a subject approaches death.

The longitudinal model is: \[ Y_i(t^*)=\beta_{\mu}(t^*)+A_i\beta_A(t^*)+X_i^T\psi_X+W_i b_i+\epsilon_i(t^*) \] where

\(t^*=D-t\), \(D\) is death time and \(t\) is the prospective observation time.

The survival component is modeled using a Cox proportional hazards model. The hazard for death: \[ \lambda_i(t)=\exp{\Big(\alpha_0(t)+A_i\alpha_A(t)\Big)}, \]

For subjects with observed death times, the retrospective timescale is directly available. For censored subjects, death times are not observed and are multiply imputed from the fitted survival model. The longitudinal model is then fit on each imputed dataset, and the final estimates are pooled across imputations.

Package installation

To install the released version from CRAN, run:

install.packages("TTM")

To install the development version from GitHub, run:

devtools::install_github("TChang-W/TTM")

Main function: RetroJM()

The input data should be in long format, with repeated longitudinal measurements for each subject. The dataset should contain subject identifiers, longitudinal measurements, observation times, treatment assignment, baseline covariates if applicable, event or censoring times, and the death indicator.

Arguments

The most important arguments are:

  • data_long: The dataset used in analysis in long format. It should contain longitudinal measurement, measurement time, ID for subjects, covariates of interest, event/censoring time, an indicator for death status, and random effects by specification.

  • id: Name of ID variable.

  • treatment: Name of parameter of interest. Have to be a binary variable coded as 0/1, such as control/treatment.

  • outcome: Name of longitudinal measurement.

  • covariates: Name of covariates. Default is NULL, meaning no covariates.

  • os_time: Name of variable for event/censoring time.

  • death: Name of indicator variable for death status.

  • t_pros: Name of observed time variables. If the input data is long format, it should be a single character; If the input is wide format, it should be a character vector. It should have same length as “outcome”.

  • random_effects: Names of the random effects to be controlled.

  • n_impute: The number of imputations performed per censored subject.

  • n_inner_knot: The number of inner knots used to approximate the trajectories beta_mu and beta_A

Examples

Dataset

The example datasets Retro_data is included in the package. They could be accessed by:

library(TTM)

data(Retro_data)

The Retro_data is a long format data set. It contains 1000 subjects with repeated measurements of one longitudinal outcome.

From the data, the name of longitudinal measurement is: Y, name of measured time are: t_pros, name of treatment assignment is: A, name of event/censoring time is OS_time, name of covariates are: X1, X2, name of ID variable is: id, name of death status is: death. We want to include a random intercept and a random slope for A:

Fitting the model

fit <- RetroJM(
   data_long = Retro_data,
   id = "id",
   treatment = "A",
   outcome = "Y",
   covariates = c("X1", "X2"),
   os_time = "OS_time",
   death = "death",
   t_pros = "t_pros",
   random_effects = c("1", "A"),
   n_impute = 10,
   n_inner_knot = 3
 )

Output

The fitted object returned by RetroJM() contains the fitted survival model, the imputed death times, the spline specification, and the pooled longitudinal model results.

Pooled results

The pooled longitudinal results can be accessed by:

fit$pooled_longitudinal_model
#> $fixed_effects
#> (Intercept)           A          X1          X2   X_spline1   X_spline2 
#>   2.1972141   0.8999984   8.4422357  18.6877301   8.2885672  15.0676142 
#>   X_spline3   X_spline4  AX_spline1  AX_spline2  AX_spline3  AX_spline4 
#>  23.7689848  11.5559050   2.4424850   7.7814608  -2.8707962 -12.4016935 
#> 
#> $fixed_effects_vcov
#> 12 x 12 Matrix of class "dsyMatrix"
#>             (Intercept)           A            X1            X2    X_spline1
#> (Intercept)  0.54008825 -0.44150088 -0.0734370415 -0.0542954146  -0.31926209
#> A           -0.44150088  0.90070404 -0.0126341167 -0.0152567790   0.37389032
#> X1          -0.07343704 -0.01263412  0.1874373451 -0.0002298076  -0.02474175
#> X2          -0.05429541 -0.01525678 -0.0002298076  0.0514298588  -0.02722257
#> X_spline1   -0.31926209  0.37389032 -0.0247417480 -0.0272225670   0.94286866
#> X_spline2   -0.76698192  0.91027693 -0.3418107489  0.1808599897  -4.17094650
#> X_spline3   -0.51179328 -3.42379184  2.7728759472  0.2398744402  14.56870979
#> X_spline4    1.66948902 -9.52965217  5.7373220739  0.3031327069  31.64148381
#> AX_spline1   0.33434387 -0.68505528 -0.0236459937  0.0124391916  -0.95137290
#> AX_spline2   0.59728218 -1.77797060  0.4955423522  0.1207033015   3.51107037
#> AX_spline3   0.76618927  2.28804467 -0.5319054623 -0.7406781322  -7.22480090
#> AX_spline4  -1.03596216 10.01228582 -1.6432440944 -1.6109237691 -16.72346493
#>                X_spline2     X_spline3     X_spline4   AX_spline1    AX_spline2
#> (Intercept)   -0.7669819    -0.5117933  1.669489e+00   0.33434387     0.5972822
#> A              0.9102769    -3.4237918 -9.529652e+00  -0.68505528    -1.7779706
#> X1            -0.3418107     2.7728759  5.737322e+00  -0.02364599     0.4955424
#> X2             0.1808600     0.2398744  3.031327e-01   0.01243919     0.1207033
#> X_spline1     -4.1709465    14.5687098  3.164148e+01  -0.95137290     3.5110704
#> X_spline2     95.3693722  -444.1566714 -9.596002e+02   5.65045415   -94.2971287
#> X_spline3   -444.1566714  3274.5609345  6.795497e+03 -22.99077495   471.0526285
#> X_spline4   -959.6002109  6795.4966483  1.415715e+04 -49.37085596  1009.7159968
#> AX_spline1     5.6504542   -22.9907749 -4.937086e+01   1.95392633    -9.3252646
#> AX_spline2   -94.2971287   471.0526285  1.009716e+03  -9.32526463   175.6074220
#> AX_spline3   252.6415463 -2663.8017831 -5.423529e+03  19.09172707  -514.8648591
#> AX_spline4   591.9560118 -5681.6023816 -1.163666e+04  44.94125527 -1185.3490672
#>                AX_spline3    AX_spline4
#> (Intercept)     0.7661893     -1.035962
#> A               2.2880447     10.012286
#> X1             -0.5319055     -1.643244
#> X2             -0.7406781     -1.610924
#> X_spline1      -7.2248009    -16.723465
#> X_spline2     252.6415463    591.956012
#> X_spline3   -2663.8017831  -5681.602382
#> X_spline4   -5423.5291970 -11636.657582
#> AX_spline1     19.0917271     44.941255
#> AX_spline2   -514.8648591  -1185.349067
#> AX_spline3   4449.6594354   9187.000798
#> AX_spline4   9187.0007983  19125.186150
#> 
#> $residual_variance
#> [1] 103.6336
#> 
#> $random_effects_covariance
#> [1] 4.512448

Baseline trajectory

To plot the estimated baseline trajectory \(\beta_\mu(t*)\) with pointwise confidence intervals:

plot_beta_mu(fit, xlim = c(0, quantile(Retro_data$OS_time, 0.7)), conf.level = 0.95)

To plot the estimated baseline trajectory \(\beta_A(t*)\) with pointwise confidence intervals:

plot_beta_A(fit, xlim = c(0, quantile(Retro_data$OS_time, 0.7)), conf.level = 0.95)

Extract trajectory estimates at specific time points

To evaluate the baseline trajectory \(\beta_\mu(t*)\) at selected retrospective times:

get_beta_mu(new_times = 1:10, results = fit, conf.level = 0.95)
#> $time_grid
#>  [1]  1  2  3  4  5  6  7  8  9 10
#> 
#> $estimate
#>  [1] 2.413251 2.849287 3.283694 3.715381 4.143257 4.566234 4.983219 5.393124
#>  [9] 5.794858 6.187330
#> 
#> $standard_error
#>  [1] 0.7075906 0.6567125 0.6128251 0.5772737 0.5510883 0.5347333 0.5279165
#>  [8] 0.5295613 0.5379745 0.5511313
#> 
#> $CI.lower
#>  [1] 1.026399 1.562154 2.082579 2.583945 3.063144 3.518176 3.948522 4.355203
#>  [9] 4.740447 5.107133
#> 
#> $CI.upper
#>  [1] 3.800103 4.136420 4.484809 4.846817 5.223371 5.614292 6.017916 6.431045
#>  [9] 6.849268 7.267528

To evaluate the baseline trajectory \(\beta_A(t*)\) at selected retrospective times:

get_beta_A(new_times = 1:10, results = fit, conf.level = 0.95)
#> $time_grid
#>  [1]  1  2  3  4  5  6  7  8  9 10
#> 
#> $estimate
#>  [1] 0.8945286 0.8838421 0.8742626 0.8665306 0.8613869 0.8595718 0.8618261
#>  [8] 0.8688902 0.8815047 0.9004102
#> 
#> $standard_error
#>  [1] 0.9065532 0.8267222 0.7571060 0.7002768 0.6584760 0.6330203 0.6237252
#>  [8] 0.6287151 0.6447995 0.6681865
#> 
#> $CI.lower
#>  [1] -0.8822830 -0.7365037 -0.6096380 -0.5059867 -0.4292025 -0.3811251
#>  [7] -0.3606528 -0.3633688 -0.3822791 -0.4092113
#> 
#> $CI.upper
#>  [1] 2.671340 2.504188 2.358163 2.239048 2.151976 2.100269 2.084305 2.101149
#>  [9] 2.145288 2.210032

These functions return the estimated trajectory, pointwise standard error, and pointwise confidence interval at each requested time point.

Interpreting the trajectories

The function \(\beta_\mu(t*)\) represents the estimated baseline retrospective trajectory, while \(\beta_A(t*)\) represents the estimated treatment effect as a function of retrospective time.

These curves help describe how the longitudinal outcome changes as subjects approach death, and how that pattern differs by treatment group.

Summary

RetroJM provides a practical framework for retrospective joint modeling when death times are missing for censored subjects. By combining survival-based imputation with mixed-effects modeling, the package allows users to estimate and visualize longitudinal trajectories on the timescale of time before death