Spearman-like correlation measure adjusting for covariates in bivariate survival data

Abstract

We propose an extension of Spearman's correlation for censored continuous and discrete data that permits covariate adjustment. Previously proposed nonparametric and semiparametric Spearman's correlation estimators require either nonparametric estimation of the bivariate survival surface or parametric assumptions about the dependence structure. In practice, nonparametric estimation of the bivariate survival surface is difficult, and parametric assumptions about the correlation structure may not be satisfied. Therefore, we propose a method that requires neither and uses only the marginal survival distributions. Our method estimates the correlation of probability-scale residuals, which has been shown to equal Spearman's correlation when there is no censoring. Because this method relies only on marginal distributions, it tends to be less variable than the previously suggested nonparametric estimators, and the confidence intervals are easily constructed. Although under censoring, it is biased for Spearman's correlation as our simulations show, it performs well under moderate censoring, with a smaller mean squared error than nonparametric approaches. We also extend it to partial (adjusted), conditional, and partial-conditional correlation, which makes it particularly relevant for practical applications. We apply our method to estimate the correlation between time to viral failure and time to regimen change in a multisite cohort of persons living with HIV in Latin America.

Keywords: Spearman's correlation; bivariate survival data; nonparametric; probability-scale residuals; semiparametric.

Figure 1

Contour plots for the absolute value of ${\rho }_{PSR}$ as a function of the proportion censored when ${\rho }_S=1$. The left and right columns represent scenarios of perfect positive and perfect negative correlations, respectively. The top and bottom rows are strict type I and unbounded bivariate censoring, respectively. Each contour represents a change of 0.016 absolute correlation value.

Figure 2

Type I error rate (top row) and power (middle and bottom rows) for unadjusted correlation and sample size of 200. The columns represent different types and proportions of censoring. The following methods are presented: 1) ${\widehat{\rho }}_{PSR}$ with Stute’s estimating equations is represented by ○; 2) ${\widehat{\rho }}_S^H$ (Eden et al., 2021) is represented by △; 3) ${\widehat{S}}_N$ (Dabrowska, 1986) is represented by +; 4) ${\widehat{T}}_N$ (Dabrowska, 1986) is represented by ×; 5) $U_N$ (Shih and Louis, 1996) is represented by ◇; 6) ${\widehat{V}}_N$ (Shih and Louis, 1996) is represented by ▽.

Figure 3

Point estimate ±SD for unadjusted correlation and sample size of 200. The columns represent different types and proportions of censoring. Solid horizontal lines represent ${\rho }_S$, and dotted horizontal lines represent ${\rho }_{PSR}$. The numbers on the bottom represent the corresponding RMSEs. The following methods are presented: 1) ${\widehat{\rho }}_{PSR}$ with Stute’s estimating equations is represented by ○; 2) ${\widehat{\rho }}_S^H$ (Eden et al., 2021) is represented by △. The horizontal dotted and solid lines are the population parameters ${\rho }_{PSR}$ and ${\rho }_S$, respectively.

Figure 4

Point estimate ±SD for partial correlation and sample size of 200. The columns represent different types and proportions of censoring. The numbers on the bottom represent the corresponding RMSEs. The following methods are presented: 1) Cox with estimating equations based on the full likelihood represented by ○; 2) Exponential survival model represented by △; 3) Log-normal survival model represented by □.

Figure 5

Top row: the population parameters for ${\rho }_S(Z)$ (in black) and ${\rho }_{PSR}(Z)$ (in gray) as functions of $Z$. Middle row: bias of ${\widehat{\rho }}_{PSR}(Z)$ for ${\rho }_S(Z)$ (in black) and ${\widehat{\rho }}_{PSR}(Z)$ for ${\rho }_{PSR}(Z)$ (in gray) as functions of $Z$. Bottom row: coverage probability of ${\widehat{\rho }}_{PSR}(Z)$ for ${\rho }_S(Z)$ (in black) and ${\widehat{\rho }}_{PSR}(Z)$ for ${\rho }_{PSR}(Z)$ (in gray) as functions of $Z$. Frank’s copula was used to model correlation. Survival probabilities were modeled using Cox proportional hazards regression (true model) with variance estimated using full likelihood score equations. For bell-shaped conditional correlation, the linear models were fit with a quadratic term. The data were simulated 1000 times with a sample size of 500. Unbounded univariate censoring of (0.3,0.3) was applied.

Figure 6

Partial-conditional correlation of PSRs between time to viral failure and time to regimen change computed as a function of CD4 count and age at the time of first ART. The left column shows results for children, the right column for adults. Cox proportional hazards regression was used to model survival probabilities. The variance was estimated using full likelihood score equations. Note that because its skewed nature, the $x$-axis of the CD4 counts is on the square-root scale.