Vignette for survRM2 package:
Comparing two survival curves using the restricted mean survival time

3.1 Unadjusted analysis and its implementation

Let μ_τ(1) and μ_τ(0) denote the RMST for treatment group 1 and 0, respectively. Now, we compare the two survival curves, using the RMST or RMTL. Specifically, we consider the following three measures for the between-group contrast.

1. Difference in RMST μ_τ(1) − μ_τ(0)
2. Ratio of RMST μ_τ(1)/μ_τ(0)
3. Ratio of RMTL {τ − μ_τ(1)}/{τ − μ_τ(0)}

These are estimated by simply replacing μ_τ(1) and μ_τ(0) by their empirical counterparts (i.e.,μ̂_τ(1) and μ̂_τ(0), respectively). For inference of the ratio type metrics, we use the delta method to calculate the standard error. Specifically, we consider log {μ̂_τ(1)} and log {μ̂_τ(0)} and calculate the standard error of log-RMST. We then calculate a confidence interval for log-ratio of RMST, and transform it back to the original ratio scale. Below shows how to use the function, rmst2, to implement these analyses.

time   = D$time
status = D$status
arm    = D$arm

rmst2(time, status, arm, tau=10)

The first argument (time) is the time-to-event vector variable. The second argument (status) is also a vector variable with the same length as time, each of the elements takes either 1 (if event) or 0 (if no event). The third argument (arm) is a vector variable to indicate the assigned treatment of each subject; the elements of this vector take either 1 (if active treatment arm) or 0 (if control arm). The fourth argument (tau) is a scalar value to specify the truncation time point ${\bf \tau}$ for the RMST calculation. Note that τ needs to be smaller than the minimum of the largest observed time in each of the two groups (let us call this the max τ). The program will stop with an error message when such τ is specified. When τ is not specified in rmst2, i.e., when the code looks like

rmst2(time, status, arm)

the max τ is used as the default τ. It is always encouraged to confirm that the size of the risk set is large enough at the specified τ in each group to make sure the stability of the KM estimates.

Below is the output with the pbc example when τ = 10 (years) is specified. The rmst2 function returns RMST and RMTL on each group and the results of the between-group contrast measures listed above.

obj = rmst2(time, status, arm, tau=10)
print(obj)

## 
## The truncation time: tau = 10  was specified. 
## 
## Restricted Mean Survival Time (RMST) by arm 
##               Est.    se lower .95 upper .95
## RMST (arm=1) 7.146 0.283     6.592     7.701
## RMST (arm=0) 7.283 0.295     6.704     7.863
## 
## 
## Restricted Mean Time Lost (RMTL) by arm 
##               Est.    se lower .95 upper .95
## RMTL (arm=1) 2.854 0.283     2.299     3.408
## RMTL (arm=0) 2.717 0.295     2.137     3.296
## 
## 
## Between-group contrast 
##                        Est. lower .95 upper .95     p
## RMST (arm=1)-(arm=0) -0.137    -0.939     0.665 0.738
## RMST (arm=1)/(arm=0)  0.981     0.878     1.096 0.738
## RMTL (arm=1)/(arm=0)  1.050     0.787     1.402 0.738

In the present case, the difference in RMST (the first row of the “Between-group contrast” block in the output) was -0.137 years. The point estimate indicated that patients on the active treatment survive 0.137 years shorter than those on placebo group on average, when following up the patients 10 years. While no statistical significance was observed (p=0.738), the 0.95 confidence interval (-0.665 to 0.939) was relatively tight around 0, suggesting that the difference in RMST would be at most +/- one year.

The package also has a function to generate a plot from the rmst2 object. The following figure is automatically generated by simply passing the resulting rmst2 object to plot() function after running the aforementioned unadjusted analyses.

plot(obj, xlab="Years", ylab="Probability")

3.2 Adjusted analysis and implementation

In most of the randomized clinical trials, an adjusted analysis is usually included in one of the planned analyses. One reason would be that adjusting for important prognostic factors may increase power to detect a between-group difference. Another reason would be we sometimes observe imbalance in distribution of some of baseline prognostic factors even though the randomization guarantees the comparability of the two groups on average. The function, rmst2, in this package implements an ANCOVA type adjusted analysis proposed by Tian et al. [4], in addition to the unadjusted analyses presented in the previous section.

Let Y be the restricted mean survival time, and let Z be the treatment indicator. Also, let X denote a q-dimensional baseline covariate vector. Tian’s method consider the following regression model, g{E(Y ∣ Z, X)} = α + βZ + γ^′X, where g(⋅) is a given smooth and strictly increasing link function, and (α, β, γ^′) is a (q + 2)-dimension unknown parameter vector. Prior to Tian et al. [4], Andersen et al. [5] also studied this regression model and proposed an inference procedure for the unknown model parameter, using a pseudo-value technique to handle censored observations. Program codes for their pseudo-value approach are available on the three major platforms (Stata, R and SAS) with detailed documentation [6, 7]. In contrast to Andersen’s method [5, 6, 7], Tian’s method [4] utilizes an inverse probability censoring weighting technique to handle censored observations. The function, rmst2, in this package implements this method.

As shown below, for implementation of Tian’s adjusted analysis for the RMST, the only the difference is if the user passes covariate data to the function. Below is a sample code to perform the adjusted analyses.

rmst2(time, status, arm, tau=10, covariates=x)

where covariates is the argument for a vector/matrix of the baseline characteristic data, x. For illustration, let us try the following three baseline variables, in the pbc data, as the covariates for adjustment.

x=D[,c(4,6,7)]
head(x)

##        age bili albumin
## 1 58.76523 14.5    2.60
## 2 56.44627  1.1    4.14
## 3 70.07255  1.4    3.48
## 4 54.74059  1.8    2.54
## 5 38.10541  3.4    3.53
## 6 66.25873  0.8    3.98

The rmst2 function fits data to a model for each of the three contrast measures (i.e., difference in RMST, ratio of RMST, and ratio of RMTL). For the difference metric, the link function g(⋅) in the model above is the identity link. For the ratio metrics, the log-link is used. Specifically, with this pbc example, we are now trying to fit data to the following regression models:

1. Difference in RMST E(Y ∣ arm, X) = α + β(arm) + γ₁(age) + γ₂(bili) + γ₃(albumin),
2. Ratio of RMST log {E(Y ∣ arm, X)} = α + β(arm) + γ₁(age) + γ₂(bili) + γ₃(albumin),
3. Ratio of RMTL log {τ − E(Y ∣ arm, X)} = α + β(arm) + γ₁(age) + γ₂(bili) + γ₃(albumin).

Below is the output that rmst2 returns for the adjusted analyses.

rmst2(time, status, arm, tau=10, covariates=x)

## 
## The truncation time: tau = 10  was specified. 
## 
## Summary of between-group contrast (adjusted for the covariates) 
##                        Est. lower .95 upper .95     p
## RMST (arm=1)-(arm=0) -0.210    -0.883     0.463 0.540
## RMST (arm=1)/(arm=0)  0.968     0.877     1.068 0.514
## RMTL (arm=1)/(arm=0)  1.035     0.806     1.329 0.786
## 
## 
## Model summary (difference of RMST) 
##             coef se(coef)      z     p lower .95 upper .95
## intercept  2.743    2.134  1.285 0.199    -1.440     6.927
## arm       -0.210    0.343 -0.613 0.540    -0.883     0.463
## age       -0.069    0.018 -3.900 0.000    -0.103    -0.034
## bili      -0.325    0.039 -8.386 0.000    -0.401    -0.249
## albumin    2.550    0.472  5.401 0.000     1.624     3.475
## 
## 
## Model summary (ratio of RMST) 
##             coef se(coef)      z     p exp(coef) lower .95 upper .95
## intercept  1.369    0.356  3.842 0.000     3.930     1.955     7.899
## arm       -0.033    0.050 -0.652 0.514     0.968     0.877     1.068
## age       -0.009    0.003 -3.410 0.001     0.991     0.985     0.996
## bili      -0.087    0.013 -6.523 0.000     0.917     0.893     0.941
## albumin    0.360    0.080  4.491 0.000     1.434     1.225     1.678
## 
## 
## Model summary (ratio of time-lost) 
##             coef se(coef)      z     p exp(coef) lower .95 upper .95
## intercept  1.992    0.695  2.865 0.004     7.332     1.876    28.655
## arm        0.035    0.127  0.272 0.786     1.035     0.806     1.329
## age        0.025    0.007  3.810 0.000     1.026     1.012     1.039
## bili       0.063    0.008  8.334 0.000     1.065     1.049     1.080
## albumin   -0.750    0.149 -5.033 0.000     0.472     0.353     0.633

The first block of the output is a summary of the adjusted treatment effect. Subsequently, a summary for each of the three models are provided.