Assessing treatment effects with adjusted restricted mean time lost in observational competing risks data

Abstract

Background

According to long-term follow-up data of malignant tumor patients, assessing treatment effects requires careful consideration of competing risks. The commonly used cause-specific hazard ratio (CHR) and sub-distribution hazard ratio (SHR) are relative indicators and may present challenges in terms of proportional hazards assumption and clinical interpretation. Recently, the restricted mean time lost (RMTL) has been recommended as a supplementary measure for better clinical interpretation. Moreover, for observational study data in epidemiological and clinical settings, due to the influence of confounding factors, covariate adjustment is crucial for determining the causal effect of treatment.

Methods

We construct an RMTL estimator after adjusting for covariates based on the inverse probability weighting method, and derive the variance to construct interval estimates based on the large sample properties. We use simulation studies to study the statistical performance of this estimator in various scenarios. In addition, we further consider the changes in treatment effects over time, constructing a dynamic RMTL difference curve and corresponding confidence bands for the curve.

Results

The simulation results demonstrate that the adjusted RMTL estimator exhibits smaller biases compared with unadjusted RMTL and provides robust interval estimates in all scenarios. This method was applied to a real-world cervical cancer patient data, revealing improvements in the prognosis of patients with small cell carcinoma of the cervix. The results showed that the protective effect of surgery was significant only in the first 20 months, but the long-term effect was not obvious. Radiotherapy significantly improved patient outcomes during the follow-up period from 17 to 57 months, while radiotherapy combined with chemotherapy significantly improved patient outcomes throughout the entire period.

Conclusions

We propose the approach that is easy to interpret and implement for assessing treatment effects in observational competing risk data.

Introduction

According to long-term follow-up data for patients with malignant tumors, individuals may face multiple possible outcomes. Therefore, assessing treatment effects requires careful consideration of competing risks [1, 2]. For instance, the clinical outcomes of cervical cancer patients might include survival, loss to follow-up, death from cervical cancer or death from another cause (such as heart disease, other cancers, circulatory system disease and so on) [3,4,5,6,7]. In most studies on the survival prognosis of cervical cancer patients, researchers tend to consider either a single endpoint (i.e., cervical cancer death, with other unobserved outcomes considered censored) or a combined endpoint (i.e., all-cause death, where several different causes of death are considered). However, when a study is designed to evaluate the effectiveness of a treatment for cervical cancer, the event of interest will be death attributed to cervical cancer. In such cases, if patients who die from other causes are simply treated as censored, this may lead to an overestimation of the occurrence rate of the events of interest, resulting in competing risks bias [8,9,10].

In the competing risks framework, the cause-specific Cox regression model is commonly employed to assess etiological treatment effects using the cause-specific hazard ratio (CHR), while the Fine-Gray model is used to explore the prognostic impact of treatment outcomes through the sub-distribution hazard ratio (SHR). The CHR measures the effect of a treatment group on the rate of dying from a cause of interest, while the SHR measures the effect of a treatment group on the risk of dying from a cause of interest [11]. Depending on the purpose of the study, both indicators have been widely used in competing risk data. However, in practical applications, both the CHR and SHR are hazard-based measures potentially subject to limitations related to practical interpretation and the proportional hazards assumption. Due to the fact that CHR and SHR are both relative indicators, it is challenging to provide clinicians and patients with tangible clinical explanations. Especially given that for patients diagnosed with this disease, the answer to "How much longer can life be extended after receiving treatment?" might be the most pressing and relevant concern. Blagoev [12] discussed considerations for assessing treatment efficacy using hazard ratios in clinical trial data. They proposed, “While a hazard ratio has some value, for the clinician caring for a patient and, more importantly, the patient, it does not convey benefit in terms that are meaningful—how much longer will the patient live or live without experiencing disease progression” [12].

To reduce the constraints of proportional hazards assumption and better explain treatment effects, in recent years, several researchers in the fields of medicine, epidemiology, and statistics have recommended incorporating restricted mean time lost (RMTL) to complement hazard-based indicators in the competitive risk setting [13,14,15,16,17,18,19]. The RMTL represent the average time lost by patients due to the event of interest from the initiation of follow-up to a particular time point [20, 21]. In practical applications, the difference of the RMTL (RMTLd) between patients receiving a certain treatment and those not receiving the treatment is usually calculated to measure the treatment effect. For example, if the time point is prespecified as 5 years and the RMTLd for cervical cancer patients in the treated group compared to the untreated group is -1 year, this can be interpreted as follows: treatment decreases life expectancy due to cervical cancer death by approximately one year.

On the other hand, in practical applications, for certain rare diseases with low incidence rates, it may be challenging to recruit participants for controlled clinical trials. Therefore, retrospective studies can be employed to explore treatment benefits. For example, small cell carcinoma of the cervix (SCCC) is a rare subtype of cervical cancer that accounts for only 0.5–1.0% of all cervical cancers [22,23,24], with a 5-year overall survival rate of less than 30% [24,25,26]. Due to the low incidence of SCCC, recommendations in clinical guidelines are unclear and researchers face challenges in conducting prospective randomized controlled studies. However, in observational studies, the distribution of baseline covariates in the treatment and non-treatment groups is likely to be unbalanced. If these covariates also predict potential survival benefits, the effect of treatment on survival will be confounded, and direct intergroup comparison methods may not result in correct inferences [27, 28]. Therefore, it is crucial to adjust covariates in observational studies. The inverse probability weighting (IPW) method is based on the propensity score and is a flexible covariate adjusting method in application, such that there have been studies that have considered adjusted RMTL using IPW [17, 18, 29]. Therefore, this study focused on obtaining accurate and robust causal effect estimates of RMTLd after adjusting for confounding factors using IPW.

Furthermore, conducting an analysis using RMTLd requires specifying the restricted time point in advance. Typically, clinicians choose time points of interest, such as 3-year RMTLd or 5-year RMTLd, or they empirically select the minimum among the maximum observation times in the two groups [16]. However, in practical analytical situations, especially in observational studies, specifying a single time point may result in partial information loss. Therefore, we considered the dynamic changes in the RMTLd at various time points, aiming to provide additional information about clinical efficacy.

In this study, for competing risks data, we used the IPW method to develop point and interval estimates of weighted RMTLd, allowing for covariate adjustments. A Monte Carlo simulation study was used to explore the statistical performance of the estimator under various situations. Furthermore, considering that treatment effects on patient prognosis between groups may vary over time, we extended our approach to develop dynamic adjusted RMTLd curves and simultaneous confidence bands. The implementation of the proposed method with a simulated example is available on Git-Hub [https://github.com/chenzgz/aRMTL]. Finally, we applied this method to patients with SCCC from the Surveillance, Epidemiology, and End Results (SEER) database to evaluate the treatment effects of different therapeutic approaches on cervical cancer.

Method

For generality, we considered a model with only one event of interest (e.g., $j=1$, death from studied cancer) and one competing event (e.g., $j=2$, death from other causes). Let $T$ be the event time, and the corresponding event type is $j\left(j=1,2\right)$. $C$ is the independent censoring time; then, the observation time is $X=min\left\{T,C\right\}$ , and the indicator of the occurrence of the event is $\delta=I\left(T\leq C\right)j\;(\delta=0$ is the censoring status, and $\delta=1,2$  are events of interest and competing events, respectively). Assuming a comparison between two groups (e.g., whether the patients received treatment), where $K=k\left(k=0,1\right)\;(k=1$indicates the treatment group and $k=0$ indicates the non-treatment group) represents the grouping of patients, we define $\mathrm Z$ as a baseline covariate matrix of dimension $n\times q$ . For patients i = 1…, n, the observed data vector is $\left(X_i,\delta_i,K_i,Z_i^T\right)^T$. Additionally, for patient i in the k-th group, we define $N_{ijk}\left(t\right)=I\left(X_i\leq t,\delta_i=j,K_i=k\right)$, which indicates whether the individual has occurred event j by time t, and $Y_{ik}\left(t\right)=I\left(X_i\geq t,K_i=k\right)$, which indicates whether any event has not occurred to the individual by time t.

In the presence of competing risks, the cause-specific cumulative incidence function (CIF) for the k-th group's event of interest is defined as $F_{1k}\left(t\right)=P\left(T\leq t,j=1,K=k\right)=\int_0^texp\left(-\wedge_{1k}\left(s\right)-\wedge_{2k}\left(s\right)\right)d\wedge_{1k}\left(s\right)$, where $\wedge_{jk}\left(t\right)$ is the cause-specific cumulative hazard function (Nelson‒Aalen estimate) for event j and group k. To compare treatment effects between groups, the measures of interest in this study are the RMTLd for the event of interest between the two groups: $\triangle\left(t\right)=\mu_1\left(t\right)-\mu_0\left(t\right)$, where $\mu_k\left(t\right)=\int_0^tF_{1k}\left(u\right)du$. In practice, it is common to choose a restricted time point, denoted as $t=\tau$. $\mu_k\left(\tau\right)$represents the RMTL for the event of interest in the k-th group up to time $\tau$, calculated by integrating the area under the CIF curve for the event of interest from 0 to$\tau$. The RMTLd is interpreted as the discrepancy in the average time lost due to the event of interest between the two groups up to time $\tau$.

Adjusted difference of the RMTL estimator using the inverse probability weighted method

In order to adjust for the covariates using IPW method, it is necessary to estimate propensity scores to construct weights. We employed a logistic regression model to estimate propensity scores, i.e., $\log\left\{p\left(\beta\right)/\left[1-p\left(\beta\right)\right]\right\}=\beta\left(1,\mathrm Z^T\right)^T$, where $p\left(\beta\right)$  represents the probability of a patient being assigned to the treatment group. The maximum likelihood estimation yields the propensity score estimate $p_i\left(\widehat{\beta}\right)$ for the i-th patient, and the inverse probability weight estimate is $w_i\left(\widehat{\beta}\right)=I\left(K_i=1\right)/{p_i\left(\widehat{\beta}\right)}+I\left(K_i=0\right)/\left[1-p_i\left(\widehat{\beta}\right)\right]$.

In this study, we construct a weighted RMTLd based on the Nelson–Aalen estimate of IPW under competing risks conditions. Following the concept of IPW, assigning weights to each individual, the IPW weighted Nelson-Aalen estimator under competing risks is given by

$$\widehat{\wedge}_{jk}^w\left(t,\widehat{\beta}\right)=\int_0^t\frac{dN_{jk}^w\left(s,\widehat{\beta}\right)}{Y_k^w\left(s,\widehat{\beta}\right)},$$

where $N_{jk}^w\left(t,\widehat{\beta}\right)=\sum_{i=1}^nw_i\left(\widehat{\beta}\right)N_{ijk}\left(t\right)$, is the weighted number of individuals who have experienced the j-event by time t, and $Y_k^w\left(t,\widehat{\beta}\right)=\sum_{i=1}^nw_i\left(\widehat{\beta}\right)Y_{ik}\left(t\right)$, is the weighted number of individuals who have not experienced any event by time t.

Then, the IPW weighted RMTLd estimate is given by $\widehat{\triangle}^w\left(t,\widehat{\beta}\right)=\widehat{\mu}_1^w\left(t,\widehat{\beta}\right)-\widehat{\mu}_0^w\left(t,\widehat{\beta}\right)$, where $\widehat{\mu}_k^w\left(t,\widehat{\beta}\right)=\int_0^t\widehat{F}_{1k}^w\left(s,\widehat{\beta}\right)ds$ and$\widehat{F}_{1k}^w\left(t,\widehat{\beta}\right)=\int_0^texp\left\{-\widehat{\wedge}_{1k}^w\left(s,\widehat{\beta}\right)-\widehat{\wedge}_{2k}^w\left(s,\widehat{\beta}\right)\right\}d\widehat{\wedge}_{1k}^w\left(s,\widehat{\beta}\right)$.

For a restricted time point $\tau$, $\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)$ is interpreted as the difference, after controlling for covariates, in the average time lost due to cervical cancer death between the entire study population receiving treatment and those not receiving treatment within $\tau$ years. To ensure the validity of the estimates, in practice, the choice of $\tau$ should neither be too early nor too late, meaning that each group should have at least one event before $\tau$ (i.e., $Pr\left(T<\tau\right)>0$) and at least one observation after $\tau$ (i.e., $Pr\left(X>\tau\right)>0$). This study refers to the RMTLd constructed based on the IPW as the adjusted RMTLd (aRMTLd).

The variance estimate of $\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)$ can be approximately derived from its large-sample properties (for the detailed process, see Web Appendix 1), which is denoted as $\widehat{\sigma}^w\left(\tau\right)^2=n^{-1}\sum_{i=1}^n\left(\widehat\xi_{i1}^w\left(\tau,\widehat{\beta}\right)-\widehat\xi_{i0}^w\left(\tau,\widehat{\beta}\right)\right)^2$, where $\widehat\xi_{ik}^w\left(t,\widehat{\beta}\right)$ is defined in formula (7) in Web Appendix 1. Then, the $\left(1-\alpha\right)$ two-sided confidence interval for $\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)$ is $\left(\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)-z_{\left(1-\alpha/2\right)}n^{-1/2}\widehat{\sigma}^w\left(\tau\right),\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)+z_{\left(1-\alpha/2\right)}n^{-1/2}\widehat{\sigma}^w\left(\tau\right)\right)$, where $z_{\left(1-\alpha/2\right)}$ is the $100\left(1-\alpha/2\right)$ percentile of the standard normal distribution and $\widehat{\sigma}^w\left(\tau\right)$ is the estimated standard deviation of $\widehat{\triangle}^w\left(\tau,\widehat{\beta}\right)$.

Dynamic curve of aRMTLd and simultaneous confidence band estimates

In some practical applications, the effect of treatment may often change over time, and a single time point alone still cannot provide a sufficient summary of the treatment effect. At this point, it is necessary to assess the variation of aRMTLd over a time interval and using simultaneous confidence bands to evaluate its variability. Due to the complexity of $\xi_{ik}^w\left(t,\beta\right)$, our study employs a perturbation-resampling method to construct confidence bands for the curve.

Assuming that the study evaluates the effect of treatment on patients from time points $t_L$ to $t_U$, the simultaneous equal-precision confidence band [30] for $\widehat{\triangle}^w\left(t,\widehat{\beta}\right)$ is estimated as $\left(\widehat{\triangle}^w\left(t,\widehat{\beta}\right)-{\widehat q}_\alpha n^{-1/2} \widehat{\sigma}^w\left(t\right),\widehat{\triangle}^w\left(t,\widehat{\beta}\right)+{\widehat q}_\alpha n^{-1/2} \widehat{\sigma}^w\left(t\right)\right)$, and the estimated value of the $100\left(1-\alpha\right)$ empirical percentile ${\widehat q}_\alpha$ for the simultaneous confidence bands is obtained using the following formula:

$$Pr\left\{{\begin{array}{c}sup\\t\epsilon[t_L,t_U]\end{array}\left|\frac{n^{-1/2}{\sum\nolimits_{i=1}^n}\left(\widehat\xi_{i1}^w\left(t,\widehat{\beta}\right)-\widehat\xi_{i0}^w\left(t,\widehat{\beta}\right)\right)R_i}{\widehat{\sigma}^w\left(t\right)}\right|}>{\widehat q}_\alpha\right\}=\alpha$$

where $R_i$ (i = 1,..,n) is a randomly sampled standard normal distribution variable. Similarly, it should be noted that, for the selection of the time interval $\left[t_L,t_U\right]$, in practice, to make the most of the follow-up information possible, one can empirically choose the maximum time interval. Specifically, set $t_L$ to be the minimum time satisfying $Pr\left(T<t_L\right)>0$, and set $t_U$ to be the maximum time satisfying $Pr\left(X>t_U\right)>0$. If one simultaneously considers a time interval of clinical interest, the final choice may be the intersection of these two time intervals [31].

Simulation studies

Simulation design

We designed a Monte Carlo simulation study to investigate the performance of aRMTLd at restricted time points. Our simulation study was based on the simulation design framework by Conner [17]. We adjusted some parameters to consider more scenarios. The simulation is conducted in two groups of comparisons, setting up an event of interest ($j=1$) and a competing event ($j=2$). The comparison methods included RMTLd without adjustment for covariates and aRMTLd with IPW after covariate adjustment. Time points $\tau=182.5,360,730$ days were selected for analysis. The simulation design is carried out under the condition that the model assumptions are correct. Four covariates are set to affect the allocation of treatment groups and are independent prognostic factors at the same time, that is, they are all confounding factors and need to be adjusted. The reported performance measures include the mean bias, empirical root mean squared error (RMSE), empirical coverage, and relative standard error (RSE).

Data generation

First, we generate four independent covariates $\mathrm Z=\left({\mathrm Z}_1,{\mathrm Z}_2,{\mathrm Z}_3,{\mathrm Z}_4\right)^T$ from the standard normal distribution. The grouping variable $K_i$ is then generated from the Bernoulli distribution with probability $p=1/\left(1+exp\left(-\boldsymbol\alpha\boldsymbol Z\right)\right)$, where $\boldsymbol\alpha\mathit=\left(\ln\left(0.6\right),\ln\left(0.5\right),\ln\left(0.7\right),\ln\left(0.8\right)\right)$. Because of $E\left[\boldsymbol\alpha\boldsymbol Z\right]=0$, the individuals assigned to the two groups are basically equal, which is a balanced sample size.

Let $\widetilde{\mathbf Z}=\left(K,\mathbf Z^T\right)^T$, For the event of interest ($j=1$), assuming that the event time obeys Gompertz distribution, the baseline risk function is $\lambda_{10}\left(t\right)=\gamma exp\left(\rho t\right)$, where $\gamma$ and $\rho$ represent the corresponding parameters of the Gompertz distribution. The sub-distribution risk function based on covariates is $\lambda_1\left(t;\widetilde{\mathbf Z}\right)=\lambda_{10}\left(t\right)exp\left\{{\boldsymbol\beta}_1\overset{\boldsymbol\sim}{\mathbf Z}+{\boldsymbol\beta}_2\overset{\boldsymbol\sim}{\mathbf Z}t\right\}$, where ${\boldsymbol\beta}_1=\left\{\beta_{11}\mathit,\ln\left(0.8\right),\ln\left(0.7\right),\ln\left(0.6\right),\ln\left(0.5\right)\right\}$, ${\boldsymbol\beta}_2=\left\{\beta_{21},0,0,0,0\right\}$, where $\beta_{11}$ and $\beta_{21}$ represent the fixed effects and time-varying effects of the grouping variables respectively. The conditional cumulative incidence function of events of interest based on covariates is $F_1\left(t;\widetilde{\mathbf Z}\right)=P\left(T\leq t,j=1\left|\widetilde{\mathbf Z}\right.\right)=1-exp\left[-\int_0^t\gamma exp\left\{{\boldsymbol\beta}_1\overset{\boldsymbol\sim}{\mathbf Z}+{\boldsymbol\beta}_2\overset{\boldsymbol\sim}{\mathbf Z}s+\rho s\right\}ds\right]$. Then when ${\boldsymbol\beta}_2\overset{\boldsymbol\sim}{\mathbf Z}+\rho\mathit<\mathit0$, the marginal probability of the event of interest is $P\left(j=1\left|\widetilde{\mathbf Z}\right.\right)=1-exp\left\{\frac{\gamma exp\left\{{\boldsymbol\beta}_1\widetilde{\mathbf Z}\right\}}{{\boldsymbol\beta}_2\widetilde{\mathbf Z}+\rho}\right\}$. Therefore, the cumulative incidence function of competing events is $F_2\left(t;\widetilde{\mathbf Z}\right)=P\left(T\leq t,j=2\left|\widetilde{\mathbf Z}\right.\right)=exp\left\{\frac{\gamma exp\left\{{\boldsymbol\beta}_1\widetilde{\mathbf Z}\right\}}{{\boldsymbol\beta}_2\widetilde{\mathbf Z}+\rho}\right\}\left[1-exp\left\{-t\;exp({\boldsymbol\beta}_3\widetilde{\mathbf Z})\right\}\right]$, where ${\boldsymbol\beta}_3=\left\{\ln\left(0.67\right),\ln\left(0.5\right),\ln\left(0.6\right),\ln\left(0.7\right),\ln\left(0.8\right)\right\}$.

For the generation of event type and event time, first generate event type j from Bernoulli distribution according to probability $P\left(j=1\left|\widetilde{\mathbf Z}\right.\right)$. Then simulate a random variable from uniform distribution $U\left(0,1\right)$ to determine the inverse of the conditional cumulative incidence function at time t to generate event time. For $j=1,2$, the corresponding conditional cumulative incidence function is $P\left(T\leq t\;\left|\;j=1,\widetilde{\mathbf Z}\right.\right)=F_1\left(t\right)/P\left(j=1\left|\widetilde{\mathbf Z}\right.\right)$ and $P\left(T\leq t\;\left|\;j=2,\widetilde{\mathbf Z}\right.\right)=F_2\left(t\right)/P\left(j=2\;\left|\widetilde{\mathbf Z}\right.\right)$ respectively.

At the same time, an independent right censoring situation is set up. By generating a censoring time C that obeys a uniform distribution $U\left(0,c\right)$, different parameters c are selected to produce censoring rates of 15% and 30%. The final observed time is the minimum value of the event time T and the censoring time C, $\delta=I\left(T\leq C\right)j$ indicates whether an observation is censored and serves as an indicator for the type of event that occurred.

Scenarios and true values

By setting different parameters $\left(\gamma,\rho,\beta_{11},\beta_{21}\right)$, we considered three scenarios, namely Scenario A $\left(1.5,-1,5,\ln\left(0.7\right),0\right)$ when the risk rate of the sub-distribution of the grouping effect is proportional, Scenario B $\left(1.5-1,5,\ln\left(0.6\right),\ln\left(1.8\right)\right)$ of the attenuation effect and Scenario C $\left(1.5,-1.2,\ln\left(1.1\right),\ln\left(0.4\right)\right)$ of the delayed effect when the risk rate of the sub-distribution of the grouping effect is disproportionate. A total of 12 situations were considered with total sample sizes of 500 and 1000 and censoring rates of 15% and 30%.

We determine the true marginal difference in RMTL by generating a sample of size n = 1,000,000. We first generate data for the treatment group ($k=1$), generating event types and corresponding event times based on the simulation settings. The process is then repeated for the treatment group ($k=0$), producing two sets of event results (event type and time) for each person. Selecting time points $\tau=182.5,360,730$ days (corresponding to 0.5, 1, and 2 years respectively), the true marginal difference of the events of interest is the RMTLd of $k=1$ and $k=0$. The real cumulative incidence function is shown in Fig. 1. The corresponding calculated true value of marginal RMTLd is shown in Supplemental Table 1.

Fig. 1

Real CIF under large samples under three scenarios in simulation study

Assessment of statistical performance

In this study, the number of simulations is set to $n_{sim}=5000$. We estimated the aRMTLd in the i-th repetition as ${\widehat\triangle}_i$($i=1,2,...,n_{sim}$) and the true value of the estimand as Δ, which was obtained from Scenarios and true values section. We examined four performance measures in estimating the true marginal RMTLd: the mean bias, $\frac1{n_{sim}}\sum_{i=1}^{n_{sim}}\left({\widehat\triangle}_i-\triangle\right)$; the RMSE, $\sqrt{\frac1{n_{sim}}\sum_{i=1}^{n_{sim}}\left({\widehat\triangle}_i-\triangle\right)^2}$; the empirical coverage rate, estimated as the proportion of 95% confidence intervals that covered $\widehat\triangle$; and the RSE is defined as the ratio of the average model standard error and Monte Carlo empirical standard error, $\frac{\sqrt{{\displaystyle\frac1{n_{sim}}}{\sum\limits_{i=1}^{n_{sim}}}\widehat{Var}\left({\widehat\triangle}_i\right)}}{\sqrt{{\displaystyle\frac1{n_{sim}-1}}{\sum\limits_{i=1}^{n_{sim}}}\left({\widehat\triangle}_i-\triangle\right)^2}}$.

Simulation results

The simulation results demonstrate that the unadjusted RMTLd exhibits substantial average bias across different scenarios, with the deviations increasing as the time points increase. Specifically, within 180.5 days, 365 days, and 730 days, the unadjusted RMTLd deviates from the true values by approximately 18 to 20 days, 46 to 50 days, and 104 to 108 days, respectively. After adjusting for covariates, the aRMTLd consistently demonstrated a small average bias across various scenarios and time points, with differences from the true values consistently occurring within 0.42 days. Similarly, the RMSE of aRMTLd was consistently lower than of unadjusted RMTLd and decreased as the sample size increased. The coverage of unadjusted RMTLd is generally less than 10%, remaining below 5% as time progresses, which due to the substantial bias. The coverage of aRMTLd approaches the nominal level of 95%, mostly remaining within a reasonable range (94.4% to 95.6%) based on 5000 simulations. Additionally, the estimated RSE for aRMTLd is close to 1, indicating that this method provides standard deviation estimates that closely resemble true variability. In summary, the adjusted RMTL estimator exhibits smaller biases compared with unadjusted RMTL in all scenarios and provides robust interval estimates. Table 1 present the simulation performance comparisons of the RMTLd and aRMTLd estimators at different time points.

Table 1 Results of the simulation study comparing the RMTLd and aRMTLd estimators

Application of the adjusted RMTLd in patients with small cell carcinoma of the cervix

The data for this study were sourced from the SEER database. We extracted information on patients diagnosed with SCCC between 2000 and 2020. The detailed data collected from the SEER database and variable classification are provided in Web Appendix 2. In this study, the event of interest for patients with SCCC was death due to cervical cancer, while death from other causes was considered a competing event.

Previous studies have indicated that surgery, radiotherapy, and chemotherapy each have certain therapeutic effects on SCCC. In this study, we investigated the impact of surgery, radiotherapy, and combined radiotherapy with chemotherapy on the prognosis of patients with SCCC. First, we examined the distribution of covariates across different treatment groups (Supplemental Tables 3–5). To estimate treatment effects, we included all variables except for the treatment group when constructing a logistic regression to estimate weights, allowing us to construct the points and intervals of the aRMTLd estimates. The unadjusted RMTLd was used for comparison. To further explore the time-varying effects of treatment, we employed aRMTLd curves and simultaneous confidence bands to depict the dynamic changes. The analysis considered fixed time points at 36 months (3 years), 60 months (5 years), and 120 months (10 years). The time interval for assessing the dynamic aRMTLd curve was chosen as [3, 120] months.

Additionally, we present the results from Fine-Gray models, employing both univariate and multivariate regression models to derive the SHR and adjusted SHR (aSHR). Notably, the SHR can be considered the weighted average of various time points before a specified time point [32]. The study also conducted tests on the proportional hazards (PH) assumption for the sub-distribution hazards of the model [33].

Surgical groups

There was an imbalance in the distributions of all the other covariates between the groups of patients who did and did not undergo surgical treatment (Supplemental Table 3). The results of Fine-Gray model showed that SHR < 1, suggesting surgery as a protective factor with statistical significance, while the aSHR value increased compared with SHR, and its confidence interval included 1 (Table 2). Nevertheless, at all three time points, surgery did not satisfy the PH assumption according to neither univariate nor multivariate analyses. This finding indicates that the SHR and aSHR may vary over time, and a single regression coefficient may not accurately indicate treatment efficacy.

Table 2 The SHR and RMTL were used to analyze the prognostic value of different treatment methods in patients with SCCC

Compared with the unadjusted CIF curves, the difference between the two CIF curves decreased after we adjusted for covariates (Fig. 2-A). The absolute values of aRMTLd are smaller than those of the unadjusted RMTLd. Taking 36 months as an example, the RMTLd is -7.56 m (95% CI: -9.65 ~ -5.48 m), while the aRMTLd is -1.95 m (95% CI: -4.66 ~ 0.77 m) (Table 2). The interpretation is that within 36 months, when other covariates are controlled, compared to patients with SCCC who did not receive surgery, patients with SCCC who did receive surgery experienced an average reduction of 1.95 months in the time lost due to death from cervical cancer. At all three time points, consistent results were obtained using the Fine-Gray model and RMTL estimates. The univariate analysis indicated significant improvements in patient prognosis after surgery, whereas after covariate adjustment, the differences were not statistically significant. This finding suggests that the unadjusted results might be influenced by confounding factors, emphasizing the importance of evaluating treatment effects accurately using results adjusted for covariates.

Fig. 2

Unadjusted and adjusted CIF curves and dynamic RMTL difference curves of three treatment groups. A CIF plots for the surgery group, both unadjusted and adjusted; B dynamic RMTLd and aRMTLd curves for the surgery group; C CIF plots for the radiotherapy group, both unadjusted and adjusted; D dynamic RMTLd and aRMTLd curves for the radiotherapy group; E CIF plots for the combined radiotherapy with chemotherapy group, both unadjusted and adjusted; F dynamic RMTLd and aRMTLd curves for the combined radiotherapy with chemotherapy group. Note: the time range for dynamic aRMTLd curve plots is selected as [3, 120] months; in Fig. 2-B, D and F, the dashed lines represent pointwise confidence intervals, and the red shading represents the simultaneous confidence band for the aRMTLd curve

The results of the dynamic aRMTLd curve analysis showed that the absolute values of the unadjusted RMTLd were relatively large, and the confidence intervals did not include 0. In contrast, the aRMTLd absolute values decreased, and the confidence intervals included 0 from 24 to 120 months, with the simultaneous confidence band including 0 from 20 to 120 months (Fig. 2-B). The findings suggest that the statistically significant impact of surgery on patient outcomes is limited to the initial 20 months of follow-up, with no significant long-term improvement in prognosis.

Radiotherapy groups

There was an imbalance in the distribution of covariates, except for tumor size and lymph node involvement, between the groups of patients who did and did not receive radiotherapy (Supplemental Table 4). The results of the Fine‒Gray model showed that SHR was < 1 while after adjusting for covariates, the aSHR increased, and its confidence interval included 1 (Table 2). Similarly, at all three time points, radiotherapy did not satisfy the PH assumption according to neither univariate nor multivariate analyses.

Compared with the unadjusted CIF curves, the difference between the two CIF curves decreased after adjusting for covariates (Fig. 2-C). The absolute values of aRMTLd were smaller than those of the unadjusted RMTLd. Taking 60 months as an example, the RMTLd was -8.07 m (95% CI: -12.19 ~ -3.95 m), while the aRMTLd was -5.65 m (95% CI: -10.26 ~ -1.05 m) (Table 2). The interpretation is that within 60 months, when other covariates are controlled, patients with SCCC who received radiotherapy experienced an average reduction of 5.65 months in the time lost due to death from cervical cancer compared to those who did not. In other words, within 60 months, radiotherapy could increase the life expectancy free of cervical cancer death by 5.65 months.

The results of the dynamic aRMTLd curve analysis showed that the absolute values of the unadjusted RMTLd were relatively large, and the confidence intervals did not include 0. In contrast, the aRMTLd absolute values decreased, and the confidence intervals excluded 0 from 3 to 90 months, with the simultaneous confidence band excluding 0 from 17 to 57 months (Fig. 2-D). In the interpretation of the simultaneous confidence band, taking 36 months as an example, the aRMTLd was -3.95 months, with a simultaneous confidence band of (-7.35 m, -0.54 m). This finding indicates that, after adjustment for multiple comparisons over the entire time interval [3, 120] using the simultaneous confidence band, radiotherapy could still significantly increase life expectancy free of cervical cancer death by 3.95 months within 36 months.

Treatment groups: combined radiotherapy with chemotherapy

Regarding the group of patients who received both radiotherapy and chemotherapy treatment or not, there was an imbalance in the distribution of covariates, except for tumor size, lymph node involvement, and whether surgery was performed, with statistically significant differences (Supplemental Table 5). The aSHR < 1 was considered to indicate a statistically significant difference, suggesting combined radiotherapy with chemotherapy as a protective factor with statistical significance after adjusting for covariates (Table 2). Similarly, at all three time points, combined radiotherapy with chemotherapy did not satisfy the PH assumption according to neither univariate nor multivariate analyses.

Compared with the unadjusted CIF curves, the difference between the two CIF curves decreased after adjusting for covariates (Fig. 2-E).The absolute values of the aRMTLd were smaller than those of the unadjusted RMTLd, with statistically significant differences observed before and after adjustment. Taking 120 months as an example, the aRMTLd was -10.90 m (95% CI: -19.23 ~ -2.58 m) (Table 2). This finding implied that within 120 months, when other covariates are controlled, patients with SCCC in the combined radiotherapy with chemotherapy group could achieve an additional 10.90 months of life expectancy.

The confidence band of the aRMTLd curve does not include 0 during the entire period (Fig. 2-F). This indicates that the improvement in prognosis for patients receiving combined radiotherapy with chemotherapy was statistically significant within the [3, 120] months duration.

Discussion

In summary, this study considered assessing treatment effects based on time scales in observational competing risk data. Depending on the purpose of the study, when the end point of interest is death from a specific cause, competing risks should be carefully considered to reduce bias (such as an overestimation of the occurrence rate of the events of interest). On the other hand, when observational study data are used to assess treatment effects, adjustment for covariates is necessary. The aRMTLd constructed in this study based on IPW, by directly weighting each individual under given covariate conditions, can be applied to observational studies with competing risk data, providing a more accurate assessment of treatment effects after adjusting for covariates. The simulation study results indicate that in the presence of confounding factors, estimating the RMTLd without covariate adjustment leads to substantial bias. In contrast, the aRMTLd demonstrates little bias, and the interval estimate is robust. In the prognostic study of patients with SCCC, the analysis of the three treatments revealed that unadjusted RMTLd indicated a significant improvement in patient prognosis for all treatments. However, after adjustment, the effects of aRMTLd were reduced, especially after surgery, where long-term efficacy was not significant. This could be attributed to the fact that patients receiving treatment may have better baseline health conditions than those not receiving treatment (e.g., younger age, smaller tumor size, or earlier stage) [34, 35]. Without controlling for the impact of covariates, the treatment effect may be exaggerated, leading to potentially erroneous conclusions. Overall, the aRMTLd can adjust for covariates and complement risk-based measures, providing a "time" perspective for estimating treatment effects.

To date, some studies have considered RMTL adjusted for covariates and have developed corresponding models. Conner and Trinquart [17] applied IPW combined with RMTL to construct estimators and variance. Notably, they used a weighted Kaplan–Meier estimator to estimate the overall survival function in the weighted RMTL, while our study used a weighted Nelson-Aalen estimator to estimate the overall survival. These two estimators are asymptotically equivalent, but the Nelson–Aalen estimator has better finite-sample properties [36]. Furthermore, previous works have treated the IPW estimator weights as fixed and derived the variance using Greenwood's formula. The proposed variance of the RMTL does not account for variability in the estimated weights [17]. In contrast, our study considers the variability in weights estimated by logistic regression when deriving the variance. Therefore, simulation results show that the coverage estimated by our method are close to 0.95. The comparative results of these two methods under various simulation scenarios are provided in Supplemental Table 2.

Globally, cervical cancer is the fourth most common cancer among women, with SCCC being a relatively rare subtype. In the National Comprehensive Cancer Network (NCCN) guidelines [37], surgery, chemotherapy and radiotherapy are advocated as the primary treatment modalities for SCCC. Previous studies have indicated that primary radiotherapy combined with aggressive chemotherapy yields better survival outcomes than surgery [35]. However, the efficacy of surgery for treating SCCC remains controversial. Several retrospective studies have reported benefits from surgery [25, 34, 38]. In contrast, our study showed that the difference in surgical efficacy was statistically significant only within the first 20 months of follow-up. In the long term, surgery does not significantly impact the prognosis of patients with SCCC, possibly because 2–3 years after surgical treatment is the period of high recurrence of cervical cancer [3]. Consequently, the prognosis of patients who undergo surgery may deteriorate due to recurrence during this period. The study utilized aRMTLd analysis to examine the common conservative treatment approach of radiotherapy combined with chemotherapy, revealing a statistically significant improvement in patient prognosis. Notably, another study on SCCC, also derived from the SEER database, reported the benefits of surgical treatment, with no significant improvement in prognosis from radiotherapy [38]. This discrepancy with our study's conclusions may be attributed to the use of overall mortality as the endpoint, while our research included consideration of competing risks, focusing more on analyzing the treatment effects of different therapies on this specific cancer, with the primary event being death due to cervical cancer.

Furthermore, in practical applications, estimating the RMTL requires specifying fixed time points for analysis in advance. However, in clinical settings, it is sometimes challenging to rely solely on determining a single time point as a sufficient summary of treatment effects [39]. Moreover, once a specific time point is determined, subsequent information may not be fully utilized [31]. From this perspective, this study further presents the dynamic aRMTLd curve over a specified time interval. For comprehensive observation of the entire curve, simultaneous confidence bands (rather than pointwise confidence intervals) are more appropriate. Compared to observing differences at a specific time, examining the entire difference curve in conjunction with confidence bands provides more information about treatment benefits. For instance, analyzing radiotherapy at fixed time points provides an estimate that, within 120 months, radiotherapy can increase the life expectancy of patients free of cervical cancer death by approximately 7.5 months. However, further analysis incorporating dynamic curves revealed that the improvement in patient prognosis with radiotherapy was not significant between 3 and 16 months. However, the effect is statistically significant starting at 17 months, while it weakens in the later period (after 57 months) and is not significant. We speculate that this may be due to the occurrence of delayed treatment effects in the early stages, and when the treatment has taken effect, the efficacy may gradually diminish over time.

Limitations and future work

Finally, there are some limitations to this study: (1) First, the method proposed in this study can only handle baseline covariates. For time-dependent confounding factors, we can consider modifying the weight construction, which will be the focus of our future research. (2) Secondly, retrospective observational studies, when selecting databases for analysis, need to consider potential selection bias; (3) Furthermore, due to the lack of detailed information on treatment regimens in the SEER database, patient survival might be indirectly influenced.

Conclusion

In conclusion, our study introduces the aRMTLd as an estimation of treatment effects in the context of competing risk data, allowing for covariate adjustment. Simulation studies demonstrate the robust statistical performance of the proposed estimator. In the illustrative analysis, we utilized long-term follow-up data from patients with SCCC to investigate the impact of different treatments on the outcome of death due to cervical cancer. The aRMTLd curves are employed to depict the changing trends of treatment effects over time. Because this rare disease lacks clear treatment guidelines, the analytical findings regarding efficacy could guide future prospective studies. Moreover, given the results of this study, we can recommend the application of the aRMTLd and dynamic curve to competing risk data from observational studies to provide a reference for clinical and public health policy decisions in practical applications.