Estimation of monotone bivariate quantile inactivity time with medical applications

1.   Introduction

Great progress has been made in the field of statistics and probability theory for interdisciplinary research. New techniques and methods have been developed to meet the challenges of data analysis. Statistical methods are becoming increasingly important in various areas of science. The increasing complexity of scientific problems requires the development of new and suitable statistical methods for interdisciplinary research. Current challenges include, ecology: Quantifying biodiversity; the epidemiology of infectious diseases: Disease outbreak detection; financial mathematics: stock option valuation; industrial engineering: Stochastic optimization; and genomics: Personalized medicine, to name a few.

For a random lifetime $T$, the conditional inactivity time is defined as ${T}_{t} = t-T|T\le t,$ $\mathrm{a}\mathrm{n}\mathrm{d}$ $t > 0$. An important measure developed based on the conditional inactivity time is the $\alpha$-quantile inactivity time ($\alpha$-QIT), which is the $\alpha$-quantile of ${T}_{t}$. Assuming that the distribution function of $T$ is denoted by $F$, the $\alpha$-QIT can be expressed by the following relationship:

where $\stackrel{-}{\alpha } = 1-\alpha$ and ${F}^{-1}\left(p\right) = \mathrm{i}\mathrm{n}\mathrm{f}\{x:F(x) = p\}$ is the inverse function of $F$. Let $T$ be the event time referring to the instances of a species. Among all instances experienced the event at a time before $t$, we expect $100\alpha \mathrm{\%}$ of these instances to have experienced the event after time $t-{q}_{\alpha }\left(t\right)$. In this sense, a smaller ${q}_{\alpha }\left(t\right)$ means larger $T$. The $\alpha$-QIT is a competitor for the mean inactivity time (MIT) function. The MIT has been intensively studied by researchers in the field of reliability theory and survival analysis, e.g., refer to Finkelstein ^[1] and Kayid and Izadkhah ^[2]. However, when the moments of the underlying model are infinite or heavily skewed to right, $\alpha$-QIT is preferred over MIT (see Schmittlein and Morrison ^[3] for a detailed justification of quantile-based than moment-based measures). The $\alpha$-QIT concept was formally defined and studied by Unnikrishnan and Vineshkumar ^[4]. Shafaei ^[5] showd that how the underlying model can be characterized by $\alpha$-QIT function. Shafaei and Izadkhah ^[6] stated some properties of a parallel system in terms of the $\alpha$-QIT measure. For a sample ${T}_{1}$, ${T}_{2}$, ..., ${T}_{n}$ of iid lifetimes, the $\alpha$-QIT can be estimated by

where ${F}_{n}\left(t\right)$ is the empirical distribution function, i.e.,

and

Then, ${q}_{\alpha, n}$ can be written as in the following.

where

Note that ${F}_{n}^{-1}\left(\stackrel{-}{\alpha }{F}_{n}\left({t}_{i}\right)\right) = {T}_{\left(i\right)}$, for every $i = \mathrm{1, 2}, \dots, k$ for some $k\le n$ ${t}_{1} = {T}_{\left(1\right)}$ and ${t}_{i}\ge {T}_{\left(i\right)}$ for $i = 2\dots..., k$. The expression $\left(1\right)$ shows that ${q}_{\alpha, n}\left(t\right)$ consists of line segments with slope 1 on intervals $({t}_{i}, {t}_{i+1})$, $i\dots, 2, ..., k-1$ and falls at each point ${t}_{i}$ by ${T}_{\left(i\right)}-{T}_{(i-1)}$, $i = \mathrm{1, 2}, ..., k$ where ${T}_{\left(0\right)} = 0$. Figure 1 shows a schematic plot of ${q}_{\alpha, n}\left(t\right)$.

For the univariate case, Mahdy ^[7] proposed the estimator (2) for the α-QIT function and investigated its asymptotic properties. Balmert and Jeong ^[8] created a nonparametric inference of the median inactivity time function for right-censored data. Balmert et al. ^[9] applied a log-linear quantile regression model to the inactivity time for right-censored data. Kayid ^[10] applied the Kaplan-Meiere survival estimator to the $\alpha$-QIT function for estimation and inference.

We can have two or more dependent events. For example, if successive events of the same person/instance are tracked, the event times depend on each other. Another example is that researchers are interested in determining the effect of a treatment on specific event times related to the eyes, ears, hands or legs. One organ was randomly selected for treatment and the other was a control organ. The events associated with these organs depend on their progression. In such cases, we need to extend the measures in question to bivariate or multivariate settings. In the following section, I refer to the authors who have implemented this idea. Basu ^[11] and Johnson and Kotz ^[12] examined the multivariate hazard rate function as a gradient vector. The mean residual lifetime was extended by Nair and Nair ^[13] to obtain a vector of dependent lifetimes. Shaked and Shanthikumar ^[14] introduced the dynamic multivariate MRL concept. Kayid ^[15] developed the multivariate MIT concept. Navarro ^[16] characterized the basic model by the bivariate hazard rate function. The concept of the $\alpha$-quantile residual lifetime ($\alpha$-QRL) was extended to the multivariate context by Shafaei and Kayid ^[17]. Shafaei et al. ^[18] discussed the multivariate $\alpha$-QRL concept in a dynamic way. Buono et al. ^[19] applied multivariate RHR for discussing reliability attributes of systems. Kayid ^[20] extended the $\alpha$-QIT concept to bivariate context and discussed its estimation.

Let $F$ be the distribution function of a random pair $\boldsymbol{T} = ({T}_{1}, {T}_{2})$. Then, the $\alpha$-QIT vector at point $\boldsymbol{t} = ({t}_{1}, {t}_{2})$ is defined to be $\left({q}_{\alpha, 1}(\boldsymbol{t}), {q}_{\alpha, 2}(\boldsymbol{t})\right)$. The first element of this vector is

where $\stackrel{-}{\alpha } = 1-\alpha$ and

is the partial inverse of $F$ in terms of the ${T}_{1}$. The second element of the $\alpha$-QIT vector is defined similarly.

where ${F}_{2}^{-1}\left(\mathrm{p}; {t}_{1}\right) = \mathrm{i}\mathrm{n}\mathrm{f}\{z:F({t}_{1}, \mathrm{z}) = p\}$ is the partial inverse of $F$ in terms of the second element. The reversed hazard rate vector of $\boldsymbol{T}$ is $\left({r}_{1}(\boldsymbol{t}), {r}_{2}(\boldsymbol{t})\right)$ and

The RHR satisfies the following relation.

Kayid ^[20] showed that if ${r}_{i}\left(\boldsymbol{t}\right)$ is decreasing (increasing) in ${t}_{i}$, then ${q}_{\alpha, i}\left(\boldsymbol{t}\right)$ is increasing (decreasing) in ${t}_{i}$. It is a surprising fact that for most of the standard bivariate models, ${r}_{i}\left(\boldsymbol{t}\right)$ is decreasing in ${t}_{i}$ (Finkelstein ^[1]). For example, bivariate Gumbel, Pareto, normal, and gamma models have decreasing reversed hazard rate functions. This implies that ${q}_{\alpha, i}\left(\boldsymbol{t}\right)$ is increasing in ${t}_{i}$. For some examples of such models, refer to Kayid ^[20]. This motivates me to introduce a new estimator of ${q}_{\alpha, 1}\left(\boldsymbol{t}\right)$ and ${q}_{\alpha, 2}\left(\boldsymbol{t}\right)$ under the assumption that they are increasing with respect to ${t}_{1}$ and ${t}_{2},$ respectively. It is expected that applying this knowledge, I have a more accurate estimator than the usual estimator defined by Kayid ^[20]. Such monotone estimators are defined and studied by Kochar et al. ^[21], Franco Pereira and Una-Alvarez ^[22], and Shafaei and Franco Pereira ^[23].

The rest of this paper is structured as follows. In Section 2, the promised increasing estimator of the bivariate $\alpha$-QIT function is proposed and its asymptotic properties are discussed. Then, the performance of the new estimator is compared with that of the usual estimator in a simulation study. In Section 4, the proposed estimator is applied to investigate the effect of laser treatment on the time to blindness. In Section 5, I summarize the results.

2.   Estimation of increasing bivariate α-QIT

Let $\dots, ..., {\boldsymbol{T}}_{n}$ be an iid random sample from bivariate distribution $F$. The empirical distribution function is defined by

and the partial inverse of ${F}_{n},$ with respect to the first and second elements, are as in the following respectively:

and

Kayid ^[20] proposed the following estimator of the bivariate $\alpha$-QIT vector.

where

with the knowledge of increasing bivariate $\alpha$-IQT, we define the natural estimator

where

Let ${t}_{2} > 0$ be fixed and define ${T}_{1}\left[{t}_{2}\right] = {T}_{1}|{T}_{2}\le {t}_{2}$, then the distribution function of ${T}_{1}\left[{t}_{2}\right]$ is

where ${F}_{2}\left({t}_{2}\right) = P({T}_{2}\le {t}_{2})$. Denote $\alpha$-QIT of ${T}_{1}\left[{t}_{2}\right]$ by ${q}_{\alpha, 1}^{\mathrm{*}}({t}_{1};{t}_{2})$, then it can be shown that

Similarly, for every fixed ${t}_{1} > 0$, we define ${T}_{2}\left[{t}_{1}\right] = {T}_{2}|{T}_{1}\le {t}_{1}$ following distribution ${F}_{2}^{\mathrm{*}}(.; {t}_{1})$. Let ${q}_{\alpha, 2}^{\mathrm{*}}({t}_{2};{t}_{1})$ be the $\alpha$-QIT of ${T}_{2}\left[{t}_{1}\right]$, then we can investigate that

Given a bivariate iid random sample $({T}_{1i}, {T}_{2i})$, $i = \mathrm{1, 2}, \dots, n$ from distribution $F$, and for every fixed ${t}_{2}$, consider the following univariate random sample which follows from ${F}_{1}^{\mathrm{*}}(.; {t}_{2})$.

I can apply this sample to estimate ${q}_{\alpha, 1}^{\mathrm{*}}({t}_{1};{t}_{2})$, as in the following.

where

and

Applying the knowledge of increasing ${q}_{\alpha, 1}^{\mathrm{*}}({t}_{1};{t}_{2})$ in terms of ${t}_{1}$, it is natural to use the following estimator.

Again, for a bivariate iid random sample $({T}_{1i}, {T}_{2i})$, $i\dots, 2, ..., n$ from distribution $F$, and for every fixed ${t}_{1}$, consider the following sample.

which follows from ${F}_{2}^{\mathrm{*}}(.; {t}_{1})$. Then, the estimator of ${q}_{\alpha, 2}^{\mathrm{*}}({t}_{2};{t}_{1})$ is defined by

where

and

In an increasing context,

It is clear that

Theorem 1. Let us assume that the following two conditions are fulfilled.

(C1). $F({t}_{1}, {t}_{2})$ be twice differentiable with respect each element.

(C2). $\frac{\partial }{\partial {t}_{1}}F\left(\boldsymbol{t}\right)$ and $\frac{\partial }{\partial {t}_{2}}F\left(\boldsymbol{t}\right)$ are bounded from zero on the intervals $(0, {F}_{1}^{-1}(\stackrel{-}{\alpha }; {t}_{2}))$ and $(0, {F}_{2}^{-1}(\stackrel{-}{\alpha }; {t}_{1}))$, respectively, for every ${t}_{1} > 0$ and ${t}_{2} > 0$.

Then, $\left(i{q}_{\alpha, 1, n}(\boldsymbol{t}), i{q}_{\alpha, 2, n}(t)\right)$ is consistent for $\left(i{q}_{\alpha, 1}(\boldsymbol{t}), i{q}_{\alpha, 2}(\boldsymbol{t})\right)$.

Proof. By Theorem 7 from Kayid ^[24], we have

and

Thus, the result follows from (2)–(4).

To state the next theorem, I need two following conditions.

(C3) $\frac{\partial }{\partial {t}_{1}}{q}_{\alpha, 1}\left(\boldsymbol{t}\right)$ and $\frac{\partial }{\partial {t}_{2}}{q}_{\alpha, 2}\left(\boldsymbol{t}\right)$ exist and there are ${c}_{1} > 0$ and ${c}_{2} > 0$ such that $\frac{\partial }{\partial {t}_{1}}{q}_{\alpha, 1}\left(\boldsymbol{t}\right) > {c}_{1}$ and $\frac{\partial }{\partial {t}_{2}}{q}_{\alpha, 2}\left(\boldsymbol{t}\right) > {c}_{2}$ for all $0 < {t}_{1} < {b}_{1}$ and $0 < {t}_{2} < {b}_{2}$ for some positive ${b}_{1}$ and ${b}_{2}$.

(C4) $\frac{{\partial }^{2}}{\partial {t}_{1}\partial {t}_{1}}{q}_{\alpha, 1}\left(\boldsymbol{t}\right)$ and $\frac{{\partial }^{2}}{\partial {t}_{2}\partial {t}_{2}}{q}_{\alpha, 2}\left(\boldsymbol{t}\right)$ exist and

Theorem 2. Assume that C1–C4 are satisfied. Then, we have

Proof. By Theorem 5 from Kayid ^[21], we have

and

Thus, the result follows from relations (2)–(4) and the concept of convergence in probability in bivariate setting.

The following lemma, which is the result of the well-known Slutsky theorem, is used in the proof of the next theorem (see Van der Vaart ^[25] for Slutsky's theorem and related results).

Lemma 1. If $\sqrt{n}({X}_{n}-{Y}_{n})\to 0$ in probability and $\sqrt{n}{X}_{n}$ converges, in distribution, to a random variable $X$ with distribution $F$, then $\sqrt{n}{Y}_{n}$ converges, in distribution, to a random variable $Y$ with the same distribution $F$.

Theorem 3. Under the conditions C1–C4, we have

where

and elements of $\mathrm{\Sigma }$ are

and

Proof. Theorem 7 of Kayid ^[20] states that under some mild conditions:

where $C$ and $\mathrm{\Sigma }$ are defined in this theorem. Thus, applying Lemma 1, the result follows immediately.

In the real world, lifetime random pairs ${\boldsymbol{T}}_{1}, {\boldsymbol{T}}_{2}, ..., {\boldsymbol{T}}_{n}$ may be censored by a random censorship ${C}_{i}$, in the sense that the observations are ${\stackrel{~}{T}}_{1i} = {T}_{1i}\wedge {C}_{i}$, ${\stackrel{~}{T}}_{2i} = {T}_{2i}\wedge {C}_{i}$, ${\delta }_{1i} = I({T}_{1i} > {C}_{i})$ and ${\delta }_{2i} = I({T}_{2i} > {C}_{i})$. Note that $a\wedge b = \mathrm{m}\mathrm{i}\mathrm{n}\{a, b\}$. Let censorship random variable ${C}_{i}$ be independent from desired lifetimes and follows from distribution $G$ and the reliability function $\stackrel{-}{G} = 1-G$, i.e., $\stackrel{-}{G}\left(t\right) = P({C}_{i} > t)$. Also, let $\stackrel{~}{R}({t}_{1}, {t}_{2}) = P({\stackrel{~}{T}}_{1i} > {t}_{1}, {\stackrel{~}{T}}_{2i} > {t}_{2})$ and $R({t}_{1}, {t}_{2}) = P({T}_{1i} > {t}_{1}, {T}_{2i} > {t}_{2})$. Then, we have

where ${t}_{1}\vee {t}_{2} = \mathrm{m}\mathrm{a}\mathrm{x}\{{t}_{1}, {t}_{2}\}$. So, we can estimate the reliability function $R$ by

Under this censoring scheme, Lin and Ying ^[26] showed that ${R}_{n}({t}_{1}, {t}_{2})$ is strongly consistent and weakly converges to a Gaussian process. Thus, when we have such censored data, the empirical distribution function could be replaced by the following estimate:

3.   Simulations

To investigate the performance of the proposed (increasing) estimator and comparing it with the usual estimator, a simulation study is conducted. The bivariate Gumbel and Pareto distributions with respectively the following reliability functions are selected for the baseline models:

and

Both models are important from practical and theoretical points of view. The Gumbel distribution was introduced by Gumbel ^[27], and the Pareto model was used by Jupp and Mardia ^[28] to analyze income data for consecutive years. Some proper values for $\beta$ and $c$ were selected. In each simulation run, $r = 1000$ replicates of bivariate samples of size $n$ were generated, where $n$ was set to $25$, $50$ or $100$. For each sample, ${q}_{\mathrm{0.5, 1}, n}\left(\right)$ and its increasing version, $i{q}_{\mathrm{0.5, 1}, n}\left(\right)$, are calculated at four appropriate time points ${t}_{1}$, ${t}_{2}$, ${t}_{3}$ and ${t}_{4}$ according to the following rules: Let ${F}_{1}$ be the marginal distribution of the first element and ${\boldsymbol{t}}_{i} = ({t}_{1i}, {t}_{2i})$. The equations ${F}_{1}\left({t}_{11}\right) = 0.25$, ${F}_{1}\left({t}_{12}\right) = 0.40$, ${F}_{1}\left({t}_{13}\right) = 0.50$ and ${F}_{1}\left({t}_{14}\right) = 0.75$ are solved to find ${t}_{11}$ to ${t}_{14}$ and given them, the equations $F({t}_{11}, {t}_{21}) = 0.2$, $F({t}_{12}, {t}_{22}) = 0.3$, $F({t}_{13}, {t}_{23}) = 0.4,$ and $F({t}_{14}, {t}_{24}) = 0.6$ are solved for ${t}_{21}$ to ${t}_{24}$. After calculating the objective functions for $r$ replicates, the bias (B) and mean squared error (MSE) were calculated and are shown in Tables 1 and 2 for the Gumbel and Pareto models, respectively. All simulations and calculations were performed in R (statistical programming language). The results show small values for B and MSE for both the conventional estimator and the proposed increasing estimator. As expected, the MSE increases with $F\left(\mathrm{t}\right)$ (see Theorem 3). The MSE values for the increasing estimator are smaller in all cases, indicating that the increasing estimator performs better than the conventional estimator. See Figures 2 and 3 for a graphicall representaion of the ratio of MSE values related to the ususal to the increasing estimator.

4.   Effect of laser treatment on blindness

In a study that began in 1971, researchers were interested in the effect of laser photocoagulation on delaying blindness in patients with DR. Patients with visual acuity ≥ 20/100 in both eyes were selected for the study. One eye of each patient was randomly selected for laser photocoagulation (treatment) and the other eye was observed without treatment (control). The time from the start of treatment to blindness is given in months. Blindness means that visual acuity fell below 5/200 on two consecutive visits. The data for this study is available in the "diabetic" dataset in the "survival" package in R. Table 3 shows part of the dataset relating to adolescents (under 20 years of age). For patient $i$, ${T}_{1i}$ and ${T}_{2i}$indicate the observed time to blindness in the control and treated eyes, respectively.

Figures 4 and 5 draw the bivariate median inactivity time functions and their increasing versions, $i{q}_{n, \mathrm{0.5, 1}}\left(\boldsymbol{t}\right)$ and $i{q}_{n, \mathrm{0.5, 2}}\left(\boldsymbol{t}\right), \ \ \mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{y}$.

The comparison of the proposed increasing median inactivity time functions $i{q}_{n, \mathrm{0.5, 1}}\left(\boldsymbol{t}\right)$ and $i{q}_{n, \mathrm{0.5, 2}}\left(\boldsymbol{t}\right)$ at different points is informative in investigating the treatment effect. To provide a simple and powerful statistics, we can consider the points on the identity line and use the following statistics

If I assume that the treatment dose not effect the time length to blindness, ${d}_{n}\left(t\right)$ should be positive or negitive values near zero, reflecting some random errors. However, if the treatment causes longer time to blindness, I expect relatively larger values for $i{q}_{n, \mathrm{0.5, 1}}(t, t)$ than $i{q}_{n, \mathrm{0.5, 2}}(t, t)$, i.e., positive values for ${d}_{n}\left(t\right)$. Figure 6 plots ${d}_{n}\left(t\right)$ in all points of the observed ${T}_{1}$ or ${T}_{2}$. The plot shows positive values that increase with $t$ and indicates that the treatment causes longer time to blindness. The effect of treatment also increases with time. The bivariate median inactivity functions are on the right side of Figure 6 to provide a better comparison of these functions.

5.   Conclusions

Assuming an increasing $\alpha$-QIT function, I define a new estimator for this function. It is proven that the proposed estimator is consistent. It is asymptotically close to the usual estimator in the sense that the difference to the usual estimator converges to zero with high probability. It is also shown that the proposed estimator converges weakly to a Gaussian process when normalized. Interestingly, none of the asymptotic results assume that the true $\alpha$-QIT function increases, which increases the applicability of the estimator in general. The simulation results show that the MSE for the proposed increasing estimator is smaller than that of the conventional estimator. When using the proposed estimator, it was found that the laser treatment causes a delay in glare.

Acknowledgements

The author thanks the two anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions.

This work was supported by Researchers Supporting Project number (RSP2024R392), King Saud University, Riyadh, Saudi Arabia.

Data availability statement

The data for this study is available in the "diabetic" dataset in the "survival" package in R.

Conflict of interest

The author declares that there is no conflict of interest.