A Quantitative Concordance Measure for Comparing and Combining Treatment Selection Markers

2.1 A quantitative concordance measure

Consider a clinical trial in which each patient receives a randomized treatment T (1 experimental; −1 control). We assume for simplicity that P(T=1)=P(T=−1)=1/2, although the main ideas extend easily to different allocation probabilities. Let Y(t) denote the potential outcome that will realize if a patient receives treatment t∈{1,−1} [22]. Assuming consistency or stable unit treatment value, the actual outcome is just Y=Y(T). Without loss of generality, we assume that higher values of Y are desirable. Let X be a vector of baseline covariates including one or more markers that may be considered in choosing between the two treatments.

We start by considering how to assess the utility of one marker, say V, for treatment selection. To this end, it is important to consider the conditional effect of t=1 versus t=−1 given V, defined as

For a useless marker, the conditional effect function is constant: δ(V)≡E{Y(1)−Y(−1)}=:Δ. A useful marker should have a large amount of variation in the function δ. We assume throughout that higher marker values are intended to predict larger treatment effects, which implies that only “positive" variation in δ is desirable. This motivates our definition of γ=E(G12) with

where V1 and V2 are the marker values of two independent patients, and sgn is the sign function: sgn(u)=1,0,−1 according as u>,=,<0.

We call γ a quantitative concordance measure because G12 involves the actual difference δ(V1)−δ(V2) and not just the sign of the difference. An analogous qualitative concordance measure might be defined as γ∗=E(G12∗) with

It is easy to see that γ∗=1 if V is continuous and δ is strictly increasing. Among such markers, it is certainly desirable to identify those with large variation in δ while considering the distribution of V. That extra information is captured by the quantitative concordance measure γ, which is therefore the focus of this article.

To understand what γ means, suppose the two patients are to receive treatments 1 and −1. If V1\gtV2, then

where the last step makes use of the independence of the two patients. If V1\ltV2, the same argument shows that

If V1=V2, then G12=0. Thus, G12 measures the impact of giving treatment 1 to the patient with a larger V-value (if there is one) versus the opposite treatment allocation, and γ represents the population average impact of treatment allocation based on V versus −V. A large value of γ is clearly desirable. Note also that γ is invariant under a strictly increasing transformation of the marker.

Further insights into γ can be gained by considering the selection impact curve of Song and Pepe [8]. Assume for the moment that V is continuously distributed on an interval. Without loss of generality, we further assume that V is uniformly distributed on (0,1). The selection impact curve is the function

which gives the mean outcome for a treatment regime that assigns treatment 1 to patients with V\gtv and treatment −1 to everyone else. Without specifying a cutoff v, a natural summary of the selection impact curve is its integral over the unit interval, which is analogous to the AUC for an ROC curve. In Appendix A, we show that

Thus, for a given pair of treatments, γ is directly related to the area under the selection impact curve.

2.2 Estimating γ for a given marker

We now consider estimation of γ for a given marker using a random sample of (X,T,Y) (recall that X includes V as a component), with individual subjects denoted by subscripts (e.g., i=1,…,n). Estimation of γ is complicated by the fact that the conditional effect function δ is unknown. However, we note that

which suggests that δ(V) may be estimated by 2TY [7, 19]. In fact, δ(V) could be estimated by 2T(Y−A) for any A=a(X), because E(TA|V)=E{E(TA|X)|V}=0 due to randomization. The inclusion of A is motivated by a potential gain in efficiency, as will become clear later.

The foregoing discussion, together with the definition of γ, suggests the following estimator:

where

Noting that γˆ(a) is a one-sample U-statistic of order 2, it follows from standard asymptotic theory e.g., van der Vaart [23] Chapter 12 that n{γˆ(a)−γ} converges to a normal distribution with mean zero and variance

A consistent estimator of σ2(a) is given by

In Appendix B, we show that the asymptotic variance σ2(a) is minimized by taking a(X) to be

To take advantage of this result, we can specify and estimate a working model E(Y|X)=m(X;α), where m is a known function and α is an unknown finite-dimensional parameter. For example, m(X;α) may be specified as a generalized linear model, and α may be estimated by (iteratively reweighted) least squares. Let αˆ denote the resulting estimate of α; then we can estimate aopt(X) by aˆopt(X)=m(X;αˆ). We show in Appendix B that, whether the model m(X;α) is correct or not, γˆ(aˆopt) is consistent for γ and asymptotically normal, and its asymptotic variance is given by σ2(m(⋅;α∞)), where α∞ is the probability limit of αˆ. If m(X;α) is correctly specified, then α∞ is the true value of α and γˆ(aˆopt) is asymptotically equivalent to γˆ(aopt). Regardless of the (in)correctness of m(X;α), the asymptotic variance σ2(m(⋅;α∞)) is consistently estimated by σˆ2(aˆopt).

2.3 Comparing two markers

We now consider the problem of comparing two markers, say Vb and Vc, both of which are included in X. Write γb=E(Gb12) with

where δb(v)=E{Y(1)−Y(−1)|Vb=v} and Vb1 and Vb2 are the values of Vb for two independent patients. For a given augmentation A=a(X), γb can be estimated by

where

Let γc and γˆc(a) be defined analogously for Vc. Our comparison of the two markers is based on the difference γb−γc, which is consistently estimated by γˆb(a)−γˆc(a). For any fixed a, it follows from standard asymptotic theory that n{γˆb(a)−γˆc(a)−(γb−γc)} converges to a normal distribution with mean zero and variance

A consistent estimator of σbc2(a) is given by

In Appendix C, we show that σbc2(a) is minimized by the same aopt defined by eq. (4). As in Section 2.2, we can estimate aopt(X) by aˆopt(X)=m(X;αˆ), where m(X;α) is a working model for E(Y|X) and αˆ is an estimate of α. The asymptotic behavior of γˆb(aˆopt)−γˆc(aˆopt) is similar to that of γˆ(aˆopt) and can be described by the last paragraph of Section 2.2 with some obvious modifications.

2.4 Combining multiple markers

Suppose V is a sub-vector of X consisting of two or more treatment selection markers. A natural question is how to combine the markers in V into a single marker S=s(V). We propose to combine V by maximizing γs, the γ-value of S. To minimize new notation, we will continue to use the notation in Section 2.1 with V replaced by V. For instance, δ(V) is the conditional effect of t=1 versus t=−1 given V. It is easy to see that γs is maximized by any s that preserves the ordering by δ:

Thus, δ(V) itself would be an ideal choice of S, as would any strictly increasing function of δ(V).

The foregoing discussion suggests that we could set S equal to an estimate of δ(V), which may be obtained under a model for E(Y|T,V) [24]. Noting that

where μ(V)=E{Y(1)+Y(−1)|V}/2, we might parameterize μ(V) and δ(V) separately and estimate the model parameters by (iteratively reweighted) least squares. The optimality of this approach generally depends on the correctness of modeling assumptions, although there may be some robustness with respect to μ(V) [7, 19]. If δ(V) is misspecified, this approach may lead to a sub-optimal choice of S whose interpretation is unclear.

In the rest of this subsection, we consider an alternative approach based on direct maximization of γs, similar in spirit to the approach of Zhang et al. [11] for estimating optimal treatment regimes. For practicality, we restrict attention to a parametric family of functions, s(V;β), indexed by a finite-dimensional parameter β. The specification of s(V;β) may be informed by a parametric model for δ(V), and may also incorporate practical considerations such as simplicity and interpretability. The parameter space for β may need to be restricted so as to ensure that s(V;β1) and s(V;β2) are not strictly increasing functions of each other if β1/=β2. For example, if s(V;β)=β′V, we might restrict β to the unit sphere or, alternatively, fix the first element of β at 1. Once such a parametric family is defined, the problem is then to estimate βopt=argmaxβγ(β), where we write γ(β)=γs(⋅;β) for convenience. The estimand βopt can be interpreted as the true value of β in the following single index model:

where ξ is an unknown, strictly increasing function. Because ξ is unspecified, this approach is more flexible than the one based on a parametric model for δ(V). Even if the single index model is misspecfied, s(⋅;βopt) remains interpretable as a local maximizer of γs (within the specified parametric family).

We propose to estimate βopt by maximizing a consistent estimate of γ(β), which may be obtained as in Section 2.2. Let γˆ(β,a) be defined by eqs (2) and (3) with V replaced by s(V;β). The objective function β↦γˆ(β,a) can be maximized using a grid search algorithm when the effective dimension of β is one or two. For a higher-dimensional β, a grid search may be impractical, and the discreteness of the sign function in eq. (3) may create difficulties in applying a more efficient algorithm. To deal with this problem, we follow the approach of Ma and Huang [18] and apply a gradient descent algorithm to a smoothed version of γˆ(β,a). Specifically, we replace the sign function with λn=2K(⋅/ϵn)−1, where K is a smooth, symmetric distribution function and (ϵn) is a sequence of positive numbers decreasing to 0. This leads to the smoothed estimator

where

It can be shown as in Ma and Huang [18] online supplement that, under regularity conditions, γ˜(β,a) is asymptotically equivalent to γˆ(β,a); see Appendix D for details. Thus, the optimal choice of a and the related discussion in Section 2.2 remain applicable to γ˜(β,a). Empirical rules for choosing ϵn are given in Ma and Huang [18] Section 3.5.

Because of the asymptotic equivalence of γˆ(β,a) and γ˜(β,a), we will for simplicity focus on the unsmoothed version in the rest of this subsection. Let βˆopt be a maximizer of γˆ(⋅,aˆ), where aˆ may be fixed or estimated. The arguments of Han [25] and Sherman [26] can be adapted to show that, under regularity conditions, βˆopt is consistent for βopt and (for a suitable parameterization) asymptotically normal. Of particular interest to us is the quantity γ(βˆopt), which measures the performance of the estimated optimal combination. Asymptotic inference on this quantity can be based on the fact that

where a∞ is the probability limit of aˆ (a sketch proof is given in Appendix D).

In finite samples, a re-substitution bias can arise from the fact that βˆopt is chosen to maximize γˆ(⋅,aˆ) [1]. In general, the maximum of an estimated objective function is an over-estimator of the true objective function evaluated at the maximizer [20, 21]. The re-substitution bias can be removed using the following cross-validation procedure. We propose to partition the study cohort randomly into a specified number, say K, of subsamples that are roughly equal in size. For each k∈{1,…,K}, we use the kth subsample as a validation sample and combine the other subsamples into a training sample. From the training sample we obtain

using the exact same method for obtaining βˆopt. The superscript (−k) in the preceding display emphasizes that each estimate is based on the training sample alone. Next, we evaluate the estimated combination s(⋅;βˆopt(−k)) in the validation sample and obtain γˆ(k)(βˆopt(−k),aˆ(k)), where the superscript (k) indicates that the corresponding estimate is based on the validation sample alone. Because βˆopt(−k) and γˆ(k)(⋅,aˆ(k)) are based on independent data, γˆ(k)(βˆopt(−k),aˆ(k)) is free of re-substitution bias. The final cross-validated estimator of γ(βˆopt) is given by

Note that γˆ(k)(βˆopt(−k),aˆ(k)) is really an estimator of γ(βˆopt(−k)), which measures the utility of s(⋅;βˆopt(−k)) based on a training sample. Because s(⋅;βˆopt) is based on a larger sample and is expected to perform better, γˆ(k)(βˆopt(−k),aˆ(k)) and hence γˆcv(βˆopt) may be biased downward for estimating γ(βˆopt), especially for small K. The choice of K thus represents a trade-off between bias and computational demand – a larger K decreases the downward bias in γˆ(k)(βˆopt(−k),aˆ(k)) at the expense of more computation. Once K is chosen, the size of a subsample is determined as roughly n/K. If this number is large enough for reliable estimation of aopt, we can set aˆ(k) to be an estimate of aopt based on the kth subsample; otherwise we can substitute a fixed a for aˆ(k).

Because the terms γˆ(k)(βˆopt(−k),aˆ(k)), k=1,…,K, are not independent of each other, it is not straightforward to derive the variance of the cross-validated estimator γˆcv(βˆopt). Nonetheless, inference on γ(βˆopt) can be based on nonparametric bootstrap standard errors and confidence intervals. A nonparametric bootstrap sample of the same size as the original sample can be created by sampling with replacement. A cross-validated estimate of γ(βˆopt) can be obtained from each bootstrap sample, repeating all steps in optimization and cross-validation. This procedure can be repeated many times to produce a collection of bootstrap estimates, from which a bootstrap standard error can be obtained as a sample standard deviation and confidence limits as empirical quantiles.

3.1 Analysis of HIV trial data

The ECHO and THRIVE studies enrolled a combined total of 1,374 treatment-naive HIV-infected adults (694 ECHO; 680 THRIVE) in multiple countries. Main inclusion criteria included a screening viral load of at least 5,000 copies/ml and confirmed viral sensitivity to background N[t]RTIs. The study subjects were randomized 1:1 to rilpivirine 25 mg once daily or efavirenz 600 mg once daily, both in conjunction with a background N[t]RTI regimen which consisted of tenofovir disoproxil fumarate (TDF) and emtricitabine (FTC) in ECHO and was based on the investigator’s choice of TDF/FTC, zidovudine/lamivudine (3TC), or abacavir/3TC in THRIVE. Randomization was stratified by screening viral load (≤105, 105–5×105, >5×105 copies/ml) and background N[t]RTI regimen (in THRIVE). A double-dummy design was used to keep patients masked from the assigned treatment. Patients were followed for up to 96 weeks, and the primary endpoint was virologic response at week 48. The observed response rates were 82.2% for rilpivirine and 79.8% for efavirenz.

Here we use the proposed methodology to evaluate, compare and combine the two markers (BVL and BCD4) for choosing between rilpivirine (T=1) and efavirenz (T=−1). For ease of handling, we work with standardized versions of −log(BVL) (minus because higher values of BVL indicate smaller benefits of rilpivirine versus efavirenz) and log(BCD4+1) (where one is added to handle zero counts). The standardization consists of subtracting by the mean and dividing by the standard deviation. We note that γ is invariant under a strictly increasing transformation. The γ-values of the two markers are estimated with A=aˆ(X) equal to 0 (i.e., no augmentation) or an estimate of E(Y), E(Y|V) or E(Y|X), where V consists of the two transformed markers and X consists of age, gender, race, body mass index, and the two markers. We estimate E(Y) with the sample mean of Y, and estimate E(Y|V) and E(Y|X) with logistic regression models that include the specified covariates as linear terms (in addition to an intercept). For brevity, we denote these augmentation terms by Eˆ(Y), Eˆ(Y|V) and Eˆ(Y|X), respectively. Our analysis is restricted to a subset of 1,314 subjects with complete covariate and outcome data.

Table 1 shows the point estimates and standard errors for the γ-values of the two markers as well as their difference. As one might expect, the standard errors do not increase when new information is added to the augmentation term. The reductions are quite large when A changes from 0 to Eˆ(Y), and negligible when A changes from Eˆ(Y|V) to Eˆ(Y|X), suggesting that the additional covariates in X (relative to V) are of limited prognostic value. The point estimates based on A=0 are visibly different from those based on non-zero augmentation terms, possibly due to the large variability for A=0. The non-zero augmentation terms lead to very similar results, which together indicate that BCD4 has a higher γ-value than BVL, although the difference is not statistically significant.

Next, we consider how to combine the two markers linearly, so s(V;β)=β′V with ∥β∥=1. Linear combinations are easy to use in practice, even though they may not attain the global optimum. Our approach is to directly maximize an estimate of γ, which may be obtained using one of the four methods described earlier (with different augmentation terms). The maximization is done using a grid search algorithm based on the representation β=(cosϕ,sinϕ)′, ϕ∈[0,2π). As a comparator, we also consider an alternative approach based on the logistic regression model

where η=(η1,ηX′,ηT,ηTV′)′ is a vector of unknown regression parameters. Fitting this logistic regression model produces a maximum likelihood estimate of ηTV, which is then normalized to provide a candidate combination of the two markers. The γ-value of the candidate combination based on logistic regression is estimated with A=Eˆ(Y|X). A 5-fold cross-validation procedure is also used to estimate the γ-value for each candidate combination (from either approach), with A=Eˆ(Y|X) in the validation step. Bootstrap standard errors are obtained from 1,000 bootstrap samples.

Table 2 shows the results for marker combination. For the proposed approach, the cross-validated γ-estimate tends to increase when new information is added to the augmentation term. Except in the no-augmentation case, the proposed approach produces higher cross-validated γ-estimates than the logistic regression approach, and both approaches suggest that BCD4 should play a much larger role than BVL in the linear combination. In the no-augmentation case, cross-validation reduces the estimate of γ dramatically for the proposed approach.

Example-based simulation

We now evaluate the performance of the proposed methodology in a simulation study mimicking the ECHO and THRIVE trials. We work with the same subset of 1,314 subjects analyzed in Section 3.1 with their original marker values, and generate T and Y randomly according to the observed pattern in the data. Specifically, we let P(T=1|V)=P(T=−1|V)=1/2 and, for each t∈{1,−1}, set P(Y=1|T=t,V)=P(Y(t)=1|V) equal to a nonparametric estimate based on tensor product splines. Because of the curse of dimensionality, this data generation mechanism does not involve the covariates in X that are not in V, which therefore will not be included in the analysis. One thousand datasets are simulated independently and analyzed as in Section 3.1 (except that X=V in this simulation study).

Table 3 shows the results for marker comparison: empirical bias, standard deviation, mean standard error and empirical coverage probability (of 95% Wald confidence intervals) for estimating the γ-values of the two markers as well as their difference with different augmentation terms (0, Eˆ(Y) and Eˆ(Y|V) but not Eˆ(Y|X)). It is clear in Table 3 that the estimators are virtually unbiased and generally become more efficient with increasing information in the augmentation, although the improvement from Eˆ(Y) to Eˆ(Y|V) is quite modest (presumably due to low prognostic value of V). In all cases, the standard errors approximate the empirical standard deviations well on average, and the resulting confidence intervals have close-to-nominal coverage.

The results for marker combination are presented in Table 4, where the proposed approach with the aforementioned augmentation terms is compared with a logistic regression approach based on the following model:

where η=(η1,ηV′,ηT,ηTV′)′. The proposed approach tends to perform better with increasing information in the augmentation, producing candidate combinations that are closer on average to the optimal combination with a higher mean of γ(βˆopt). By the same measures, logistic regression outperforms the proposed approach with no augmentation but not with non-zero augmentations. For the proposed approach, γˆ(βˆopt,aˆ) (without cross-validation) tends to over-estimate γ(βˆopt), especially in the no-augmentation case, and the re-substitution bias is effectively corrected by cross-validation. Although γˆcv(βˆopt) tends to under-estimate γ(βˆopt), it is on average closer to γ(βˆopt) than γˆ(βˆopt,aˆ) is. It is also worth noting that, for the proposed approach, the variability in βˆopt and γ(βˆopt) tends to decrease with increasing information in the augmentation.

3.3 Additional simulations

Here we report additional simulation experiments with a continuous outcome Y, two markers in V, and a separate baseline variable W such that X=(V′,W)′. The three components of X are independently distributed as standard normal. As before, treatment is randomized 1:1, so that P(T=1|X)=P(T=−1|X)=1/2. Given (X,T), we generate Y from the following single index model:

where α1=0, αX=(0.5,−0.5,1)′, βV=(1,1)′ or (1,2)′, ξ may be identity, expit (i.e., ξ(u)=1/{1+exp(−u)}) or sgn, and ϵ∼N(0,1) independently of (X,T). Thus, all three components of X are prognostic, and both components of V are effect-modifying. The value of βV determines the relative strength of the two markers, and the choice of ξ has a direct impact on δ(V), which describes the heterogeneity of the treatment effect. One thousand samples of size n=500 are generated independently and analyzed for marker comparison and combination.

The results for marker comparison are presented in Table 5, where γj (j=1,2) denotes the γ-value for the jth marker in V. Here, βV=(1,1)′ represents the null case (γ1=γ2) and βV=(1,2)′ represents the alternative case (γ1<γ2). The γ’s are estimated with the same choices of A=aˆ(X) as before (0, Eˆ(Y), Eˆ(Y|V) and Eˆ(Y|X)), expect that E(Y|V) and E(Y|X) are now estimated under linear regression models with V and X, respectively, as linear terms (in addition to an intercept). The results in Table 5 are consistent with those in Table 3 in the sense of small bias in point estimation and variance estimation as well as good coverage of confidence intervals. In addition, Table 5 demonstrates that the prognostic ability of X translates into improved precision when A changes from Eˆ(Y) to Eˆ(Y|V) and further to Eˆ(Y|X). Note that there is no improvement from A=0 to A=Eˆ(Y) because E(Y)=0 in this simulation study.