Bivariate median residual life with BMeRL

3.1 Complete data

Let \((X_{1i}, X_{2i})\), \(i = 1, \ldots, n\), be a random sample. The empirical joint survival function is

\[\bar F_n(x_1, x_2) = \frac{1}{n} \sum_{i=1}^{n} I(X_{1i} > x_1,\ X_{2i} > x_2),\]

and the empirical bivariate median residual life is obtained by replacing \(\bar F\) with \(\bar F_n\) in the definition above (Kayid et al., 2020).

The population quantity \(q_1(x)\) is the median of the conditional distribution of \(X_1 - x_1\) given \(X_1 > x_1\) and \(X_2 > x_2\), because

\[\bar F^{(1)}(t;\ x) = \frac{\bar F(x_1 + t,\ x_2)}{\bar F(x_1, x_2)} = \frac12 \quad \Longleftrightarrow \quad \bar F(x_1 + t,\ x_2) = \frac12\, \bar F(x_1, x_2).\]

The estimator applies the same idea to the data. Let

\[R(x) = \{i :\ X_{1i} > x_1,\ X_{2i} > x_2\}, \qquad N(x) = |R(x)|,\]

be the subjects beyond the point and their number, and let \(Y_i = X_{1i} - x_1\), \(i \in R(x)\), be their remaining times of the first component, ordered as \(Y_{(1)} \le Y_{(2)} \le \cdots \le Y_{(N)}\) with \(N = N(x)\). For \(t \ge 0\), a subject is beyond \((x_1 + t, x_2)\) exactly when it is in \(R(x)\) and \(Y_i > t\), so

\[\bar F_n^{(1)}(t;\ x) = \frac{\bar F_n(x_1 + t,\ x_2)}{\bar F_n(x_1, x_2)} = \frac{1}{N} \sum_{i \in R(x)} I(Y_i > t).\]

That is, the empirical conditional survival function is the ordinary empirical survival function of the remaining times \(Y_i\) of the subjects in \(R(x)\); the other \(n - N\) subjects play no role.

Because \(\bar F_n^{(1)}\) is a step function, the equation \(\bar F_n^{(1)}(t;\ x) = 1/2\) often has an interval of solutions or none at all. The package uses the two crossing points

\[t_{\le} = \inf\left\{t :\ \bar F_n^{(1)}(t;\ x) \le \tfrac12\right\}, \qquad t_{<} = \inf\left\{t :\ \bar F_n^{(1)}(t;\ x) < \tfrac12\right\},\]

and sets \(\hat q_1(x) = (t_{\le} + t_{<})/2\). At \(t = Y_{(k)}\) (ignoring ties) exactly \(N - k\) remaining times exceed \(t\), so

\[\bar F_n^{(1)}(Y_{(k)};\ x) = \frac{N - k}{N} \le \frac12 \iff k \ge \frac N2, \qquad \frac{N - k}{N} < \frac12 \iff k > \frac N2 .\]

Hence \(t_{\le} = Y_{(\lceil N/2 \rceil)}\) and \(t_{<} = Y_{(\lfloor N/2 \rfloor + 1)}\):

• if \(N\) is odd, both equal \(Y_{((N + 1)/2)}\), the middle value;
• if \(N\) is even, they are \(Y_{(N/2)}\) and \(Y_{(N/2 + 1)}\), and their average is the usual sample median.

With ties the same argument holds, and the result is again the sample median. Since subtracting \(x_1\) does not change the order, the estimator is the conventional sample median of the subjects beyond the point:

\[\hat q_1(x) = \mathrm{median}\{X_{1i} :\ X_{1i} > x_1,\ X_{2i} > x_2\} - x_1,\]

\[\hat q_2(x) = \mathrm{median}\{X_{2i} :\ X_{1i} > x_1,\ X_{2i} > x_2\} - x_2.\]

Using only \(t_{\le}\), as the infimum in the definition would suggest, gives the lower median \(Y_{(N/2)}\) when \(N\) is even, which is biased downwards in small samples; the midpoint removes this bias. Note that both components use the same subjects \(R(x)\): the dependence between \(X_1\) and \(X_2\) enters through the choice of this set. If \(R(x)\) is empty, both estimates are set to 0.

3.2 Right-censored data

Suppose each subject has one censoring time \(C_i\), independent of \((X_{1i}, X_{2i})\), that applies to both lifetimes. We observe \(T_{ji} = \min(X_{ji}, C_i)\) and \(\delta_{ji} = I(X_{ji} \le C_i)\) for \(j = 1, 2\). Following Lin and Ying (1993), the joint survival function is estimated by

\[\tilde F_n(x_1, x_2) = \frac{n^{-1} \sum_{i=1}^{n} I(T_{1i} > x_1,\ T_{2i} > x_2)} {\hat G\left(\max(x_1, x_2)\right)},\]

where \(\hat G\) is the Kaplan-Meier estimator of the censoring survival function \(P(C > t)\), computed from the times \(\max(T_{1i}, T_{2i})\) with censoring indicators \(1 - \delta_{1i}\delta_{2i}\). The numerator counts the subjects still under observation beyond the point, and the denominator corrects for those lost to censoring.

As for complete data, the estimator is built from the conditional survival function of the remaining time. Let

\[\tilde N(x_1, x_2) = \sum_{i=1}^{n} I(T_{1i} > x_1,\ T_{2i} > x_2)\]

be the number of subjects whose observed times both lie beyond the point \((x_1, x_2)\). With this notation the Lin and Ying estimator reads

\[\tilde F_n(x_1, x_2) = \frac{\tilde N(x_1, x_2)}{n\, \hat G\left(\max(x_1, x_2)\right)}.\]

The estimated conditional survival function of the remaining time of the first component is the ratio of \(\tilde F_n\) at the shifted point \((x_1 + t, x_2)\) to \(\tilde F_n\) at the point itself. The factor \(1/n\) cancels, and the ratio splits into two factors:

\[\tilde F_n^{(1)}(t;\ x) = \frac{\tilde F_n(x_1 + t,\ x_2)}{\tilde F_n(x_1, x_2)} = \underbrace{\frac{\tilde N(x_1 + t,\ x_2)}{\tilde N(x_1, x_2)}}_{\text{(a) observed share}} \times \underbrace{\frac{\hat G\left(\max(x_1, x_2)\right)}{\hat G\left(\max(x_1 + t,\ x_2)\right)}}_{\text{(b) censoring correction}}, \qquad t \ge 0.\]

Factor (a) is the share of the subjects observed beyond \((x_1, x_2)\) whose first observed time also exceeds \(x_1 + t\). It is exactly the complete-data quantity \(\bar F_n^{(1)}(t;\ x)\) of Section 3.1, computed from the observed times \(T_{1i}\) instead of the lifetimes \(X_{1i}\). Under censoring it is too small: a subject censored between \(x_1\) and \(x_1 + t\) leaves the count although its event has not occurred.

Factor (b) corrects this. Since \(\hat G\) is a non-increasing function and \(\max(x_1 + t,\ x_2) \ge \max(x_1, x_2)\), factor (b) is at least 1. Its inverse is the estimated conditional probability of remaining uncensored from the point to the shifted point, so dividing by it inflates factor (a) by the right amount.

The median residual life is then found exactly as in Section 3.1: it is the value of \(t\) at which \(\tilde F_n^{(1)}(t;\ x)\) falls to one half. Because \(\tilde F_n^{(1)}\) is a step function, we take the two crossing points

\[\tilde t_{\le} = \inf\left\{t :\ \tilde F_n^{(1)}(t;\ x) \le \tfrac12\right\}, \qquad \tilde t_{<} = \inf\left\{t :\ \tilde F_n^{(1)}(t;\ x) < \tfrac12\right\},\]

that is, the first time the curve reaches one half and the first time it falls strictly below one half, and define the estimator as their midpoint,

\[\tilde q_1(x) = \frac{\tilde t_{\le} + \tilde t_{<}}{2}.\]

The second component \(\tilde q_2(x)\) is obtained in the same way from

\[\tilde F_n^{(2)}(t;\ x) = \frac{\tilde F_n(x_1,\ x_2 + t)}{\tilde F_n(x_1, x_2)}.\]

Because \(\hat G\) can drop at a censoring time, \(\tilde F_n^{(1)}\) need not be monotone; taking the first crossing, as in the definition, handles this.

Without censoring, \(\hat G \equiv 1\), factor (b) equals 1, \(\tilde F_n^{(1)} = \bar F_n^{(1)}\), and \(\tilde q_1(x)\) reduces to the sample median of Section 3.1. The estimates then agree exactly with those of BMeRL_Complete().

Under censoring an estimate can be undefined (NA): if \(\tilde F_n^{(1)}\) never falls to \(1/2\) before the largest observed time, the median residual life is longer than the follow-up can show.

3.3 Two groups with an order restriction

Now suppose there are two groups. Let \((X_{1i}, X_{2i})\), \(i = 1, \ldots, n_1\), be a random sample from the group expected to have the smaller median residual life, with joint distribution function \(F\), and let \((Y_{1j}, Y_{2j})\), \(j = 1, \ldots, n_2\), be a random sample from the other group, with joint distribution function \(G\). Let \(q_F(x) = \left(q_{F,1}(x),\ q_{F,2}(x)\right)\) and \(q_G(x) = \left(q_{G,1}(x),\ q_{G,2}(x)\right)\) be their bivariate median residual life functions. We assume

\[q_F(x) \le q_G(x) \qquad \text{for all } x,\]

component by component. The unrestricted estimates \(\hat q_F\) and \(\hat q_G\) of Section 3.1 or 3.2 may violate this order, especially where few subjects remain. Khan and Ghebremichael (2026) combine them through the pooled estimate

\[\hat q_R(x) = w_1\, \hat q_F(x) + w_2\, \hat q_G(x), \qquad w_1 = \frac{n_1}{n_1 + n_2}, \quad w_2 = 1 - w_1,\]

which always lies between the two, and define the restricted estimators

\[q_F^{*}(x) = \min\left(\hat q_R(x),\ \hat q_F(x)\right), \qquad q_G^{*}(x) = \max\left(\hat q_R(x),\ \hat q_G(x)\right),\]

where the minimum and maximum are taken component by component. If the unrestricted estimates already respect the order, they are kept unchanged; otherwise both groups receive the pooled value. The restricted estimators therefore always satisfy \(q_F^{*} \le q_G^{*}\), and they change the unrestricted estimates only where the order is violated.

Estimation with complete data. At a point \(x = (x_1, x_2)\) the estimators are computed in three steps, separately for each component \(j = 1, 2\).

1. Unrestricted estimates in each group. Apply Section 3.1 to each sample on its own. In the first group, take the subjects beyond the point, \[R_F(x) = \{i :\ X_{1i} > x_1,\ X_{2i} > x_2\},\] and compute the sample medians of their remaining times, \[\hat q_{F,1}(x) = \mathrm{median}\{X_{1i} :\ i \in R_F(x)\} - x_1, \qquad \hat q_{F,2}(x) = \mathrm{median}\{X_{2i} :\ i \in R_F(x)\} - x_2.\] In the second group, \(\hat q_{G,1}(x)\) and \(\hat q_{G,2}(x)\) are computed in the same way from \(R_G(x) = \{j :\ Y_{1j} > x_1,\ Y_{2j} > x_2\}\).
2. Pooled estimate. Combine the two groups with the weights \(w_1\) and \(w_2\): \[\hat q_{R,j}(x) = w_1\, \hat q_{F,j}(x) + w_2\, \hat q_{G,j}(x).\] The weights depend on the full sample sizes \(n_1\) and \(n_2\), not on the numbers of subjects beyond the point, so they are the same at every \(x\).
3. Restriction. Set \[q_{F,j}^{*}(x) = \min\left(\hat q_{R,j}(x),\ \hat q_{F,j}(x)\right), \qquad q_{G,j}^{*}(x) = \max\left(\hat q_{R,j}(x),\ \hat q_{G,j}(x)\right).\] If \(\hat q_{F,j}(x) \le \hat q_{G,j}(x)\), then \(\hat q_{F,j} \le \hat q_{R,j} \le \hat q_{G,j}\), and both estimates are kept. If \(\hat q_{F,j}(x) > \hat q_{G,j}(x)\), the order is violated, and both groups receive the pooled value \(\hat q_{R,j}(x)\).

Because the restriction is applied to each component separately, it can change \(q_1\) at a point while leaving \(q_2\) unchanged, or the reverse.

Estimation with censored data. Each group now has its own censoring times. In the first group we observe \((T^{F}_{1i}, T^{F}_{2i}, \delta^{F}_{1i}, \delta^{F}_{2i})\), \(i = 1, \ldots, n_1\), and in the second group \((T^{G}_{1j}, T^{G}_{2j}, \delta^{G}_{1j}, \delta^{G}_{2j})\), \(j = 1, \ldots, n_2\). The three steps are the same, except for step 1.

1. Unrestricted estimates in each group. Apply Section 3.2 to each sample on its own. The Kaplan-Meier estimator of the censoring distribution is computed separately in each group, \(\hat H_F\) from the first sample and \(\hat H_G\) from the second, because the two groups may be followed for different lengths of time. The Lin and Ying estimator of the first group is \[\tilde F_{n_1}(x_1, x_2) = \frac{n_1^{-1} \sum_{i=1}^{n_1} I(T^{F}_{1i} > x_1,\ T^{F}_{2i} > x_2)} {\hat H_F\left(\max(x_1, x_2)\right)},\] and \(\tilde G_{n_2}\) is defined in the same way from the second sample with \(\hat H_G\). The unrestricted estimates \(\tilde q_{F,j}(x)\) and \(\tilde q_{G,j}(x)\) are then obtained from \(\tilde F_{n_1}\) and \(\tilde G_{n_2}\) with the crossing-point rule of Section 3.2.
2. Pooled estimate and 3. Restriction are exactly as for complete data, with \(\tilde q\) in place of \(\hat q\).

If the unrestricted estimate of either group is undefined at a point (NA, see Section 3.2), the pooled estimate is undefined there too, and so are both restricted estimates. Without censoring, \(\hat H_F \equiv 1\) and \(\hat H_G \equiv 1\), and the censored-data estimators reduce exactly to the complete-data estimators.

3.4 Simultaneous confidence intervals

The asymptotic variances of the estimators are hard to estimate, and for the restricted estimators the limiting distribution is not normal where the two groups meet. The package therefore computes simultaneous confidence intervals by the percentile bootstrap method of Mandel and Betensky (2008).