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.1041130   0.9906772   8.5050098  18.6916119   8.0662918  14.7108156 
#>   X_spline3   X_spline4  AX_spline1  AX_spline2  AX_spline3  AX_spline4 
#>  15.7106057  -4.3123072   2.6708466   2.6560630  20.8271074  39.7151475 
#> 
#> $fixed_effects_vcov
#> 12 x 12 Matrix of class "dsyMatrix"
#>             (Intercept)            A           X1           X2   X_spline1
#> (Intercept)  0.57338869 -0.466627745 -0.078730075 -0.053956812 -0.29015082
#> A           -0.46662775  0.936020613 -0.008781805 -0.006307536  0.36436795
#> X1          -0.07873007 -0.008781805  0.190849010  0.001004553 -0.03148702
#> X2          -0.05395681 -0.006307536  0.001004553  0.048436019 -0.01287468
#> X_spline1   -0.29015082  0.364367949 -0.031487017 -0.012874675  0.91216824
#> X_spline2   -0.99182041  0.454598753  0.026968853  0.053049960 -2.41839230
#> X_spline3   -0.42402250  0.806831156  0.494819879  0.153136816  2.95097141
#> X_spline4    2.32349464 -0.985562378  0.938907927  0.244149801  7.72412426
#> AX_spline1   0.31177909 -0.670371673 -0.014942818 -0.008222001 -0.94832920
#> AX_spline2   0.69045389 -1.517104724  0.184841317  0.144208877  1.97269730
#> AX_spline3   2.49392766  0.324407676 -0.877509889 -0.260937782  2.95956790
#> AX_spline4   1.98639839  6.167324309 -1.971490320 -0.641472843  4.40877269
#>                 X_spline2    X_spline3    X_spline4   AX_spline1   AX_spline2
#> (Intercept)   -0.99182041   -0.4240225    2.3234946  0.311779085    0.6904539
#> A              0.45459875    0.8068312   -0.9855624 -0.670371673   -1.5171047
#> X1             0.02696885    0.4948199    0.9389079 -0.014942818    0.1848413
#> X2             0.05304996    0.1531368    0.2441498 -0.008222001    0.1442089
#> X_spline1     -2.41839230    2.9509714    7.7241243 -0.948329205    1.9726973
#> X_spline2     32.68825328  -44.8208939 -125.4151283  2.874907328  -29.9346193
#> X_spline3    -44.82089387  201.8816071  432.7209654 -2.756970184   32.8801597
#> X_spline4   -125.41512835  432.7209654  979.5676397 -8.030669839   99.2403027
#> AX_spline1     2.87490733   -2.7569702   -8.0306698  2.003329373   -5.6610058
#> AX_spline2   -29.93461934   32.8801597   99.2403027 -5.661005757   71.4605799
#> AX_spline3   -11.87287401  -86.1458820 -130.5438748  4.313303342 -110.5845504
#> AX_spline4    10.49885291 -192.7099968 -358.1177468 12.891373369 -293.6549081
#>               AX_spline3   AX_spline4
#> (Intercept)    2.4939277    1.9863984
#> A              0.3244077    6.1673243
#> X1            -0.8775099   -1.9714903
#> X2            -0.2609378   -0.6414728
#> X_spline1      2.9595679    4.4087727
#> X_spline2    -11.8728740   10.4988529
#> X_spline3    -86.1458820 -192.7099968
#> X_spline4   -130.5438748 -358.1177468
#> AX_spline1     4.3133033   12.8913734
#> AX_spline2  -110.5845504 -293.6549081
#> AX_spline3   812.3717396 1668.1721319
#> AX_spline4  1668.1721319 3532.5260608
#> 
#> $residual_variance
#> [1] 103.695
#> 
#> $random_effects_covariance
#> [1] 4.840688

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.400674 2.849023 3.295309 3.738291 4.176726 4.609371 5.034985 5.452325
#>  [9] 5.860149 6.257214
#> 
#> $standard_error
#>  [1] 0.7168270 0.6612098 0.6138815 0.5765042 0.5503106 0.5357460 0.5322212
#>  [8] 0.5381386 0.5511982 0.5688176
#> 
#> $CI.lower
#>  [1] 0.9957189 1.5530756 2.0921237 2.6083636 3.0981368 3.5593282 3.9918504
#>  [8] 4.3975924 4.7798200 5.1423520
#> 
#> $CI.upper
#>  [1] 3.805629 4.144970 4.498495 4.868218 5.255315 5.659414 6.078119 6.507057
#>  [9] 6.940477 7.372076

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.9594689 0.9128807 0.8680434 0.8260117 0.7878401 0.7545832 0.7272958
#>  [8] 0.7070323 0.6948473 0.6917955
#> 
#> $standard_error
#>  [1] 0.9082543 0.8257111 0.7542294 0.6966069 0.6552445 0.6314461 0.6247760
#>  [8] 0.6329207 0.6522129 0.6784763
#> 
#> $CI.lower
#>  [1] -0.8206767 -0.7054833 -0.6102189 -0.5393128 -0.4964156 -0.4830284
#>  [7] -0.4972428 -0.5334695 -0.5834665 -0.6379936
#> 
#> $CI.upper
#>  [1] 2.739615 2.531245 2.346306 2.191336 2.072096 1.992195 1.951834 1.947534
#>  [9] 1.973161 2.021585

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.

Introduction to RetroJM