--- title: "Bivariate median residual life with BMeRL" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Bivariate median residual life with BMeRL} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 3.6) ``` ```{r} library(BMeRL) ``` ## 1. Introduction Survival and hazard functions describe the risk of an event from the time of diagnosis. In many clinical settings, however, patients and physicians ask a different question: given that the event has not yet happened, how much longer is it likely to take? This question is answered by the *residual life* function. When a new treatment is started part of the way through follow-up, the residual life is often more informative than the survival function, because it describes the remaining time from the current state of the patient. The mean residual life is the most studied measure of residual life, but it has practical drawbacks. It may not exist for heavy-tailed distributions, it is sensitive to extreme observations. The *median residual life* (MeRL) avoids these problems: it always exists, it is robust to outliers and skewness. Many studies record two dependent event times on the same subject. In HIV research, for example, treatment success is judged by both virological suppression and immunological restoration, and the two times are corelated. Analysing each time separately ignores this dependence. The *bivariate* median residual life describes both remaining times jointly, given that neither event has happened by a chosen point. Often two groups are also expected to be ordered. For example, patients with a coexisting disease are expected to have shorter residual lives than patients without it (or longer recovery times, depending on the event). Khan and Ghebremichael (2026) proposed estimators of the bivariate median residual life of two groups that respect such an order, and showed that they have smaller mean squared error than the unrestricted estimators. The **BMeRL** package implements these methods: | Estimator | Surface plot | Data | Groups | |---|---|---|---| | `BMeRL_Complete()` | `BMeRL_Surface_complete()` | complete | one | | `BMeRL_Censored()` | `BMeRL_Surface_censored()` | right-censored | one | | `BMeRL_Complete_ordered()` | `BMeRL_Surface_complete_ordered()` | complete | two, ordered | | `BMeRL_Censored_ordered()` | `BMeRL_Surface_censored_ordered()` | right-censored | two, ordered | Every estimator can also return bootstrap simultaneous confidence intervals (SCI), computed by the method of Mandel and Betensky (2008). ## 2. The bivariate median residual life Let $X = (X_1, X_2)$ be a non-negative random vector with joint survival function $$\bar F(x_1, x_2) = P(X_1 > x_1,\ X_2 > x_2).$$ For a point $x = (x_1, x_2)$ at which no event has occurred, the conditional survival function of the remaining time of component $i$ is $$\bar F^{(i)}(t_i;\ x) = P(X_i - x_i > t_i \mid X_1 > x_1,\ X_2 > x_2), \qquad i = 1, 2.$$ The bivariate median residual life function at $x$ is the vector $q(x) = \left(q_1(x),\ q_2(x)\right)$, where $$q_i(x) = \inf\left\{t_i :\ \bar F^{(i)}(t_i;\ x) = \tfrac12\right\}, \qquad i = 1, 2.$$ That is, $q_i(x)$ is the median of the remaining time $X_i - x_i$, given that neither event has occurred by $x$. Written in terms of the joint survival function, $$q_1(x) = \inf\left\{t_1 :\ \bar F(x_1 + t_1,\ x_2) = \tfrac12\,\bar F(x_1, x_2)\right\},$$ $$q_2(x) = \inf\left\{t_2 :\ \bar F(x_1,\ x_2 + t_2) = \tfrac12\,\bar F(x_1, x_2)\right\}.$$ So $q_1(x)$ is the additional time, measured from $x_1$ along the first axis with $x_2$ held fixed, at which the joint survival falls to half its value at $x$; the second component is defined in the same way along the second axis. For example, if $X_1$ and $X_2$ follow the bivariate Gumbel distribution with $\bar F(x_1, x_2) = \exp(-x_1 - x_2 - \theta x_1 x_2)$, then $$q_1(x) = \frac{\log 2}{1 + \theta x_2}, \qquad q_2(x) = \frac{\log 2}{1 + \theta x_1}.$$ The remaining time of one component depends on how long the *other* component has survived. This is exactly the information that two separate univariate analyses cannot show. ## 3. Estimation ### 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). ## 4. Example data The package includes two **artificial** data sets. They were simulated for this package to illustrate its use and are not real data. Two data sets are enough to illustrate all four estimators: * `complete` holds the lifetimes (in months) of 120 subjects, without censoring. It is used by `BMeRL_Complete()` and `BMeRL_Complete_ordered()`. * `censored` holds the same subjects observed under right censoring, with event indicators `status1` and `status2` (1 = event observed, 0 = censored). It is used by `BMeRL_Censored()` and `BMeRL_Censored_ordered()`. The lifetimes were drawn from the bivariate exponential distribution of Gumbel, with joint survival function $\bar F(z_1, z_2) = \exp(-z_1 - z_2 - \theta z_1 z_2)$ and dependence parameter $\theta = 0.5$, scaled to a mean of 12 months in group `"Yes"` and 18 months in group `"No"`. For the censored data, each subject was given one censoring time $C$ from the uniform distribution on 0 to 60 months, applied to both lifetimes: the observed times are the minima of the lifetimes and $C$, and a status is 0 when the lifetime exceeds $C$. About 42% of the subjects (38% in group `"Yes"` and 46% in group `"No"`) have at least one censored time. All times are rounded to 0.1 month. The column `disease` splits the subjects into two groups, `"Yes"` (50 subjects) and `"No"` (70 subjects), with shorter lifetimes in group `"Yes"`. Only the ordered functions use this column; the other functions ignore it and treat all 120 subjects as one sample. ```{r} head(complete, 3) head(censored, 3) table(complete$disease) ``` In all functions, the column names can be given unquoted, as in the examples below, or as character strings. ## 5. Complete data `BMeRL_Complete()` estimates the bivariate median residual life at one or more points $(x_1, x_2)$. The two vectors `x1` and `x2` are paired, so the call below evaluates the estimator at $(0, 0)$, $(5, 5)$ and $(10, 5)$. ```{r} BMeRL_Complete(x1 = c(0, 5, 10), x2 = c(0, 5, 5), dataset = complete, time1 = time1, time2 = time2, SCI = FALSE) ``` The second row is read as follows: among subjects with neither event by 5 months, half have their first event more than `q1` months later, and half have their second event more than `q2` months later. **In practice.** Section 3.1 showed that the estimate is the sample median of the remaining times of the subjects beyond the point. This can be checked directly at $(5, 5)$. First find the set $R(x)$ and its size $N$: ```{r} beyond <- subset(complete, time1 > 5 & time2 > 5) N <- nrow(beyond) N ``` The two crossing points of $\bar F_n^{(1)}$ are the order statistics $Y_{(\lceil N/2 \rceil)}$ and $Y_{(\lfloor N/2 \rfloor + 1)}$ of the remaining times, and their midpoint is $\hat q_1$: ```{r} Y <- sort(beyond$time1 - 5) t_le <- Y[ceiling(N / 2)] t_lt <- Y[floor(N / 2) + 1] c(t_le = t_le, t_lt = t_lt, q1 = (t_le + t_lt) / 2) ``` This is the conventional sample median, for both components: ```{r} c(q1 = median(beyond$time1) - 5, q2 = median(beyond$time2) - 5) ``` Both agree with the second row of the table above. With `SCI = TRUE` (the default for this function), the simultaneous confidence intervals of Section 3.4 are added. Use `seed` for reproducible results and `verbose = FALSE` to suppress the progress messages. The examples in this vignette use `B = 200` to keep the running time short; the default `B = 1000` is recommended in practice. ```{r} out <- BMeRL_Complete(c(0, 5), c(0, 5), complete, time1, time2, B = 200, seed = 1, verbose = FALSE) out ``` ## 6. Right-censored data `BMeRL_Censored()` has the same form, with two extra arguments for the event indicators. By default it returns point estimates only. ```{r} BMeRL_Censored(c(0, 5, 10), c(0, 5, 5), censored, time1, time2, status1, status2) ``` **In practice.** Section 3.2 showed that without censoring the censoring correction equals 1 and the estimator reduces to the sample median. To check this, treat every time in `complete` as an observed event and apply `BMeRL_Censored()`; the result is identical to that of `BMeRL_Complete()`: ```{r} all_events <- transform(complete, status1 = 1, status2 = 1) fit_cens <- BMeRL_Censored(c(0, 5, 10), c(0, 5, 5), all_events, time1, time2, status1, status2) fit_comp <- BMeRL_Complete(c(0, 5, 10), c(0, 5, 5), complete, time1, time2, SCI = FALSE) all.equal(fit_cens[, c("q1", "q2")], fit_comp[, c("q1", "q2")]) ``` With censoring the estimates differ from the naive sample median of the observed times beyond the point, which ignores that some of those subjects were censored before their events: ```{r} beyond_c <- subset(censored, time1 > 5 & time2 > 5) c(naive_q1 = median(beyond_c$time1) - 5, naive_q2 = median(beyond_c$time2) - 5) ``` The naive values are too small, because censored times are shorter than the true event times; the correction factor of Section 3.2 accounts for this. To add simultaneous confidence intervals, set `SCI = TRUE`: ```{r} BMeRL_Censored(c(0, 5), c(0, 5), censored, time1, time2, status1, status2, SCI = TRUE, B = 200, seed = 1, verbose = FALSE) ``` ## 7. Two groups with an order restriction The ordered functions need two more arguments: `group`, the column that defines the two groups, and `smaller`, the value of that column for the group expected to have the smaller median residual life. Here the subjects with `disease = "Yes"` are expected to have shorter residual lives. ```{r} est <- BMeRL_Complete_ordered(c(0, 5, 20), c(0, 5, 1), complete, time1, time2, disease, smaller = "Yes", SCI = FALSE) est ``` The unrestricted estimates of Section 3.1 are stored as an attribute: ```{r} attr(est, "unrestricted") ``` Comparing the two tables shows where the restriction changed anything: where the unrestricted estimates already satisfy the order they are kept, and where they do not, both groups receive the pooled estimate. For censored data the call is the same, with the event indicators added. With `SCI = TRUE` (the default for the ordered functions), the intervals are simultaneous over $q_1$ and $q_2$ of both groups at all points: ```{r} BMeRL_Censored_ordered(c(0, 5), c(0, 5), censored, time1, time2, status1, status2, disease, smaller = "Yes", B = 200, seed = 1, verbose = FALSE) ``` Where the two groups are close, the restricted estimators are not smooth functions of the data, and the percentile bootstrap may cover less than the nominal level there. A simulation study is advisable before relying on the intervals in such regions. ## 8. Surface plots A table of estimates shows the median residual life at a few points; a surface plot shows the whole function. The surface functions evaluate an estimator on a grid and draw $\hat q_1$ and $\hat q_2$ side by side. Supplying `x1.grid` and `x2.grid` focuses the plot on the region with enough data; grid points without an estimate are left blank. ```{r} BMeRL_Surface_complete(complete, time1, time2, x1.grid = 0:30, x2.grid = 0:30) ``` ```{r} BMeRL_Surface_censored(censored, time1, time2, status1, status2, x1.grid = 0:30, x2.grid = 0:30) ``` For two groups there is one row per group, the smaller group first. Both groups share one height scale and one colour scale per component, so the surfaces can be compared directly, and the restricted surface of the smaller group never lies above that of the larger group. ```{r, fig.height = 7} BMeRL_Surface_complete_ordered(complete, time1, time2, disease, "Yes", x1.grid = 0:30, x2.grid = 0:30) ``` ```{r, fig.height = 7} BMeRL_Surface_censored_ordered(censored, time1, time2, status1, status2, disease, "Yes", x1.grid = 0:30, x2.grid = 0:30) ``` The view can be changed with `theta` and `phi`, and further arguments are passed to `persp()`, for example `ticktype = "detailed"` for numbered axes. With `plot = FALSE` the functions only return the grid of estimates. ## References Khan, R. A. and Ghebremichael, M. (2026). Bivariate quantile residual life functions for estimating immunological and virological recovery in HIV patients. *Biometrical Journal*. doi:10.1002/bimj.70112 Kayid, M., Shafaei Noughabi, M. and Abouammoh, A. M. (2020). A nonparametric estimator of bivariate quantile residual life model with application to tumor recurrence data set. *Journal of Classification*. doi:10.1007/s00357-018-9300-z Mandel, M. and Betensky, R. A. (2008). Simultaneous confidence intervals based on the percentile bootstrap approach. *Computational Statistics and Data Analysis*, 52(4), 2158-2165. doi:10.1016/j.csda.2007.07.005 Lin, D. Y. and Ying, Z. (1993). A simple nonparametric estimator of the bivariate survival function under univariate censoring. *Biometrika*, 80(3), 573-581. doi:10.1093/biomet/80.3.573