Nonparametric estimation of median survival times with applications to multi-site or multi-center studies

We propose a nonparametric shrinkage estimator for the median survival times from several independent samples of right-censored data, which combines the samples and hypothesis information to improve the efficiency. We compare efficiency of the proposed shrinkage estimation procedure to unrestricted estimator and combined estimator through extensive simulation studies. Our results indicate that performance of these estimators depends on the strength of homogeneity of the medians. When homogeneity holds, the combined estimator is the most efficient estimator. However, it becomes inconsistent when homogeneity fails. On the other hand, the proposed shrinkage estimator remains efficient. Its efficiency decreases as the equality of the survival medians is deviated, but is expected to be as good as or equal to the unrestricted estimator. Our simulation studies also indicate that the proposed shrinkage estimator is robust to moderate levels of censoring. We demonstrate application of these methods to estimating median time for trauma patients to receive red blood cells in the Prospective Observational Multi-center Major Trauma Transfusion (PROMMTT) study.


Introduction
Multi-site and multi-center studies have become very popular in clinical and translational sciences in the past two decades. Although multi-site studies allow for increased enrollment rate and improved generalizability to the target population [1], they introduce additional statistical challenges in the study design and analysis of data from these studies. For  studies could lead to non-homogeneous sub-samples due to the differences in the study sites. This is particularly relevant to trauma and emergency care research.
Having served as the Data Coordinating Center for the Prospective Observational Multicenter Major Trauma Transfusion (PROMMTT) study [1], we have identified a number of challenges in analysis of time to event data from PROMMTT. For example, in blood transfusion research, time to receive the first unit of red blood cells (RBCs) is one of important surrogates to measure how rapidly trauma patients receive blood transfusion. However, such data were often collected from multi-sites. Due to the fact that trauma centers may have their own guideline and practice to manage trauma patients, different sites may not only contribute different number of patients but also have different distributions of time to event of interest with different levels of censoring and variability. As a result, analysis of data from multi-site studies requires an exploratory step prior to pooling samples from different sites.
In some medical research, particularly in multi-site randomized clinical trials (e.g., for trauma/blood transfusion), when the main outcome of interest is time to observe a certain outcome of interest (e.g., time to receive the first unit of RBCs), complete observations in all patients are not usually available due to death before receiving a fixed number of units of RBCs, which results in right-censored data. [2] introduced nonparametric estimation of the survival curve, based on right-censored data. Nonparametric estimation of the mean survival time has been studied by many investigators including [3,4] and [5,6], and [7] have extensive discussions on nonparametric estimation of the mean and quantiles of the survival function. [8] introduced a nonparametric procedure for testing the equality of median survival times from k-independent samples of right censored data. [9] introduced an alternative to the BC test, which may perform better than the BC test under certain situations. To deal with the inflated Type I error rates when sample sizes are small, more recently, [10] extended Mood's median test for uncensored data to the setting of survival data.
The main aim of this research is to propose an improved nonparametric method for estimation of the median survival time from right-censored data from k-independent samples when there is uncertainty regarding the homogeneity of the k-population medians. We compare the proposed estimator to other two commonly used estimators asymptotically, and through extensive simulations. In addition, we demonstrate application of this method to data from the PROMMTT study ( [11,12]).
The remainder of the article is organized as follows. In Section 2, we propose the estimation strategies for median survival times and make comparisons among different estimators. In Section 3, we present the results from our simulation study comparing the performance of the proposed estimator against the other two commonly used estimators based on mean square errors (MSE). In Section 4, we demonstrate an application of our proposed method to data from PROMMTT ([1, 11]). Section 5 is devoted to concluding remarks with some discussion.

Materials and methods
We consider the nonparametric estimation of median survival times based on right-censored data in the presence of uncertain prior information in a k-independent sample situation. We assume here that T i1 ; Á Á Á ; T in i ; ði ¼ 1; Á Á Á ; kÞ are random samples selected from k populations with n i observations taken from the i-th population. We will refer to T ij as to the survival time of j-th subject in the i-th population. Due to censoring at time C ij , the survival time or time-toevent may not be observable in some subjects. Therefore, for each subject, the data are recorded in the form (Y ij , δ ij ), j = 1, Á Á Á, n i , i = 1, Á Á Á, k, where Y ij = min(T ij , C ij ), δ ij = I(T ij C ij ), and I(Á) is the indicator function. We assume that random variables T ij and C ij are independent with continuous survival distributions F i (x) = P(T ij > x) and G i (y) = P(C ij > y), respectively.

Unrestricted estimator (UE)
A straightforward way to estimate the median-parameter vector can be defined as, whereŷ i ¼F À 1 i ð0:5Þ andF i ðpÞ ¼ inf ft :F i ðtÞ ! pg for i = 1, Á Á Á, k. We call this estimator as unrestricted estimator (UE) of Θ. This estimator is usually used when no hypothesis information is available on Θ. For example, in a multi-site study with k sites one can provide site-specific estimates for median survival time without combing the data from all sites. All the estimators for the k components of the vector are independent.
For each i = 1, Á Á Á, k, ffiffiffiffi n i p fF À 1 i ðpÞ À F À 1 i ðpÞg converges in distribution to a normal random variable with mean 0 and variance C i ðpÞ ¼ pð1 À pÞ=ff i ðF À 1 i ðpÞÞ 2 g. Hence operationally, F À 1 i ðpÞ is asymptotically normal with mean F À 1 i ðpÞ and variance C i (p)/n i .

Combined estimator (CE)
Θ can also be estimated by combining the sample and hypothesis information under the assumption of homogeneity of the k medians given by We can use this additional information together with the sample information to obtain improved estimators. Under the null hypothesis (2.2), we consider the combined/restricted estimator (CE) of Θ defined bŷ This estimator is expressed as a linear combination of theŷ i 's, i.e.,ŷ CE n ¼ i Þ=ðn 1 =ŝ 2 1 þ Á Á Á þ n k =ŝ 2 k Þ, whereŝ 2 i can be estimated by Var ðF i ðŷ i ÞÞ=f 2 i ðŷ i Þ for i = 1, Á Á Á, k. However, in order to avoid the difficulty in estimating the density function, we used the bootstrap to estimate the variance of the estimated median, following [9].
For the preliminary test on H 0 in (2.2), we consider the following test statistic that is defined by the normalized distance ofŶ fromŶ CE : We assume that ω i = lim n!1 ω i,n is fixed for i = 1, Á Á Á, k, and G ¼ lim n!1Ĝn exists and is nonsingular. It is shown that under the null hypothesis for large n, Λ n follows the central χ 2 distribution with k − 1 degrees of freedom [9]. For given α, the critical value of Λ n may be approximated by w 2 kÀ 1;a , the upper 100α% point of the chi-square distribution with k − 1 degrees of freedom. More details can be found in [9].

Positive-part shrinkage estimator (PP)
The combined estimator works well only when the null hypothesis (2.2) holds. If the null hypothesis (2.2) is rejected, we propose another estimator based on the James-Stein type shrinkage estimator (SSE) [12] which is defined bŷ The Stein-type estimator in (2.6) is not sensitive to departure from H 0 , and will provide uniform improvement in terms of efficiency for the entire parameter space of Θ. It is, however, not a convex combination ofŶ CE andŶ. Also, this estimator may not remain nonnegative. To avoid this strange behavior, we truncateŶ JS at positivity boundary by adding an extra term to (2.6), which leads to a convex combination ofŶ andŶ CE , namely, the positive-part shrinkage estimator (PP). When k ! 4, the positive-part shrinkage estimator is defined as follows: Y PP ¼Ŷ À ðk À 3ÞL À 1 n ðŶ ÀŶ CE Þ À f1 À ðk À 3ÞL À 1 n gIðL n < k À 3ÞðŶ ÀŶ CE Þ : ð2:7Þ whereŶ,Ŷ CE and Λ n were defined in Sections 2.1 and 2.2.

Comparison ofŶ UE ,Ŷ CE andŶ PP
In this Section, we compare the performance ofŶ UE ,Ŷ CE andŶ PP by asymptotic distribution quadratic risk function [13]. For an estimator Θ Ã , define the weighted quadratic loss function of the form is called the risk function. The performance of the estimators can be evaluated by comparing the risk functions and an estimator with a smaller risk is preferred.
Since the test statistic in (2.2) is consistent for fixed Θ when Θ = 2 H 0 ,Ŷ PP is asymptotically equivalent toŶ UE for fixed alternatives, this makes it difficult to compare their performance [14]. Alternatively, we may evaluate the asymptotic performance of each estimator under the following contiguous sequence of alternatives: where φ is a fixed vector and Then, considering the asymptotic distribution of ffiffiffi n p ðY Ã À YÞ, we can define the asymptotic distribution quadratic risk (ADQR) as R(Θ Ã , Θ) = tr(WΓ) where Γ is the asymptotic covariance matirx. To facilitate the numerical computation and general discussion, we consider the particular case with W = Γ − 1 . Then, following similar arguments in [14], and ð2:11Þ and , if the ADQR of Θ Ã is smaller for at least some value of Θ, and the ADQR does not exceed that of Θ 0 for any value of Θ. Further, . At Δ Ã = 0, that is, under the null, the dominance of the estimators is usually observed asŶ CE 1Ŷ UE , where the notation 1 stands for dominance in terms of risk performance. For all Δ Ã and k ! 4,Ŷ PP 1Ŷ UE is satisfied, that is,Ŷ PP asymptotically dominatesŶ UE under local alternatives. Thus, we conclude thatŶ PP consistently performs better thanŶ UE in the entire parameter space induced by Δ Ã . The gain in risk overŶ UE is substantial when Δ Ã = 0 or near 0.

Simulation studies
We conducted extensive simulation studies to examine the performance of the proposed estimators in situations with different degrees of departure from the assumption of homogeneity and censoring rates.
In order to evaluate the effect of the departure from the null hypothesis, we generated samples with median where η is a fixed value to achieve a desired level of censoring (e.g., p = 0.3). We set Θ 0 = (6, 6, 6, 6), a = (2, 2, 2, 2), and c = (1, 2, 1, 2). The simulation procedure was repeated 1000 times for k = 4 independent samples with size of 100 for each group.
The simulation results of REs and the comparative plots are presented in Table 1 and Fig 1. The performance of Θ Ã was measured by the relative efficiency (RE), i.e., comparing its MSE with that of Θ UE , defined as REðŶ UE ;Ŷ Ã Þ ¼ MSEðŶ UE Þ=MSEðŶ Ã Þ, whereŶ Ã is one of the estimators (Ŷ CE orŶ PP ) considered in this study. The amount by which a RE is larger than 1 indicates the degree of superiority of the estimatorŶ Ã overŶ UE . As highlighted in [9], to avoid difficulty in estimating the densitity function, we used boostrapping to estimate the variance of the estimated median. We also compare the performance ofŶ Ã by asymptotic distribution quadratic risk described in Section 2.4. Table 1 shows that when d = 0 or near zero,Ŷ CE outperforms all other estimators. However, as d moves away from 0, the risk ofŶ CE become unbounded, making it very inefficient.
Overall,Ŷ PP maintain its superiority over other estimators for a wide range of d . This clearly suggests thatŶ PP is preferred as there always remains uncertainty about level of heterogeneity between survival medians. In situations where the assumed model is grossly wrong,Ŷ PP is expected to be as good or equal toŶ UE . The asymptotic behavior of ADQR over d is illustrated in Fig 1 under different significance alpha levels. It shows thatŶ PP may has smaller risks and thus be more efficient uponŶ UE under the null or near. The ADQR ofŶ CE approaches infinity as d grows. Overall,Ŷ PP has a better performance in the entire parameter space.
Next, we study scenarios with different distributions and censoring rates. In each scenario, 50 or 100 samples were generated from k = 4 subpopulations. The distributions considered in Tables 2 and 3 are (i) a uniform distribution with median Θ 0 = (6,6,6.5,6), (ii) a log-normal distribution with Θ 0 = (403.43, 665.14, 403.43, 665.14), and (iii) an exponential distribution with Θ 0 = (6.69, 6.58, 7.0, 6.43). The distributions considered in Tables 4 and 5 are mixed distributions with (iv) two subpopulations had a uniform distribution and the other two had a log-normal distribution with median Θ 0 = (6, 6, 4.48, 7.39), (v) two uniform and two exponential distributions with median Θ 0 = (6, 6, 6.15, 6.15), and (vi) two log-normal and two exponential distributions with median Θ 0 = (6, 6, 6.89, 6.89). For each scenario, simulations were repeated 1000 times under 0%, 30%, and 50% censoring schemes and the results are summarized in Tables 2 to 5. We also considered the scenarios where different sites had different censoring rates (0% to 70%) in Table 6, using three mixed distributions and a sample size of 50. In general, the results in Tables 2 and 6 show that survival distribution and censoring rate affect the deviation from the true value for different estimators. The RE ofŶ CE ranges from

Application to the PROMMTT study
The PROMMTT study was a multi-site prospective observational cohort study in a severely injured transfused trauma patients, conducted at 10 level 1 trauma centers in the United States [1,11]. The original objectives of PROMMTT study were to accurately describe when some blood components were infused and to assess the association between in-hospital mortality and the timing and amount of blood products. Understanding current blood product usage patterns and their impact on patient outcomes among a severely injured and substantially hemorrhaging cohort is critically important. In our analysis, we applied our method on 698 patients from four major medical centers to examine difference of their median times of receiving the first unit of RBC infusion. The number of patients in the 4 sites are 303, 137, 133, and 125. The median age of patients are 34, 37, Table 5. Mean of 1000 estimated medians and REs for 0%, 30% and 50% censoring rates for (2 uniform+2 lognormal), (2 uniform+2 exponential), (2 lognormal +2 exponential) distributions when sample size is 100. Since no patients dropped out of study before the first unit of RBC infusion, the censoring rate for each site is zero. We apply the proposed method to this data in two steps. First, we formally evaluate homogeneity of medians using randomly right censored data. Second, based on the results of the test of homogeneity our proposed method evaluates whether the unrestricted estimator (UE), combined estimator (CE), or the Positive Part shrinkage estimator (PP) should be used to obtain a more efficient estimator. The primary outcome is time to the first unit of RBC infusion measured as the number of minutes from ED admission. This outcome may reflect patient's status and may differ across different medical institutions because management of severely hemorrhaging patients may differ. As we found in the data, the UE median times to the first unit of RBC infusion are 18, 55, 65, and 24 for four sites, indicating potential heterogeneity of survival times across 4 sites. Fig 2 shows the estimated survival curves of time to the first unit of RBCs for each site, along with p-values from (2.2) and the log-rank test. It shows that median times are significantly different, suggesting that caution should be exercised when merging data sets from different sites. All of these motivate us to apply the proposed Positive Part shrinkage estimator (PP) method.   unreliable since the test for homogeneity of medians to first RBC infusion is rejected (p-value<0.0001). The PP median times are almost identical to the UE median times, with slightly narrower confidence intervals than UE median times at site 2, 3, and 4, indicating a slight efficiency gain.

Concluding remarks
In this paper, other than the unrestricted estimator and simply combined estimator, we have presented a shrinkage nonparametric approach for estimating the median survival vector in a k-sample problem. Other than asymptotical comparison on ADQR, extensive simulations have been done to assess the performance of these estimators, considering various scenarios allowing varying levels of censorship and different level of departure from homogeneity of the survival medians.
The performance of the combined estimator heavily depends on the strength of homogeneity. When homogeneity holdsŶ CE is more efficient compared toŶ UE andŶ PP . However,Ŷ CE becomes inconsistent and the efficiency of the CE decreases significantly when homogeneity fails. On the other hand,Ŷ PP seems to be robust to the non-homogeneity and different levels of censoring. Though the relative efficiency againstŶ UE decreases as we deviate from the quality of the survival medians, it keeps being greater or equal to 1. Like any other shrinkage estimation procedures, there is a bias-variance tradeoff forŶ PP . The proposed procedures have applications in epidemiologic and health care research. For example, for estimating survival median based on data from a multi-site study one always faces with the challenge of whether to pool data from all sites or not pool such data. In this study using PROMMTT data we have demonstrated the utility of various estimators as well as how one can make a decision as to choose the most appropriate estimation procedure. On the other hand, if the distributions of median times are similar, greater efficiency of the estimators by using shrinkage-type methods may be gained, depending on the distribution of event time.
Supporting information S1 File. Revised PROMMTT publication committee guidelines. (DOCX) S1 Dataset. R code for generating simulated data and analysis. (ZIP)