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.2379358   0.9502794   8.4775547  18.7046724   7.8569400  13.8507494 
#>   X_spline3   X_spline4  AX_spline1  AX_spline2  AX_spline3  AX_spline4 
#>  17.7277904   1.1318662   2.9829673   0.4866215  13.1230378  26.4769189 
#> 
#> $fixed_effects_vcov
#> 12 x 12 Matrix of class "dsyMatrix"
#>             (Intercept)            A           X1           X2   X_spline1
#> (Intercept)  0.52801144 -0.438975901 -0.073205355 -0.054190751 -0.30655957
#> A           -0.43897590  0.878109639 -0.010557121 -0.006464373  0.35264758
#> X1          -0.07320536 -0.010557121  0.192595732 -0.001605569 -0.04386407
#> X2          -0.05419075 -0.006464373 -0.001605569  0.049404047 -0.02265760
#> X_spline1   -0.30655957  0.352647576 -0.043864073 -0.022657595  0.93740268
#> X_spline2   -0.39977367  0.052345973  0.165188974  0.113323656 -1.80286469
#> X_spline3   -1.13098457  1.498361452 -0.089157894  0.085953162  1.38840621
#> X_spline4    0.12912211  0.733485773 -0.273892932  0.047343819  3.83256305
#> AX_spline1   0.35931105 -0.725229052  0.002079987  0.001515364 -0.93019050
#> AX_spline2   0.11307262 -0.264012719 -0.063753185  0.064522465  1.68375000
#> AX_spline3   0.58095462 -1.254562648  0.505083281 -0.187122075 -1.23726840
#> AX_spline4  -0.99916341  1.675213898  0.885586401 -0.407033541 -3.55036478
#>                X_spline2     X_spline3     X_spline4   AX_spline1    AX_spline2
#> (Intercept)  -0.39977367   -1.13098457    0.12912211  0.359311052    0.11307262
#> A             0.05234597    1.49836145    0.73348577 -0.725229052   -0.26401272
#> X1            0.16518897   -0.08915789   -0.27389293  0.002079987   -0.06375319
#> X2            0.11332366    0.08595316    0.04734382  0.001515364    0.06452246
#> X_spline1    -1.80286469    1.38840621    3.83256305 -0.930190498    1.68375000
#> X_spline2    22.06042020  -35.97869079  -90.20514968  1.648529846  -19.58869525
#> X_spline3   -35.97869079  180.34830460  366.25838170 -0.856549563   27.36597629
#> X_spline4   -90.20514968  366.25838170  777.41904830 -2.538573722   71.16228583
#> AX_spline1    1.64852985   -0.85654956   -2.53857372  1.965613632   -3.71263852
#> AX_spline2  -19.58869525   27.36597629   71.16228583 -3.712638518   38.69285258
#> AX_spline3   28.46588853 -165.41624729 -338.76757242  1.176740407  -38.59308803
#> AX_spline4   72.64777389 -326.06451312 -701.23481234  4.677940399 -111.71226641
#>               AX_spline3   AX_spline4
#> (Intercept)    0.5809546   -0.9991634
#> A             -1.2545626    1.6752139
#> X1             0.5050833    0.8855864
#> X2            -0.1871221   -0.4070335
#> X_spline1     -1.2372684   -3.5503648
#> X_spline2     28.4658885   72.6477739
#> X_spline3   -165.4162473 -326.0645131
#> X_spline4   -338.7675724 -701.2348123
#> AX_spline1     1.1767404    4.6779404
#> AX_spline2   -38.5930880 -111.7122664
#> AX_spline3   321.7446305  641.5283694
#> AX_spline4   641.5283694 1346.5859863
#> 
#> $residual_variance
#> [1] 103.6262
#> 
#> $random_effects_covariance
#> [1] 4.489927

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.395930 2.835564 3.273521 3.708569 4.139473 4.564999 4.983915 5.394986
#>  [9] 5.796979 6.188661
#> 
#> $standard_error
#>  [1] 0.7062446 0.6537787 0.6090359 0.5735593 0.5485056 0.5343317 0.5305742
#>  [8] 0.5358575 0.5481438 0.5650970
#> 
#> $CI.lower
#>  [1] 1.011716 1.554181 2.079833 2.584413 3.064422 3.517729 3.944009 4.344725
#>  [9] 4.722637 5.081091
#> 
#> $CI.upper
#>  [1] 3.780144 4.116946 4.467210 4.832725 5.214524 5.612270 6.023821 6.445248
#>  [9] 6.871321 7.296230

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.9346453 0.8915381 0.8499800 0.8111107 0.7760701 0.7459981 0.7220345
#>  [8] 0.7053192 0.6969919 0.6981927
#> 
#> $standard_error
#>  [1] 0.9055441 0.8243860 0.7554781 0.7018774 0.6660550 0.6490792 0.6499975
#>  [8] 0.6658984 0.6926568 0.7258591
#> 
#> $CI.lower
#>  [1] -0.8401884 -0.7242287 -0.6307299 -0.5645437 -0.5293737 -0.5261738
#>  [7] -0.5519372 -0.5998178 -0.6605904 -0.7244650
#> 
#> $CI.upper
#>  [1] 2.709479 2.507305 2.330690 2.186765 2.081514 2.018170 1.996006 2.010456
#>  [9] 2.054574 2.120850

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

Introduction to RetroJM