ROC curve inference for best linear combination of two biomarkers subject to limits of detection

The use of biomarkers to assist medical decision making, the diagnosis and prognosis of individuals with a given disease, is increasingly common in both clinical settings and epidemiological research. This has spurred an increase in exploration for and development of new biomarkers. All biomarkers, established and emerging, are limited by the sensitivity of the measurement instrument which can lead to censored observations, often more frequently in emerging biomarkers as laboratory methods are underdeveloped. This censoring process is due to a limit of detection (LOD) or the inability of an instrument to reliably measure samples below (or in some cases above) some threshold. Polychlorinated biphenyl (PCB) levels fit this scenario well as they have been linked with several adverse outcomes and their measurements are affected by LODs. A common statistical tool used to evaluate the utility of a potential biomarker such as PCBs is the receiver operating characteristic (ROC) curve. Consider populations of diseased and non-diseased people with levels of a specific biomarker denoted by independent random variables X and Y , respectively, with cumulative distribution functions F ( x ) and G ( y ). Suppose, the biomarker is utilized as an indicator of disease status where a level above some cut point, c , indicates a positive test for the disease and a level below c corresponds to a negative test. The sensitivity (the true positive test rate) and specificity (the true negative test rate) of the biomarker for a given c are q ( c ) = 1− F ( c ) and p ( c ) = G ( c ), respectively. The ROC curve is then a mapping of {1 − p ( c ), q ( c )} across all possible c . Proposed uses for the ROC curve ( Zhou et al., 2002 ; Pepe, 2003 ) include assessing discriminatory ability over all c , over a specific range of q ( c ) or p ( c ) and the maximum ability to differentiate between the populations. The area under the ROC curve, denoted here by AUC, is the most commonly used summary measure ( Zhou et al., 2002 ; Pepe, 2003 ) and has been shown to be P( X > Y ) for continuously measured biomarkers. AUC tends to range from 0.5 to 1 with larger values indicating greater separation between diseased and non-diseased biomarker levels. As a result, given multiple markers for the same disease, a researcher would be inclined to choose the one with the highest AUC to aid in decision making. Another approach would be to utilize a linear combination of these multiple biomarkers in lieu of selecting one. The attractiveness of a linear combination lies in being able to better discriminate, achieve a higher AUC, than if using any single biomarker alone. Now consider the case where two biomarkers, X⃗ = ( X 1 , X 2 ) T and Y⃗ = ( Y 1 , Y 2 ) T , measured in individuals with and without a disease, respectively, are independent bivariate normally distributed. Clearly, each biomarker considered individually would be normally distributed and AUC’s for each could be calculated and contrasted. The alternative described previously would be to use a linear combination, say U = β⃗ T X⃗ = β 1 X 1 + β 2 X 2 and V = β⃗ T Y⃗ = β 1 Y 1 +β 2 Y 2 , as a composite “biomarker” for decision making. Conveniently, U and V are also normally distributed with an ROC curve that is now a function of the choice of β⃗ T and the corresponding AUC can be denoted by AUC β . Linear combinations for binary outcomes are often estimated using logistic regression for a given set of covariates. However, given this scenario of having several biomarkers following a multivariate normal distribution, Su and Liu (1993) developed a best linear combination (BLC), β → 0 T , leading to AUC 0 , that is “best” with respect to maximizing AUC β over all real β⃗ T directly rather than maximizing the logistic regression model which, can be concordant or discordant depending on the scenario ( Pepe et al., 2006 ). This would allow for better discrimination than from any individual biomarker. When normal parameter values are unknown, we can use random samples to calculate estimates β → ^ 0 T and AÛC 0 via maximum likelihood techniques. One common complication in biomarker evaluation is that for a variety of reasons, random samples X and Y of biomarkers are often evaluated with non-detects or missing data below some LOD, quantified as d , effectively censoring the data below d . Omitting these values and proceeding with a complete case analysis has been shown to lead to biased AÛC for univariate normal and gamma distributed biomarkers ( Perkins et al., 2007 ). Categorizing the missing biomarker levels as ties at the lowest score, the standard non-parametric ROC curve and empirical AUC yield unbiased estimates of the biomarkers’ effectiveness given the measurement limitations. However, when underlying discriminatory ability is of interest, parametric methods can be used to estimate the ROC curve below the LOD and thus an AUC for a latent variable that might be measured completely. Generally, substituting a replacement value such as 0, d /2, d / 2 and d for unobservable data have been shown as a simple method to lessen biased parametric estimation (e.g. mean, variance and potential AUC) but the magnitude and direction of the remaining bias is highly dependent on the value chosen as well as the parameter being estimated ( Hornung and Reed, 1990 ). It has been shown ( Haas and Scheff, 1990 ; Lyles et al., 2001 ; Singh and Nocerino, 2002 ; Lynn, 2001 ; Perkins et al., 2007 ) that all of these methods lead to biased estimation of mean and standard deviation parameters, regression coefficients, odds ratio and AUC for a single biomarker following common distributions. In the logistic regression framework which has similarities here, Lynn (2001) demonstrated that the bias of these replacement values lessens the closer they are to the expected value below d , which performed similarly to more sophisticated methods of estimation. Lyles et al. (2001) proposed a more thoughtful solution of constructing a likelihood function for two censored bivariate normally distributed random variables. Based on this likelihood and under the assumption of normality, subsequent maximum likelihood estimators (MLEs) are efficient and asymptotically unbiased point estimators. In this paper, we propose in Section 2 to obtain MLEs for normal distribution parameters using Lyles et al.’s (2001) bivariate normal developments for parameter estimates based on a sample with left censoring in order to construct β → ^ 0 T and AÛC 0 to estimate that underlying potential of the BLC of two biomarkers. In Section 3, we extend Lyles et al. (2001) and consider the asymptotic distributions of parameter estimates by constructing the Fisher Information matrix for this case. Also in Section 3, these developments are subsequently used in finding the asymptotic distributions of point estimates β → ^ 0 T and AÛC 0 , leading to accompanying confidence intervals (CIs). This procedure allows for the estimation of the BLC’s underlying discriminatory ability or the potential that could be realized if it were possible to eliminate the LOD and censoring in the tail. This would be especially useful when researchers are exploring biomarkers for a given outcome using more cost-effective but less-sensitive assays, with the idea of further measurements or re-measurement being conducted on a narrowed set of promising biomarkers using a more sensitive and costly, “gold standard” assay. When a more sensitive assay is lacking this would help identify biomarkers with potential that are worth additional resources in refining a measurement process. However, the discriminatory ability of the BLC of these biomarkers as measured with LOD is better estimated using the AUC for traditional empirical ROC curve, denoted AŨC 0 here, which appropriately accounts for the censored values by essentially treating them as ties because they are indiscernible from one another. In Section 4, simulation is used to assess β → ^ 0 T and point estimators AÛC 0 and AŨC 0 as well as CIs for AUC 0 . An example in Section 5 using levels of the environmental toxicants PCBs to classify women with endometriosis is used to illustrate empirical and maximum likelihood techniques. We end with a brief discussion of issues surrounding estimation based on two biomarkers affected by LODs.

Suppose that pairs of biomarkers’ levels X⃗ = ( X 1 , X 2 ) T and Y⃗ = ( Y 1 , Y 2 ) T for cases and controls, respectively, are independent and have bivariate normal distributions f ( X⃗ ; μ⃗ X , Σ X ) and g ( Y⃗ ; μ⃗ Y , Σ Y ), respectively. These distributions can be written in a matrix form or explicitly expanded g ( W → ; μ → , ∑ ) = 1 2 π ∣ ∑ ∣ 1 / 2 exp { - 1 2 ( W → - μ → ) T ∑ - 1 ( W → - μ → ) } = 1 2 π σ 1 σ 2 1 - ρ 2 exp { - 1 2 ( 1 - ρ 2 ) ( η 1 2 - 2 ρ η 1 η 2 + η 2 2 ) } , where η l = ( w l − μ l )/σ l , l = 1,2, μ⃗ = (μ 1 , μ 2 ) T is the mean vector and the covariance matrix Σ consists of variances, ∑ l l = σ l 2 and covariance terms, Σ 12 = Σ 21 = ρσ 1 σ 2 . For ease of development, we will exclusively work with the latter, explicit form. Under this assumption Su and Liu (1993) showed that the formulae for coefficients leading to the BLC’s ROC curve is (1) β → 0 T ∝ ( μ → X - μ → Y ) T ( ∑ X + ∑ Y ) - 1 , which in turn leads to, Now suppose that the biomarker levels are measured with fixed LODs d⃗ = ( d 1 , d 2 ) T . Let the measured observations, Z⃗ X and Z⃗ Y , be the componentwise transformation of X⃗ and Y⃗ , respectively, such that Z X l = { X l ; X l ≥ d l , N / A ; X l < d l , and Z Y l = { Y l ; Y l ≥ d l , N / A ; Y l < d l , where for a fixed d l , l = 1, 2, the l- th biomarker level is either quantified or not. Without loss of generality we assumed that both diseased and non-diseased biomarker measurements were affected by the same point of censoring, say d Xl = d Yl = d l . Lyles et al. (2001) considered the case of two censored bivariate normally distributed random variables and developed the likelihood (3) L ( μ → , ∑ ; Z → W ) = ∏ Q → = ( 1 , 1 ) 1 σ 1 σ 2 1 - ρ 2 ϕ ( η 1 j ) ϕ ( η 2 ∣ 1 j ∗ ) × ∏ Q → = ( 1 , 0 ) 1 σ 1 ϕ ( η 1 j ) Φ ( η d 2 ∣ 1 j ∗ ) × ∏ Q → = ( 0 , 1 ) 1 σ 2 ϕ ( η 2 j ) Φ ( η d 1 ∣ 2 j ∗ ) × ∏ Q → = ( 0 , 0 ) Φ ( η d 1 j ) Φ ( η d 2 ∣ 1 j ∗ ) , where j = 1, …, n , η 2 ∣ 1 j ∗ = ( w 2 j - μ 2 - ρ σ 2 / σ 1 ( w 1 j - μ 1 ) ) / ( σ 2 1 - ρ 2 ) , Q⃗ is a vector of indicator functions with Q l = 1 if w l ≥ d l and Q l = 0 otherwise, and ϕ and Φ are the univariate standard normal pdf and cdf , respectively. Given random samples, Z X and Z Y , of two biomarkers levels measured with LODs in n X and n Y individuals, respectively, we can maximize Eq. (3) in order to generate MLEs for underlying normal parameters, θ⃗ = (μ⃗, Σ ) = (μ 1 , σ 1 , μ 2 , σ 2 , ρ), for cases and for controls. Substituting these estimators for the appropriate parameters in Eqs. (1) and (2) , the MLE’s β → ^ 0 and AÛC 0 are formed based on samples with multiple biomarkers affected by LODs.

Previously, the MLE AÛC , for a single biomarker affected by an LOD ( Perkins et al., 2007 ), was developed along with 1−α level CI formed using the asymptotic properties of AÛC . Applying those developments to A U ^ C 0 = AUC 0 ( θ → ^ X , θ → ^ Y ) results in n X + n Y ( A U ^ C 0 - AUC 0 ) ∼ . N ( 0 , σ A 2 ) where ~˙ denotes the asymptotic distribution. Again similar to the univariate case, the variance σ A 2 is obtained by the standard delta method (4) σ A 2 = λ - 1 ( ∂ AUC 0 ∂ θ → X ) T ( lim n X → ∞ n X Cov ( θ → ^ X ) ) ( ∂ AUC 0 ∂ θ → X ) + ( 1 - λ ) - 1 ( ∂ AUC 0 ∂ θ → Y ) T ( lim n Y → ∞ n Y Cov ( θ → ^ Y ) ) ( ∂ AUC 0 ∂ θ → Y ) with λ = n X / ( n X + n Y ) ∈ (0, 1) as n X → ∞ and n Y → ∞ and (∂AUC 0 /∂θ⃗) ij = ∂AUC 0 /∂θ⃗ ij being the ij- th element of (∂AUC 0 /∂θ⃗). The covariance matrices of MLE’s of unknown vectors of parameters are evaluated by the inverse of the Fisher information matrices, say I X , I Y , i.e. Cov(θ⃗ X ) = [ I X ] −1 and Cov ( θ → ^ Y ) = [ I Y ] - 1 . The asymptotic properties of β → ^ 0 mirror those for AÛC 0 with distribution n X + n Y ( β → ^ 0 T - β → 0 T ) ∼ . N p ( 0 → , ∑ β ) , where 0⃗ = (0, 0) T and the standard delta method is used to obtain the 2×2 covariance matrix (5) ∑ β = λ - 1 ( ∂ β → 0 T ∂ θ → X ) T ( lim n X → ∞ n X Cov ( θ ^ X ) ) ( ∂ β → 0 T ∂ θ → X ) + ( 1 - λ ) - 1 ( ∂ β → 0 T ∂ θ → Y ) T ( lim n Y → ∞ n Y Cov ( θ ^ Y ) ) ( ∂ β → 0 T ∂ θ → Y ) with ( ∂ β → 0 T / ∂ θ → ) i j = ∂ β → 0 j T / ∂ θ → i being the ij -th element of ∂ β → 0 T / ∂ θ → . To calculate CI’s for the elements of β → ^ 0 and AÛC 0 , we must find Σ β and σ A 2 , respectively, which consist of the covariance of the MLE’s θ → ^ X and θ → ^ Y and the partial derivatives of β → 0 T and AUC 0 with respect to each parameter. Here we are considering the case where two biomarkers’ levels are bivariate normally distributed with unknown parameters μ⃗ Y and Σ Y . The covariances, Cov ( θ → ^ X ) and Cov ( θ → ^ Y ) , can be determined by (6) Cov ( θ → ^ ) = Cov ( μ ^ 1 , σ ^ 1 , μ ^ 2 , σ ^ 2 , ρ ^ ) = [ I ] - 1 = [ I 11 ⋯ I 15 ⋮ ⋱ ⋮ I 51 ⋯ I 55 ] of bivariate normal distributions with LODs, where The details of Eq. (6) and the partial derivatives of β → 0 T and AUC 0 for Eqs. (4) and (5) are found in the Supporting Information. When all or a portion of the parameters are unknown, the MLE’s μ → ^ X , Σ̂ X , μ → ^ Y and Σ̂ Y are substituted for the appropriate parameters and generate approximate variances Σ β and σ A 2 for β → ^ 0 T and the AÛC 0 , respectively. We can then use this distributional information to assess the variability in our BLC coefficients, β → ^ 0 T , and to approximate α-level CI’s for AUC 0 by A U ^ C 0 ± z α / 2 σ A 2 / n X + n Y .

Utilizing the R programming language, we simulated B = 2000 data sets of biomarker values via the mvtnorm package ( n Y = n X = 50, 100, 200) for non-diseased and diseased individuals from various bivariate normal distributions with θ⃗ Y = (μ⃗ Y , Σ Y ) = (μ Y 1 , σ Y 1 , μ Y 2 , σ Y 2 , ρ Y ) = (0, 1, 0, 1, 0) and θ⃗ X = (μ⃗ X , Σ X ) = (μ X 1 , σ X 1 , μ X 2 , σ X 2 , ρ X ), ( σ X 1 2 , σ X 2 2 ) = ( 1 , 1 ) , ( 0.5 , 0.5 ) , ( 1 , 0.5 ) , mean μ⃗ X corresponding to AUC 0 =0.6, 0.7, 0.8, 0.9 and ρ X =0, 0.2, 0.5, 0.8. For each simulated group of samples, μ → ^ X , Σ̂ X , μ → ^ Y and Σ̂ Y were calculated by maximizing Eq. (3) via the optim function with method= “L-BFGS-B” allowing for bounded results necessary for the variance/covariance terms. Using Eqs. (1) and (2) , these estimators then lead to β → ^ 0 T and subsequently AÛC 0 . Additionally, we calculated the empirical AŨC 0 using the same β → ^ 0 T and replacing missings with the often used replacement values a l = d l , l =1, 2, for each biomarker although any value 0 ≤ a ≤ d will result in the same AŨC 0 for these positive β → ^ 0 T . Using methods laid out in Section 3 and bootstrapped quantiles (1000 resamplings using the boot function), 95% CI’s for AUC 0 accompany AÛC 0 and AŨC 0 , respectively. We calculated these estimators first with no LOD, no missingness, to appreciate the baseline properties of the estimators under these various conditions and then applied increasing levels of censoring, where values of d were chosen to be the 20, 40, 60 and 80th quantiles of marginal distributions of the non-diseased population. The β → ^ 0 T themselves were assessed and, while bias increased slightly as missingness increased, the proportionality of the coefficients’ bias remained consistent with the proportion of the true coefficients, e.g. β̂ 01 /(β̂ 01 +β̂ 02 ) ≈ β 01 /(β 01 +β 02 ). This is an excellent result as Eq. (1) only requires proportionality in the coefficients. As expected, the estimated relative bias, (AÛC 0 − AUC 0 )/AUC 0 , and root mean squared error (RMSE) of AÛC 0 decrease as sample sizes increase. The ranges of relative bias were 0.0027 to 0.1463, 0.0019 to 0.0880, 0.0003 to 0.0446 and were 0.0245 to 0.1504, 0.0174 to 0.0973, 0.0121 to 0.0618 for RMSE, respectively for n Y = n X =50, 100, 200. Figure 1 depicts the relative bias and RMSE for ( σ X 1 2 , σ X 2 2 ) = ( 1 , 1 ) ; the ( σ X 1 2 , σ X 2 2 ) = ( 0.5 , 0.5 ) and (1, 0.5) cases display similar relations in direction and magnitude regarding increasing sample size and percent missing. Relative bias in Fig. 1 clearly increases as ρ X and percent missing increase and decreases as a function of sample size and AUC 0 . Changes in relative bias are slight for all levels of ρ X and from 0 to 40% missing, increasing moderately at 60% and only substantially from 60 to 80% missing. Generally, AÛC 0 generated from only n x = n y =100, with over half of those values missing, has bias that is less than two percent of AUC 0 , which is generally equivalent to the relative bias of AÛC 0 based on a full data set. The second row of plots in Fig. 1 show that RMSE of AÛC 0 has generally the same relationship with sample size, correlation, percent missing and level of discrimination as with relative bias. Again, AÛC 0 calculated using only n x = n y =100, with over half of those values missing, has relative RMSE around one percent, which is comparable to that of AÛC 0 based on all the data. The empirical estimator AŨC 0 , is obviously unbiased for the discriminatory ability of the BLC of the biomarkers as measured with ties for missing values below the LOD. However, Fig. 2 displays a plotting of pairs of relative bias of AÛC 0 and AŨC 0 with regard to the underlying latent AUC 0 for the various scenarios with n x = n y =50. This depiction shows, relatively, how consistently AŨC 0 can estimate AUC 0 with the β → ^ 0 T and the current measurements. Toward this end, scenarios of relative bias of AŨC 0 close to the horizontal dashed line indicating consistent estimates range from higher levels of AUC 0 , (triangles=0.8, cross=0.9) and for 0, 20 and 40% missingness and also for lower AUC 0 with 40 and 60% missingness. In the latter instance, AŨC 0 displayed better consistency for AUC 0 than AŨC 0 for these small sample sizes. Figure 2 also shows the potential bias of using AÛC 0 in lieu of AŨC 0 when measurement cannot be improved. The coverage probability of 95% CI’s that accompanied AÛC 0 were nominal or near nominal coverage for all AUC 0 ranging from 0.890 to 0.964, 0.913 to 0.965 and 0.933 to 0.960 for n x = n y =50, 100 and 200, respectively. No discernable patterns exist regarding ρ X =0, 0.2, 0.5, 0.8 or the percent missing due to the LOD. It should be noted that unlike Reiser and Faraggi’s (1997) CI for AUC using complete data, our CI’s do not take into account the bounded nature of the AUC and coverage may be adversely affected near these bounds. However, our simulation did not reflect any difficulties for the sample sizes and levels of discrimination shown. Figure 3 displays the width of MLE and bootstrapped CI’s, open and solid points, respectively, versus their coverage probability for our smallest sample size, n x = n y =50. Generally, CI’s increase in width to reflect uncertainty from increasing missingness, where larger point size corresponds to a larger proportion missing. Additional plots for the larger sample sizes displaying increased coverage and tighter widths are included in the Supporting Information.

To illustrate our method, we use PCBs, environmental toxicants, as potential indicators of endometriosis ( Louis et al., 2005 ). Endometriosis is a gynecological disease exclusive to species that menstruate such as humans and other primates, occurring predominantly in women of reproductive age. Data from experimental studies in animals and observational human studies suggest an association between dioxin and PCBs and endometriosis. In our data, PCBs 153 and 180 were measured in 28 women with and 51 women without endometriosis. The biomarkers were measured jointly and were expected to be correlated. However, the sensitivity of the measurement process differed for each biomarker with PCB 153 having an LOD of d 153 =0.2, resulting in 64% of the cases and 74% of the controls having unobservable levels and PCB 180 having an LOD of d 180 =0.034, resulting in no missing cases but 11% of the controls having unobservable levels ( Whitcomb et al., 2005 ). Analyzing the biomarkers univariately, empirical methods led to AŨC 153 = 0.564 and AŨC 180 = 0.609. As measured, PCB 180 and 153 appear to have some discriminatory ability for women with endometriosis, greater for PCB 180. To investigate what potential might lie below the LODs, we can employ parametric techniques. Univariate normal distributions were assessed on the log transformed biomarkers similar to Lyles et al. (2001) via q–q plots. Figure 4 displays q–q plots where both biomarkers fit log normal distributions well, as evidenced by diagonal points corresponding to the points above the LODs. The horizontal points are essentially quantile place holders for the observations below the LOD. Assuming that the data we don’t see follow the data we do see, univariate normal likelihoods led to the MLE’s AÛC 153 = 0.511 and AÛC 180 = 0.630 for PCB 153 and PCB 180, respectively ( Perkins et al., 2007 ). Again, PCB 180 displays a moderate ability to differentiate but the potential of PCB 153 seems to have been diminished. The result AŨC 153 >AÛC 153 is due to a large “hook” in the ROC curve ( Pesce and Metz, 2007 ) resulting from disparate variances. This behavior is rarely reflective of actual etiology but rather is likely to be a result of highly variable parameter estimates given the small sample sizes and large amount of missings. Regardless, PCB 153 might be discarded as potentially lacking discriminatory ability for endometriosis. However, these two biomarkers are closely linked biologically and are suspected as being highly correlated as a result of being in a mixture of environmental exposures. Using a variety of methods, we considered the joint distribution of PCB 153 and 180 to investigate whether PCB 153 might have some potential as a contributor in a BLC. First, Eq. (3) was maximized using log transformed PCB levels, MLEs for distribution parameters were calculated as well as β → ^ 0 T = { β 153 , β 180 } = { - 1.73 , 1.65 } via Eq. (1) . Using Eq. (2) , we estimated that the potential from this BLC if measurements could be improved to be AÛC 0 = 0.747 with 95% CI (0.580, 0.915). This result shows that the BLC of PCBs 153 and 180 could be far superior to either biomarker alone and is potentially a good discriminator of women with and without endometriosis if the measurement sensitivity of PCB 153 could be improved. However, the sample sizes here have led to a significant but wide CI which does not rule out AUC 0 <0.630. Table 1 compares AÛC 0 to several alternative estimators of the potential AUC 0 . Naïve methods to estimate potential AUC 0 in Table 1 shows the results of using replacement values (e.g. imputing a l = d l /2 for values below the LODs) and standard parametric methods to obtain naïve BLC, β → ^ 0 R T . Imputing values in this fashion clearly violates the assumption of bivariate normality. However, we generated such estimates for the sake of comparison since standard methodology is often applied to data sets with imputed values. While these results are a slight improvement over the univariate biomarkers, it is vastly less than the proposed AÛC 0 . These naïve cases also point to widely varying linear combinations, as in column 3 of Table 1 . The second set of estimates in Table 1 begins with this BLC’s current diagnostic ability, where AŨC 0 was estimated empirically using various replacement values ( d , d/ 2 and a l =E[ X l | X l < d l ], the expected values below the LODs) and the mles β → ^ 0 T estimated previously. These estimates vary but all are less than AŨC 180 and even AŨC 153 . Alternatively, we could use these replacement values in conjunction with simple logistic regression as a means to identify the “best” coefficients to estimate AŨC 0 as measured. Lynn (2001) showed that in some cases these naïve methods can provide logistic regression coefficient estimates comparable to those using a full maximum likelihood approach for logistic regression accounting for an LOD. These “logistic” AŨC 0 in Table 1 are on par with those using the β → ^ 0 T , “mle,” likely due to the logistic coefficients having the same signs and not grossly different proportionality. However, they do shift some of the focus to PCB 180 as evidenced by the slightly larger proportion of the linear combination, making the subsequent AŨC 0 better than AŨC 153 . The ordering AŨC 153 <AŨC 0 <AŨC 180 is due to the difference in signs for the coefficients and in percent below the LOD for the two biomarkers. If both coefficients were positive, the order of PCB 180 would be preserved when ties from PCB 153 are introduced to a BLC and thus AŨC 0 ≥ AŨC 180 . As measured, the BLC of PCBs 153 and 180 would not be of use beyond what PCB 180 can already achieve, however, Fig. 2 shows us that we need only to reduce the missingness in PCB 153 to between 20 and 40% to realize most of the difference between the potential latent and as measured effectiveness.

While multiple biomarkers are often available or obtainable for the same outcome, investigators often ignore all but one and in doing so essentially throw out potential ability to discriminate beyond the biomarker they have chosen. For this reason and where biomarkers can be assumed to be normally distributed, the proposed BLC should be used in lieu of a single marker, especially in cases where multiple biomarkers come at no additional cost, as is the case with multiplex assays. Using two biomarkers measured with LODs can result in one or both measured with censored values for an individual. The methods developed here allow us to properly account for the missing values while simultaneously taking advantage of the benefits in discriminatory ability realized by a BLC. While replacement has been shown previously to be useful in certain situations as an ad hoc approach for the AUC in the univariate case, we’ve shown here that standard parametric methods with replacement values can yield estimates of the potential AUC 0 that vary greatly from the MLE ( Perkins et al., 2007 ). For this reason, replacing missing values by a constant is not recommended when estimating, even naïvely, the potential of a BLC of biomarkers. The method proposed here is related to the approaches for handling censored covariates in logistic regression setting. Lynn (2001) discusses an approach for estimating the linear combination of subject-to-LOD covariates for prediction of a binary outcome. Due to the invariance of the ROC curve with respect to order-preserving transformations and because of the assumed multivariate normality of the covariates, the potential AUC 0 for Lynn’s linear combination can be computed using the standard formulation for binormal AUC. The novelty of the approach developed here is that the proposed BLC targets maximization of the AUC 0 rather than maximization of the likelihood under the logistic regression model. Depending on the parameters of the distributions the results can be either similar or substantially different ( Pepe et al., 2006 ). While the likelihood ratio has been demonstrated to achieve the highest possible AUC for a given set of biomarkers, in the case of two bivariate normally distributed biomarkers the absolute difference in estimated AUC can be minimal. However, this is not always the case as likelihood ratio ROC curves are assuredly “proper” or completely concave while conventional “binormal” ROC curves and those based on linear combinations of normal biomarkers can have minimal to severe non-concave portions or “hooks” (e.g. σ X 2 ≫ σ Y 2 , (β⃗ T Σ X β⃗)≫(β⃗ T Σ Y β⃗) or b ≪1 from Pesce and Metz (2007) ). The linear combination here has the benefit of being simple to interpret and can incorporate the empirical estimate of the biomarkers’ effectiveness as currently measured where the likelihood ratio does not. More importantly though is that in the context of the latent discriminatory ability, both would work similarly in identifying biomarkers with potential that could be measured more precisely and more extensively in the future. While, Section 3 considers the asymptotic properties of AÛC 0 for the BLC of two multivariate normal biomarkers, Section 4 shows small sample performance that is nearly unbiased and achieves nominal coverage probability for sample sizes as little as 100 with 20, 40 and 60% missing values due to LODs. The need for further investigation regarding magnitude and potential causes of differences between the likelihood- ratio and conventional “binormal” ROC curves in the context here remains. As with all parametric estimation, this BLC of two biomarkers corrected for LODs is dependent on distributional assumptions specifically that the biomarkers or their transforms are multivariate normally distributed. The literature ( Molodianovitch et al., 2006 ; Perkins et al., 2007 ) has shown MLEs of AÛC for the univariate case are robust to minor deviations from the normal assumption and unfavorable results when assumptions are grossly violated. Estimates of the AUC based on replicates of a single biomarker, a special case of the multivariate normal developments here, have shown favorable robustness by generating biomarker levels from gamma distributions ( Perkins et al., 2009 ). While the bias was on average less than one percent of the AUC being estimated, caution and diligence should still be taken when employing this BLC resulting in AÛC 0 . A biomarker may provide differential discrimination based on some associated factor. The inclusion of this covariate information can of course be incorporated by conducting a stratified analysis for each level of such a factor. A better approach might be to directly include the covariate information in the likelihood function. In our example we applied several methods to empirically estimate the as measured discriminatory ability, AŨC 0 , by using simple replacement values in standard logistic regression or with β → ^ 0 T . Lynn (2001) showed that in some cases, naïve replacement could yield similar logistic regression coefficient estimates to those estimated using maximum likelihood accounting for the LOD properly. We found a similar concordance here in AŨC 0 by virtue of the naïve logistic coefficients and our β → ^ 0 T yielded similarly proportional estimated coefficients. However, instances where this is not the case surely exist and care should be taken to investigate potential discordance before applying either BLC in practice. These methods for estimating AUC 0 while properly accounting for missing data censored below LODs will provide an asymptotically unbiased view of the potential discriminatory ability of the BLC of a set of two biomarkers. We must note again that the AÛC 0 in this scenario reflect potential discriminatory ability and cannot be realized in practice until observations can be measured and thus allow for differentiation below the LODs. This is often the case where a less sensitive, but cheaper assay could be conducted to narrow down potential biomarkers before conducting more expensive but more precisely measured assays. We have also demonstrated the usefulness of the empirical estimate of what is possible “as measured” and how the majority of two biomarkers’ latent discriminatory ability might be realized by improving measurement, reducing missingness to between 20 and 40% rather than to zero. Estimating this potential is increasingly important as it allows researchers to focus limited resources to improve the measurement process in biomarkers that display promise and subsequently lead to improved diagnostic care.