Analyzing differences between restricted mean survival time curves using pseudo-values

Hazard ratios are ubiquitously used in time to event analysis to quantify treatment effects. Although hazard ratios are invaluable for hypothesis testing, other measures of association, both relative and absolute, may be used to fully elucidate study results. Restricted mean survival time (RMST) differences between groups have been advocated as useful measures of association. Recent work focused on model-free estimates of the difference in restricted mean survival through follow-up times, instead of focusing on a single time horizon. The resulting curve can be used to quantify the association in time units with a simultaneous confidence band. In this work a model-based estimate of the curve is proposed using pseudo-values allowing for possible covariate adjustment. The method is easily implementable with available software and makes possible to compute a simultaneous confidence region for the curve. The pseudo-values regression using multiple restriction times is in good agreement with the estimates obtained by standard direct regression models fixing a single restriction time. Moreover, the proposed method is flexible enough to reproduce the results of the non-parametric approach when no covariates are considered. Examples where it is important to adjust for baseline covariates will be used to illustrate the different methods together with some simulations. Supplementary Information The online version contains supplementary material available at (10.1186/s12874-022-01559-z).


An illustrative example
For illustration of the pseudo-values method, we use the classical Freireich data [1]. We compare the results obtained using pseudo-values with two other approaches, namely the flexible regression method of Royston and Parmar [2], and the direct regression method of Tian et al. [3], based on weighted estimating equations. Freireich's data deal with a trial in childhood leukemia comparing length of remission in (paired) groups treated with 6-MP or placebo. The data are available in [4] and were used also in [5] for illustrative purposes. In the placebo group there was no censoring.
Note that the non parametric estimate of the restricted mean survival time at τ = 23 is τ 0 S(t)dt = 13.1, as the average of the pseudo-values.   The RMST for the placebo group is estimated equal to 8.38 while γ 1 estimates the difference of RMST(τ =23) between 6-MP and placebo groups.
The analysis can be performed using the approach of Royston and Parmar [2] using flexible parametric survival models. The estimate of the restricted mean survival time at time τ is obtained using a regression model on the log cumulative hazard function. This estimate is then transformed into the survival function and integrated over the interval (0, τ ). Standard errors can be estimated using the bootstrap or the delta method. In the following code the function stpm2, [7], is used to estimate the log cumulative hazard function, then the predict function provides the estimate of the RMST. The methods of Zucker [8] is implemented in the R function restricted.residual.mean from the timereg package [9]. The result is reported below: leuk2$treat<-1*(leuk2$treatment=="6-MP") out<-cox.aalen(Surv(time,status)~prop(treat),data=leuk2,max.timepoint.sim=NULL,resample.ii coxrm <-restricted.residual.mean(out, tau=23, x=rbind(0, 1),iid=1) summary(coxrm) mean se 8.99 1.21 17.90 1.46 The last method illustrated here is that of Tian and colleagues [3] based on inverse probability weighted estimating functions. The implementation is in package survRM2 [10]. When there are no covariates the function uses the non parametric RMST estimator τ 0 S(t)dt, having estimates equal to those reported above.
ipcw.model <-rmst2(leuk2$time, leuk2$status, leuk2$treatment=="6-MP", tau=23) ipcw.model The SAS proc rmstreg can perform the analysis using the method of Tian [3] and also using pseudo-values. For the sake of comparison of the results across different software, we report below the analysis with proc rmstreg with the use of pseudo-values. We use as reference the placebo group for direct comparison with the result obtained in R.   [3] obtained using SAS proc rmstreg, the flexible regression of Royston and Parmar [2] using stpm2 R function, and the method of Zucker [8] from the R function cox-aalen.
Considering pseudo-values, it is possible for example, to do the calculations at τ = 15 and τ = 23. This will result in two pseudo-values for each subject. Stacking the data, it is possible to perform the analysis simultaneously for the two restriction times. Below you can see the code for pseudo-values calculation and the first 20 rows of the stacked data.
bv <-data.frame() for (j in c(15, 23) The regression model with the estimates obtained are reported in the code below.

Quasi-likelihood function for pseudo-observations
The presented application of pseudo-values models for the RMST tries to evaluate the difference between RMST curves through follow-up times. To obtain a smooth function we propose to use regression splines. The spline complexity, i.e. the number of knots, may be chosen on a clinical/application ground. For example if only mild non linearities are expected a spline with few knots may be sufficient. To select the complexity on the basis of a goodness-of-fit statistic we may try to borrow from some criteria specifically designed for generalized estimating equations (GEE).
The pseudo-values are calculated at a set of M selected time points, the pseudo times (τ 1 , . . . , τ M ), as In fact, according to chapter 4 of [12], the Gaussian quasi likelihood is given by − 1 2 (y − µ) 2 . A goodness of fit measure was proposed for GEE by Pan [13], a quasi information criterion (QIC): where N is the naïve variance estimate, considering independent values, while V is the robust variance. This criterion was proposed to select from different correlation structures, while if the goal is selecting covariates, or, such as in our case, the number of knots of the spline, the QIC u can be used, [12], substituting the trace with 2p, where p is the number of parameters in the model, like in the original AIC. This is a pragmatic solution to the problem of model selection while other approaches, more tailored to the specificities of pseudo-values, are emerging in literature, [14], and will hopefully bridge the gap in this framework.