Nonparametric estimation of Spearman's rank correlation with bivariate survival data

Abstract

We study rank-based approaches to estimate the correlation between two right-censored variables. With end-of-study censoring, it is often impossible to nonparametrically identify the complete bivariate survival distribution, and therefore it is impossible to nonparametrically compute Spearman's rank correlation. As a solution, we propose two measures that can be nonparametrically estimated. The first measure is Spearman's correlation in a restricted region. The second measure is Spearman's correlation for an altered but estimable joint distribution. We describe population parameters for these measures and illustrate how they are similar to and different from the overall Spearman's correlation. We propose consistent estimators of these measures and study their performance through simulations. We illustrate our methods with a study assessing the correlation between the time to viral failure and the time to regimen change among persons living with HIV in Latin America who start antiretroviral therapy.

Keywords: HIV; Spearman's correlation; bivariate survival; nonparametric; viral failure.

FIGURE 1

Illustration of bivariate distributions underlying the three population parameters. Left: Original distribution over [0, ∞) × [0, ∞), which has Spearman’s correlation ρ_S. Middle: Conditional distribution over Ω_R = Ω, which has Spearman’s correlation ${\rho }_{S\mid {\Omega }_R}$. Right: Mixture-like distribution S^H over region Ω ∪ [0, τ_X] × τ_Y ∪ τ_X × [0, τ_Y], which has Spearman’s correlation ${\rho }_S^H$

FIGURE 2

Restricted Spearman’s correlation, ${\rho }_{S\mid {\Omega }_R}$ (left panel), and highest rank Spearman’s correlation, ${\rho }_S^H$ (right panel), for Frank’s copula family for different restricted regions defined by τ_X and τ_Y (τ_X = τ_Y): 0.5th (50% censored), 0.6th (40% censored), 0.8th (20% censored) quantiles. A diagonal gray line is added for reference. Although the plots are generated based on data under strict type I censoring, the population parameters are the same for generalized type I censoring and are invariant to the rate of censoring within the restricted region

FIGURE 3

Example of bivariate distributions and their population parameters ρ_S with no censoring, and ${\rho }_S^H$ and ${\rho }_{S\mid {\Omega }_R}$ with strict type I censoring with Ω_R = [0, τ_X) × [0, τ_Y). The proportions of observed double events in the left, middle, and right panels are 25%, 43%, and 7%, respectively. Drawn are 1000 points randomly selected from the underlying distributions. Although the plots are based on strict type I censoring, the population parameters are the same for generalized type I censoring and are invariant to the rate of censoring within Ω_R

FIGURE 4

Point estimates (X-axis) vs population parameters (Y-axis) under different bivariate censoring scenarios. The top and second rows are ${\widehat{\rho }}_S^H$ and ${\widehat{\rho }}_{IMI}$ as estimators of the overall Spearman’s correlation, ρ_S. The third row is ${\widehat{\rho }}_S^H$ as an estimator of ${\rho }_S^H$. The bottom row is ${\widehat{\rho }}_{S\mid {\Omega }_R}$ as an estimator of ${\rho }_{S\mid {\Omega }_R}$. The columns represent Clayton’s and Frank’s copulas. The population parameters for Clayton’s family are 0, 0.2, and 0.6 for all estimates. For Frank’s family, the population parameters of ρ_S are −0.6, −0.2, 0.2, and 0.6; the population parameters of ${\rho }_S^H$ are −0.512, −0.173, 0.180, and 0.545; the population parameters of ${\rho }_{S\mid {\Omega }_R}$ are −0.098, −0.042, 0.058, and 0.261. The dots are the mean point estimates based on 1000 simulations. The shaded areas represent the 0.025th and 0.975th quantiles. For generalized type I censoring, the restricted region, Ω_R, was defined by the median survival times

FIGURE 5

Upper left: Kaplan–Meier curves for time to viral failure and time to regimen change, where time is measured in years. Upper right: bivariate probability mass function for the mixture-like distribution, ${\widehat{S}}^H(dx,dy)$. Lower left: conditional bivariate probability mass function for 15-year follow-up. Lower right: conditional bivariate probability mass function for 10-year follow-up. For probability mass functions, the bars on the left and on the bottom represent histograms of the univariate survival mass for each event. The probability mass function was computed from the Dabrowska’s survival surface and then aggregated over half-year bivariate time periods. After aggregation, any negative values were set to 0. Lighter shade represents smaller values