A nonparametric test for equality of survival medians using right-censored prevalent cohort survival data

Abstract

The median is a robust summary commonly used for comparison between populations. The existing literature falls short in testing for equality of survival medians when the collected data do not form representative samples from their respective target populations and are subject to right censoring. Such data commonly occur in prevalent cohort studies with follow-up. We consider a particular case where the disease under study is stable, that is, the incidence rate of the disease is stable. It is known that survival data collected on diseased cases, when the disease under study is stable, form a length-biased sample from the target population. We fill the gap for the particular case of length-biased right-censored survival data by proposing a large-sample test using the nonparametric maximum likelihood estimator of the survivor function in the target population. The small sample performance of the proposed test statistic is studied via simulation. We apply the proposed method to test for differences in survival medians of Alzheimer’s disease and dementia groups using the survival data collected as part of the Canadian Study of Health and Aging.

1 Introduction

In studies examining a survival time outcome for multiple groups/strata, it is typically of interest to the researcher to compare quantities related to the survival distribution such as the median. For example, Rahbar et al. compared the median waiting times for trauma patients to receive red blood cells using data from the Prospective Observational Multi-Center Major Trauma Transfusion study and, Brookmeyer and Crowley compared the median survival times for subjects receiving four types of cancer treatment using data from a Phase III colorectal cancer clinical trial.^1,2 Hypothesis testing procedures based on comparing the median survival times between groups, for a variety of data types, have been well studied in the literature. For right-censored failure time data, Brookmeyer and Crowley proposed a testing procedure for comparing $k$ sample medians using the asymptotic properties of the Kaplan–Meier product-limit estimator.^1,3 Although their procedure was based on comparing median point estimates, they implicitly assumed that the failure time distributions are members of the same location family. Brookmeyer and Crowley defined a weighted Kaplan–Meier estimate to determine the sample median, however, the non-parametric estimate assumes the group survival distributions are equivalent. Given that members of a location family can be identified by their medians, these authors essentially test the equality of the underlying survival distributions. Rahbar et al. highlight this testing dichotomy by explaining that equality of distributions implies equality of medians whereas the converse statement is not necessarily true. ⁴ Reid validated the bootstrapping procedure of Efron applied to the sample median for right-censored data and discussed the associated confidence interval properties.^5,6

Recently, Rahbar et al. proposed a hypothesis testing procedure for comparing population medians using right-censored failure time data which does not require the location family assumption. ⁴ Using the asymptotic distribution of the Kaplan–Meier estimator and an application of the functional delta method, Rahbar et al. derived a quadratic test statistic with an asymptotically chi-squared distribution. ⁷ In this article, we propose an extension to the method of Rahbar et al. for the case in which the failure time data are both length-biased and right-censored using the nonparametric maximum likelihood estimator (NPMLE) of the unbiased survival function. ⁸ We note that when the observed failure times are length-biased, the censoring mechanism is informative yielding an NPMLE which places mass on both the observed failure/censoring times of the prevalent cohort study.^8–11 Thus, the unbiased survival function estimator is fully defined over the entire real line establishing the existence of the median point estimate. While other product-limit estimators for left-truncated right-censored failure time data may be applied (see, e.g. Huang and Qin, ¹² Luo and Tsai, ¹³ and Tsai et al. ¹⁴ ), the survival function estimates may not yield a median point estimate when censoring is high. Although there may be instances in which the median point estimate exists when utilizing the bootstrapping procedure of Efron, some bootstrap samples may not yield median estimates requiring them to be discarded. ⁵ In addition to ensuring the existence of the median point estimate, our use of the NPMLE of the unbiased survival function is motivated by its efficiency compared to all other nonparametric estimators. ⁹

In Section 2., we define the notation for prevalent cohort studies with follow-up resulting in length-biased right-censored failure time data and describe how the asymptotic representation of the NPMLE of the unbiased survival function may be utilized to define a statistic for testing differences between multiple population medians. In Section 3., we generate simulated failure time data from prevalent cohort studies with follow-up to assess the performance of the proposed testing procedure and in Section 4., we apply our testing methodology to the Canadian Study of Health and Aging (CSHA) dataset to compare median survival between various dementia/Alzheimer’s disease groups. We provide some extensions and future research directions in Section 5.

2 Notation and methodology

In a general prevalent cohort study with follow-up, subjects are screened and only diseased cases are followed until failure/censoring for which their respective onset dates are recorded retrospectively. ¹⁵ Let the fixed constant $R$ denote the date at which subjects are tested for the prevalence of the condition/disease under study (prevalence day) and let $O<R$ denote a generic onset time generated from a stationary Poisson point process. Let $T$ be the underlying failure random variable with distribution function $F_U(\cdot )=1-S_U(\cdot )$ and mean $\mu$ . A subject is included in the prevalent cohort study if $T>R-O=A$ (i.e. they survive past prevalence day) where they are subsequently followed until the end of the study or they are lost to follow-up. Let $C$ denote the censoring time measured from $R$ with distribution function $G(\cdot )$ . Thus, the prevalent cohort consists of the triples $\{(A_i,Y_i,{\delta }_i)=(R-O_i,min(T_i,C_i),1_{\{T_i\le C_i\}}):i=1,2,\dots ,n\}$ made up of the observed truncation times, failure/censoring times and associated indicator functions. When the additional stationarity assumption is made for the underlying onset process of the prevalent cohort data (i.e. the underlying onset dates are drawn from a stationary Poisson point process) the observed failure times are typically referred to as being “length-biased.” For details on testing the validity of the stationary onset process assumption, see the literature.^16–19 Let $F_{LB}(\cdot )=1-S_{LB}(\cdot )$ denote the length-biased distribution function for which it may be shown ⁸ :

$$S_U(y)=\frac{{\int }_y^{\mathrm{\infty }}\frac{1}{x}d{S}_{LB}(x)}{{\int }_0^{\mathrm{\infty }}\frac{1}{x}d{S}_{LB}(x)};y>0$$

(1)

In this setting, it is sufficient to only use the observed failure/censoring times without the observed truncation times in order to make inferences regarding the unbiased survival function. For a set of length-biased right-censored failure time data, Asgharian and Wolfson established uniform consistency and weak convergence for the NPMLE of ${\hat{S}}_{LB}$ and showed

$$\sqrt{n}({\hat{S}}_{LB}-S_{LB})\stackrel{D}{\to }U={\mathcal{F}}^{-1}(V)$$

where $U$ is a Gaussian process. ⁹ Using the relationship between the unbiased and length-biased survival functions given in equation (1), Asgharian et al. established uniform consistency and weak convergence for the NPMLE of the unbiased survival function ${\hat{S}}_U$ . ⁸ In particular, they established the following weak convergence result:

$$\sqrt{n}({\hat{S}}_U(y)-S_U(y))\stackrel{D}{\to }{U}^{*}(y)=\mu {\int }_0^{\mathrm{\infty }}{\mathcal{L}}_y(x)dU(x)$$

where ${\mathcal{L}}_y(x)=I_{[y,\mathrm{\infty }]}(x)-S_U(y)/x$ . As the operator acting on the limiting Gaussian process of ${\hat{S}}_{LB}$ is a bounded linear operator, it may then be shown that $U^{*}(y)$ is also a Gaussian process. ²⁰ ${}^{(\mathrm{p}.377)}$

Using similar notation as Rahbar et al., consider the setting in which $k\ge 2$ independent populations are sampled for which the observed samples consist of length-biased failure/censoring times. ⁴ We utilize the subscript $j$ , $j=1,\dots ,k$ , to differentiate the parameters/estimators between the $k$ populations/samples, respectively. For $j=1,2,\dots ,k$ , we define the median as:

$${\theta }_{U,j}=inf\{t:S_{U,j}(t)\le 0.5\}$$

(2)

for which the sample median estimator, ${\hat{\theta }}_{U,j}$ is determined via plugging in the estimator of ${\hat{S}}_{U,j}$ in equation (2). ⁴ We note that the median estimators for the different groups will always exist as the corresponding unbiased survival function estimators, ${\hat{S}}_{U,j}$ , under the stationarity assumption, are fully defined over the entire real line. This is a direct result of the informative censoring mechanism which forces non-zero probability mass to fall on the observed censoring times. ²¹ Moreover, as the observed length-biased data are collected from a prevalent cohort study with follow-up through a cross-section at prevalence day, the subjects’ onset dates may occur any time arbitrarily far before prevalence day. This implies that if the follow-up period is short, resulting in a large proportion of the observed failure times being right-censored, the NPMLE of the $S_{U,j}$ functions will be unaffected. For details on the NPMLE of $S_{U,j}$ without the stationarity assumption, see Tsai et al. ¹⁴ For group $j$ , $j=1,2,\dots ,k$ , it may be shown that the sample median estimator converges to a Gaussian random variable:

$$\sqrt{n_j}\{{\hat{\theta }}_{U,j}-{\theta }_{U,j}\}\to Z_j\sim N(0,{\sigma }_j^2)$$

(3)

We summarize the necessary conditions of this result with a sketch of its proof in Theorem 1 and Corollary 4 in the Appendix. Estimates for ${\sigma }_j^2$ , $j=1,2,\dots ,k$ may be obtained by applying the simple bootstrap procedure of Efron to each of the $k$ samples by resampling the pairs $(Y_i,{\delta }_i)$ , $i=1,2,\dots ,n_j$ , $j=1,2,\dots ,k$ . ⁵

Let ${\theta }_{U,1},\dots ,{\theta }_{U,k}$ be the population medians and let $n_1,\dots ,n_k$ be the sample sizes for each of the $k$ samples. Consider the null hypothesis of equal medians for each of the $k$ groups:

$$H_0:{\theta }_{U,1}=\cdots ={\theta }_{U,k}(={\theta }_0)$$

Following the approach of Rahbar et al., we define the combined median estimator as:

$${\hat{\mathit{\theta }}}^c=({\hat{\theta }}_N^c,\dots ,{\hat{\theta }}_N^c)\mathrm{where}{\hat{\theta }}_N^c=\sum _{\mathrm{j}=1}^k\frac{n_j/{\hat{\sigma }}_j^2}{({\mathrm{n}}_1/{\hat{\sigma }}_1^2)+\cdots +({\mathrm{n}}_{\mathrm{k}}/{\hat{\sigma }}_{\mathrm{k}}^2)}{\hat{\theta }}_j$$

which can be re-expressed as

$${\hat{\mathit{\theta }}}^c=\sum _{j=1}^k{\hat{L}}_j{\hat{\theta }}_j$$

where

$${\hat{L}}_j=\frac{n_j/{\hat{\sigma }}_j^2}{(n_1/{\hat{\sigma }}_1^2)+\cdots +(n_k/{\hat{\sigma }}_k^2)}$$

for which $\sum _{j=1}^k{\hat{L}}_j=1$ . ⁴ Denoting the combined sample size $N=n_1+\cdots +n_k$ , Rahbar et al. proposed the following quadratic test statistic:

$${\hat{\mathrm{\Lambda }}}_N=N(\hat{\mathit{\theta }}-{\hat{\mathit{\theta }}}^c{)}^{\prime }{\hat{\mathrm{\Gamma }}}_N^{-1}(\hat{\mathit{\theta }}-{\hat{\mathit{\theta }}}^c)$$

(4)

where ${\hat{\mathrm{\Gamma }}}_N$ is a $k\times k$ matrix with diagonal entries given by

$${\hat{\mathrm{\Gamma }}}_{Ni}={\hat{\sigma }}_i^2{w}_{i,N}^{-1}(1-{\hat{L}}_i{)}^2+\sum _{j\ne i}^k{\hat{\sigma }}_j^2{w}_{j,N}^{-1}{\hat{L}}_j^2;i=1,2,\dots ,k$$

and off-diagonal entries given by

$${\hat{\mathrm{\Gamma }}}_{Nij}=-{\hat{L}}_i{w}_{i,N}^{-1}{\hat{\sigma }}_i^2-{\hat{L}}_j{w}_{j,N}^{-1}{\hat{\sigma }}_j^2+\sum _{m=1}^k{\hat{\sigma }}_m^2{w}_{m,N}^{-1}{\hat{L}}_m^2;i\ne j=1,2,\dots ,k$$

for which $w_{i,N}=n_i/N$ for $i=1,2,\dots ,k$ . ⁴ As remarked by Rahbar et al., the matrix ${\mathrm{\Gamma }}_N$ is not necessarily invertible and so the generalized inverse may be utilized instead of expression (4). Using the asymptotic normality of the sample median estimators, ${\hat{\theta }}_{U,j}$ , $j=1,2,\dots ,k$ , ${\mathrm{\Lambda }}_N$ may be shown to weakly converge to a central ${\chi }^2$ distribution with $k-1$ degrees of freedom (see Theorem 2 of Rahbar et al. ⁴ ). The percentiles of the central ${\chi }^2$ distribution may be utilized to test for significance when conducting the equality of medians hypothesis test.

2.1 Length-biased/unbiased survival function median comparisons

Although the proposed testing procedure is based on determining differences between the unbiased survival function medians, a question that may be posed is whether the correct conclusion will be reached when testing for equality of the length-biased population medians instead. Using the weak convergence of the length-biased survival function NPMLE ${\hat{S}}_{LB}(\cdot )$ , an equivalent form of the test statistic in equation (4) may be formed to test the null hypothesis:

$$H_0:{\theta }_{LB,1}={\theta }_{LB,2}=\cdots ={\theta }_{LB,k}(={\theta }_0^{\prime })$$

where ${\theta }_{LB,j},j=1,2,\dots ,k$ are the population medians of the length-biased survival functions. For simplicity of exposition, consider the case with $k=2$ in which the medians of two survival distributions are being compared. Addona examined different types of group comparisons between unbiased and length-biased survival curves and found that in the two-sample setting, there exist examples in which ${\theta }_{LB,1}<{\theta }_{LB,2}$ , however, ${\theta }_{U,1}>{\theta }_{U,2}$ (i.e. the ordering of medians is “reversed”). ²² In Figure 1, we illustrate a similar phenomenon using Weibull failure time data and show that the unbiased survival curves for two groups can have identical medians whereas the medians of the corresponding length-biased survival curves are clearly different. Therefore, without restriction on the family of underlying failure time distributions, it is not sufficient to test for differences in medians of the length-biased distributions and then subsequently make conclusions regarding the unbiased population medians. For related discussions on testing for equality of survival curves using length-biased failure time and adjusting for time-dependent selection mechanism/covariates, see Ning et al. ¹⁰ and Huiart and Sylvestre. ²³

Figure 1.

Illustration comparing unbiased and length-biased survival curves for two groups. The two unbiased survival curves were generated from Weibull distributions with (shape, scale) parameters of $(2,2)$ and $(1,2\log (2{)}^{(1/2)}/\log (2))$ yielding equal unbiased medians indicated by the single vertical line in panel (a). The length-biased medians, computed using the length-biased survival functions, are indicated by the two distinct vertical lines in panel (b).

3 Simulations

To assess the performance of the proposed hypothesis test, we generated samples of length-biased right-censored failure time data from various distributions to control for the values of the unbiased population medians. Specifically, we considered three failure time distributions: uniform distribution with support $(0,2\theta )$ (median equal to $\theta$ ), the LogNormal distribution with parameters $\mu$ and $1$ (median equal to $e^{\mu }$ ) and the exponential distribution with rate parameter $\log (2)/\lambda$ (median equal to $\lambda$ ). To generate a sample of length-biased right-censored failure time data, we sampled an onset time, $O$ , from a uniform distribution over $(0,R)$ where $R$ was chosen arbitrarily large such that $P(T>R)\approx 0$ . We then sampled a failure time $T$ and retained both the onset, failure time pair if $T>R-O$ , otherwise both sampled quantities were discarded. This procedure of accepting/rejecting the sampled (onset, failure time) pairs simulates the effect of left-truncation on the underlying failure time random variables. We censored the forward failure times $T-(R-O)$ by drawing an exponentially distributed censoring time, $C$ , with rate parameter $\varphi$ where we recorded $min(C,T-(R-O))$ and whether the minimum was a failure/censoring time. The rate parameter $\varphi$ was varied to allow for either $15\text{\%}$ or $30\text{\%}$ censoring.

We generated multiple length-biased right-censored failure time data sets and computed the empirical test sizes and empirical powers when the medians of the distributions were equal and different, respectively. We considered the case in which $k=3$ where we allowed the underlying failure time distributions to be equal or different from each other. We set the sample sizes of the three groups equal to $150$ , $200$ , and $250$ , respectively, and performed the same simulations when all three samples were of size $50$ . We ran $500$ simulation runs and for each run, we estimated the sample median variances using $1000$ bootstrapped samples. We reported the empirical test sizes in Table 1 and the empirical powers in Table 2. We find that for the majority of cases when the sample sizes are large, the empirical test sizes are similar to the true test significance level of $\alpha =0.05$ . When the sample sizes of the three groups were small, the empirical sizes increased in most cases. In the cases involving the exponential distribution, we attribute the larger empirical test sizes to the difficulty in drawing a sufficient number of observations arbitrarily close to time $t=0$ . In Table 2, we find that when the sample sizes are large, the empirical powers are in most cases above $0.50$ and when the sample sizes are small, the empirical powers are diminished. This phenomenon is to be expected as with smaller sizes, the variability of the estimators’ increases and it becomes more challenging to detect any differences between the underlying population medians.

Table 1.

Empirical test sizes for simulated length-biased right-censored failure times of three groups drawn from underlying survival distributions with equal medians of $\theta =10$ . Large sample sizes of three groups: $n_1=150$ , $n_2=200$ , $n_3=250$ . Small sample sizes of three groups: $n_1=n_2=n_3=50$ . Simulation runs: 500. Number of bootstrapped samples per run: 1000.

	Large samples		Small samples
	Censoring percentage
Group distribution types	15%	30%	15%	30%
All Uniform	0.050	0.040	0.110	0.088
All LogNormal	0.060	0.052	0.066	0.064
All Exponential	0.100	0.100	0.122	0.116
2 Uniform, 1 LogNormal	0.056	0.048	0.106	0.096
2 Uniform, 1 Exponential	0.080	0.070	0.120	0.108
2 LogNormal, 1 Uniform	0.048	0.058	0.066	0.076
2 LogNormal, 1 Exponential	0.084	0.086	0.088	0.082
2 Exponential, 1 Uniform	0.086	0.078	0.134	0.126
2 Exponential, 1 LogNormal	0.086	0.082	0.126	0.112
1 Uniform, 1 LogNormal, 1 Exponential	0.084	0.082	0.106	0.104

Table 2.

Empirical powers of a hypothesis test for simulated length-biased right-censored failure times of three groups drawn from underlying survival distributions with varying medians. Large sample sizes of three groups: $n_1=150$ , $n_2=200$ , and $n_3=250$ . Small sample sizes of three groups: $n_1=n_2=n_3=50$ . Simulation runs: 500. Number of bootstrapped samples per run: 1000.

	Large samples		Small samples
	Censoring percentage
Group distribution types and group medians $({\theta }_1,{\theta }_2,{\theta }_3)$	15%	30%	15%	30%
All Uniform $(10,10,6)$	0.728	0.690	0.718	0.656
All LogNormal $(10,10,e^2)$	0.364	0.368	0.102	0.110
All Exponential $(10,10,20)$	0.690	0.674	0.342	0.314
2 Uniform, 1 LogNormal $(10,10,e^2)$	0.544	0.518	0.220	0.200
2 Uniform, 1 Exponential $(10,10,20)$	0.742	0.732	0.354	0.338
2 LogNormal, 1 Uniform $(10,10,6)$	0.660	0.588	0.636	0.560
2 LogNormal, 1 Exponential $(10,10,20)$	0.736	0.726	0.352	0.306
2 Exponential, 1 Uniform $(10,10,6)$	0.664	0.620	0.658	0.598
2 Exponential, 1 LogNormal $(10,10,e^2)$	0.428	0.408	0.212	0.210
1 Uniform, 1 LogNormal, 1 Exponential $(10,e^2,20)$	0.896	0.896	0.496	0.462

To assess how the power of the hypothesis test varied depending on the group sample sizes, we simulated length-biased right-censored failure time data from three groups with different medians (Uniform: 10, LogNormal: $e^2$ , Exponential: 20). We varied the censoring variable’s rate parameter, $\varphi$ , so that either $15\text{\%}$ or $30\text{\%}$ of the sampled length-biased failure times were right-censored. We computed the empirical powers over $500$ simulation runs with $1000$ bootstrapped samples per run where all three sample sizes were equal to $50$ , $100$ , $200$ , $300$ , $400$ , $500$ , and $1000$ . The empirical power curves are plotted in Figure 2. Generally, as the sample size of the individual groups’ increases, the power of the test converges to $1$ independently of the two different proportions of censoring between the three groups.

Figure 2.

Empirical power curves of the median comparison hypothesis test for simulated length-biased right-censored failure data of three groups drawn from Uniform, LogNormal, and Exponential distributions with respective medians of $10$ , $e^2$ , and $20$ . The empirical powers were computed when the sample sizes of the three groups $n_1=n_2=n_3$ were equal to $50$ , $100$ , $200$ , $300$ , $400$ , $500$ , and $1000$ , where the proportion of censoring in the three groups was equal to $15\text{\%}$ (black) or $30\text{\%}$ (red). Simulation runs: $500$ . Number of bootstrapped samples per run: $1000$ .

4 Application

The CSHA was a multicenter study of dementia and related geriatric conditions for individuals over the age of 65 in Canada. ²⁴ In 1991, in the first stage of the CSHA, approximately 10,000 subjects living in a community or institutional residences were screened for various forms of cognitive impairment. At the second stage of the CSHA in 1996, the subjects who were still alive were screened a second time for cognitive impairment through a clinical evaluation. For the subjects who screened positive at the first stage of the CSHA, their onset dates were obtained through the recollections of their family members and caregivers and the dates of death were recorded for those subjects who died between the first and second stages of the CSHA. Through the screening tests and the clinical evaluations, the prevalent cohort of subjects ( $n=823$ ) consisting of $240$ males and $583$ females who screened positive at the first stage of the CSHA were classified as having either vascular dementia ( $n_1=173$ ), possible Alzheimer’s disease ( $n_2=252$ ) or probable Alzheimer’s disease ( $n_3=396$ ). Subjects’ survival times were (administratively) right-censored if they were still alive at the time of the second stage of the CSHA. Approximately $19\text{\%}$ , $24\text{\%}$ , and $21\text{\%}$ of the subjects screened with vascular dementia, possible Alzheimer’s disease, or probable Alzheimer’s disease, respectively, were right-censored at the second stage of the CSHA.

Using the observed failure/censoring times for the three dementia/Alzheimer’s disease groups of the CSHA, we calculated the unbiased survival function estimates and the associated unbiased survival function median estimates (see Figure 3). The medians for vascular dementia, possible Alzheimer’s disease, and probable Alzheimer’s disease groups (in months) were $39.9$ , $50.9$ , and $47.7$ , respectively. For other median estimates and estimates of the hazard ratio for death, stratifying by sex, education level, and age at onset of dementia, see Wolfson et al. ²⁴ Using the simple bootstrap, we generated 10,000 bootstrap samples for each group and determined estimates of the sample median variances. Setting the significance level $\alpha =0.05$ , using the bootstrapped estimates, we calculated the test statistic from equation (4), ${\mathrm{\Lambda }}_n=2.988$ , and did not reject the null hypothesis that individual dementia/Alzheimer’s disease group medians are equal.

Figure 3.

Nonparametric maximum likelihood survival function estimates for vascular dementia, possible Alzheimer’s disease, and probable Alzheimer’s disease groups in the Canadian Study of Health and Aging. Corresponding pointwise 95% confidence intervals are indicated with dotted lines and the median estimates for each group are indicated by dashed lines.

5 Discussion

The median estimator is a key summary statistic for determining the center of a failure time distribution and can be utilized to compare the failure time distributions of multiple independent groups. Using length-biased right-censored failure time data collected from a prevalent cohort study with follow-up, we proposed an extension to the Rahbar et al. testing procedure using the NPMLE of the unbiased survival function to compare the unbiased medians of $k\ge 2$ groups/populations. We applied our testing procedure to simulated prevalent cohort data sets and to data collected from the CSHA.

Throughout this article, we proposed a testing procedure for comparing the medians of multiple groups/populations using length-biased right-censored failure time; however, a closely related function to the median survival time is the quantile residual lifetime function (qrlt function). The qrlt function ${\theta }_{\alpha }(t)$ satisfies the equality $P(T-t\ge {\theta }_{\alpha }(t)|T\ge t)=\alpha$ for all $t\ge 0$ and some $\alpha \in [0,1]$ and determines the conditional $\alpha$ quantile when $T\ge t$ . The median survival time may be regarded as a single point of the qrlt function when $\alpha =0.5$ and $t=0$ . Recently, Liu et al. utilized the asymptotic properties of the unbiased survival function NPMLE to develop confidence intervals for two test statistics based on the ratio and difference of two qrlt functions. ¹¹ It is not immediately clear how to utilize the proposed test statistics in practice when performing a hypothesis test of equality of qltr functions. It remains an open problem how the proposed methodologies may be extended to the $k>2$ sample setting in which multiple population qrlt functions are being compared.