Generalized ROC curve inference for a biomarker subject to a limit of detection and measurement error

Similar to previous asymptotic developments, the MLE's μ ^ X , μ ^ Y , σ ^ X 2 , σ ^ Y 2 , σ ^ ε 2 , and AÛC constructed here are asymptotically normally distributed and unbiased [ 9 , 15 ]. As ultimately the AUC is our parameter of interest we will focus on the distribution of AÛC, specifically N ( A U ^ C − AUC ) ∼ ⋅ Normal ( 0 , σ A 2 ) , where N = n X + n Y and σ A 2 is the variance parameter. The development of σ A 2 here will mirror those in the appendix of Perkins et al. differing only in that Fisher's information matrix, denoted by I , for μ ^ X , σ ^ X and μ ^ Y , σ ^ Y , pair wise [ 9 ]. The derivatives and subsequent expected values necessary for I here can be achieved by accounting for the specific covariate structure in those developed for the general lower censored bivariate normal [ 15 ]. The details of these developments are left to the appendix . Utilizing these distributional properties, α level two-tailed CI for AUC is constructed from A U ^ C ± z α ∕ 2 σ A 2 N As is often the case, when σ A 2 is unknown, an approximate α level CI for AUC is formed by using the estimator σ ^ A 2 in lieu of σ A 2 .

A simulation study assesses the estimators developed in Sections 2 and 3, MLE for AUC that corrects for both an LOD and random additive ME. The AÛC was assessed using bias and root mean-squared error, RMSE, over B = 2000 iterations. Coverage probabilities of 95 per cent CIs for AUC's were recorded and examined as well. Sample sizes of varying magnitude, n X = n Y = 50, 100, 1000, were used to show the asymptotic nature of these estimators. Simulated true non-diseased values were generated from a standard normal distribution, μ Y = 0 and σ Y = 1. True diseased values were generated from a normal distribution with parameters σ X 2 = 1.0 or 0.5, equal and unequal variance cases, respectively, and mean μ X corresponding to AUC = 0.6, 0.7, 0.8, 0.9. ME was then introduced from independent mean zero normally distributions with variance σ ε 2 . Error variability was linked to inherent biomarker variability, σ Y 2 and σ X 2 , through the ‘reliability index’ extended to two populations, [ 16 ] R = σ X 2 + σ Y 2 σ X 2 + σ Y 2 + 2 σ ε 2 R is an index commonly used in the social science literature to reflect the proportion of observed variability that is due to the actual population [ 17 ]. Of course, 1 – R is the amount of observed variability that is due to error. Several scenarios of increasing error variability were considered R = 0.90, 0.75, 0.50, 0.25. As the method in Section 3.2 also corrects for an LOD, we set d based on the desired per cent, 0, 10, 25, 50, 75 per cent, of the non-diseased population, observed with ME, that would be unquantifiable, missing. Tables I – VI show the results of our simulations. Bias of the equal and unequal variances cases are displayed in Tables I , II for every combination of sample size, true AUC, R and LOD. Tables III , IV and V, VI are for the same parameters but display RMSE and coverage probability of 95 per cent CIs, respectively. Bias and RMSE levels were generally consistent between Tables I , II and III , IV , respectively, or for equal and unequal variance cases. The ranges of bias were –0.0175 to 0.0143, –0.0157 to 0.0103 and –0.0048 to 0.0089 for n X = n Y = 50, 100, 1000, respectively, with mean bias across scenarios being approximately 0.25 per cent of the AUC being estimated. Overall RMSE levels ranged from 0.0009 to 0.0346, 0.0005 to 0.0156 and 0.00003 to 0.0012 for n X = n Y = 50, 100, 1000, respectively, with corresponding mean RMSEs 0.0055, 0.0028 and 0.0003. As expected, bias and RMSE decrease with increases in sample sizes. Other relations between variations in distributional scenarios and bias, RMSE are better seen via Figures 1 and 2 , respectively. Figures 1 and 2 display bias and RMSE results as a function of missingness due to LOD (0, 10, 25, 50, 75 per cent) where different lines correspond to differing levels of ME quantified by R (0.90 to 0.25 corresponding to low to high levels). Figures are for n X = n Y = 100 only, however, patterns are consistent over sample sizes and are excluded only for clarity of exposition. Figure 1 shows that the bias of AÛC is no more than marginally dependent for levels of R = 0.90–0.50, up to half of the observed variability originating from ME. Only when we consider the extreme case of R = 0.25, where three-fourths of overall observed variability is due to random error, do we see an increase in the magnitude of bias. Bias is also relatively unchanged as a function of per cent missing (flat lines in figures) until after correcting for 50 per cent missing. Figure 2 displays a decrease in RMSE as a function of AUC. Similar to bias, RMSE increases slightly as a function of R , until R = 0.25, exemplified by the different lines increasing in overall magnitude. RMSE increases slightly as a function of per cent missing (left to right) up to 50 per cent after which a sharp increase in RMSE is observed. The increase in RMSE of estimators correcting for R = 0.90 and 75 per cent missing substantially from the level at 50 per cent. This is driven by the bias of correcting for 75 per cent missing data at these levels and not variance. The vast majority of coverage probabilities of 95 per cent CI's for AUC in Tables V and VI are nominal or with in normal variation for a simulation with 2000 iterations. CI's seem to be conservative for low levels of discrimination with large per cent missing (75 per cent) and under covers for AUC = 0.90 with 75 per cent. Again, increasing the number of replicates per individual may help improve coverage probability. Overall, point and CI estimators perform well while correcting for various levels of LOD and ME.

Pre-eclampsia affects about 5 per cent of pregnancies, resulting in substantial maternal and neonatal morbidity and mortality [ 18 ]. Identifying women suffering from pre-eclampsia is a very important step in the management of the disease. Vascular endothelia growth factor (VEGf) and its soluble receptor 1 (sFLT-1) are promising biomarkers for such classification. Apriori it was suspected that sFLT-1 and VEGf had substantial levels of ME. Subsequently, replicate samples were taken and analyzed. An additional source of error for both biomarkers were the manufacturer established LODs that each biomarker was subject to. In a nested case control study, using two replicates per n X = 344 the and n Y = 313 women, we evaluated discriminatory abilities of these biomarkers. In order to apply our method we assessed the normality of biomarker levels with Q–Q plots. As Figure 3 demonstrates, we found that a log transformation resulted in a distribution that reasonably reflects our assumption of normally distributed biomarker levels. After the log transformation, we maximized the likelihood, equation (1) , with respect to the normal parameters and subsequently generated the mean and covariance matrix estimators for sFLT-1 levels in of our transformed cases and controls μ ^ X = 7.137 , μ ^ Y = 6.844 Σ ^ X = [ 0.652 0.628 0.628 0.638 ] and Σ ^ Y = [ 0.431 0.404 0.404 0.410 ] and μ ^ X = 2.065 , μ ^ Y = 2.16 Σ ^ X = [ 0.428 0.338 0.338 0.354 ] and Σ ^ Y = [ 0.295 0.219 0.219 0.233 ] for VEGf levels, respectively. The covariance matrices lead to Reliability coefficient ≈0.97 for sFLT-1 and Reliability coefficient ≈0.77 for VEGf which corresponds to error variance making up roughly 3 and 23 per cent of the total biomarker variability, respectively. Fortunately in this instance, all sample measurements of sFLT-1 and VEGf registered above the LODs set by the test kit manufacturer resulting in no censored observations. The individual MLEs above lead to an AÛC = 0.613 and 95 per cent CI (0.583–0.643) for the true AUC for sFLT-1 and AÛC = 0.554 with 95 per cent CI (0.522–0.585) for VEGf as a discriminator of women with and without pre-eclampsia. Neither biomarkers have outstanding capabilities but both are significantly better than chance when corrected for ME. While no observations were censored, Figure 4 displays point (solid) and confident interval (dashed) estimators for AUC as if biomarker measurements had been affected at various levels of censoring do to the LODs. The bold lines in Figure 4 are based on the methods developed in Sections 2 and 3 utilizing both samples per woman. Point estimators for both sFLT-1 and VEGf are nearly unchanged correcting for increasing levels of missingness up to approximately 50 per cent missing as demonstrated by the nearly horizontal bold line. From 50 all the way to 80 per cent missing, AÛC becomes only slightly more variable. Bold dashed lines gradually widen as per cent missing increases to allow for additional variability due to estimation being based on fewer observed measurements but at a rate such that VEGf remains a significant discriminator up to about 35 per cent missing and sFLT-1 is significant for up to 75 per cent missing. Figure 4 also displays similar lines (non-bold) based on a single sample per woman and univariate LOD methods [ 9 ]. These lines illustrate that ignoring ME results in wider CIs and reduced point estimates AÛC across increasing levels of the LOD in both sFLT-1 and VEGf. The degree to which estimators are affected by ignoring ME is obviously a function of the level of error being ignored. Using a single measurement of sFLT-1 results in minimal if any detriment in estimation due to the low level of error being corrected for, 3 per cent of sFLT-1 variability. However, ignoring ME in VEGf measurements would result in substantially reduced AÛC and wider CI across all levels of per cent missing. These single measurement CIs for VEGf never display a discriminatory ability significantly better than chance. For sFLT-1 investigators would have achieved adequate results using only one sample per woman even though using two replicates resulted in a more consistent estimator and narrower CI across levels of per cent missing. The low level ME present for this to be possible was unknown apriori . VEGf levels illustrate the other possibility, where the level of error is a significant contributor in the overall variability of biomarker levels. Here a single sample yields substantially wider CIs and a point estimate displaying less discriminatory ability than estimators based on replicates and corrects for ME.

Diagnosis and prognosis of diseases are increasingly dependent on the effectiveness biological markers to discriminate individuals with and without disease. The receiver operating characteristic (ROC) curve is a statistical device commonly used to evaluate the utility of a potential biomarker. Consider populations of individuals with and without a particular disease. Let independent random variables X and Y , respectively, denote the levels of a specific biomarker of individuals within these populations which are governed by cumulative distribution functions F ( x ) and G ( y ). Suppose, the biomarker is thought to have utility as an indicator of disease status with levels above some cut point, c , indicating individuals testing positive for the disease and levels below c corresponding those individuals testing negative. Two quantities useful in assessing biomarker utility, sensitivity ( q ( c )) and specificity ( p ( c )), are then 1– F ( c ) and G ( c ), respectively, for a given c . The ROC curve is a mapping of ( p ( c ), q ( c )) across all possible c . Proposed uses for the ROC curve include assessing overall discriminatory ability, or over a specific range of q ( c ) or p ( c ) and the maximum ability to differentiate between the populations [ 1 , 2 ]. The area under the ROC curve, denoted here by AUC, is the most commonly used summary measure and is equal to P ( X > Y ) [ 1 , 2 ]. The AUC measures separation between diseased and non-diseased populations based on biomarker levels. Ranging from 0.5 to 1, larger AUC values indicate greater separation and thus greater utility as a diagnostic tool. A major complication in the development of emerging biomarkers is an investigator's ability to adequately assess biomarker utility from data with measurement error (ME). Conventional ROC methods assume that the biomarker levels measured are the true levels of those individuals. However, most biomarkers are subject to random additive ME or random variability from this true level of an individual. This variation could be due to intra-individual or measurement variability among other causes. Regardless of the origin, biomarker levels are often measured with error what is commonly assumed to be mean zero normally distributed error. Subsequently, measurements may be subject a limit of detection (LOD) where due to instrumental limitations measurements are undetectable below some level for that particular biomarker. Formally, in the presence of ME the ‘true’ value of the biomarker, X and Y , is not available. Instead we assume that actual observations are the true biomarker levels plus some additional error, ε X i , ε Y j ∼ N ( 0 , σ ε 2 ) , x i = X i + ε X , i , i = 1 , … , n x y i = Y i + ε Y , j , j = 1 , … , n y with X i , Y j , ε X,i and ε Y,j independent of each other. The equal variance assumption of ε X,i and ε Y,j is not necessary, but in many instances, justified because the variations leading to ME are due to factors connected with the laboratory or measurement process and do not depend on the individual's risk or disease status or the true biomarker value. When the true levels are normally distributed, adding independent normally distributed error results in levels that are also normally distributed with mean μ X or μ Y and variance σ X 2 + σ ε 2 or σ Y 2 + σ ε 2 , respectively, depending on disease status. Naively assuming that these are the true biomarker levels would result in an ROC curve and associated AUC estimate that is biased toward the null [ 3 ]. These negatively biased results hinder an investigator's ability to identify a biomarker as potentially useful and might lead to it being mislabeled useless for disease diagnosis. Instead, replicate measurements x ik and y jk where k = 1,..., p , are commonly used in order to correct for ME by estimating the amount of intra-individual variability present [ 3 , 4 ]. Estimation is often further complicated by the inability of laboratory methods or machinery to detect biomarker levels below some LOD, quantified as d , effectively censoring the data below d . That is to say that we observe z X , i = { X i + ε X , i if X i + ε X , i ≥ d N ∕ A if X i + ε X , i < d } and z Y , j = { Y j + ε Y , j if Y j + ε Y , j ≥ d N ∕ A if Y j + ε Y , j < d } Obviously, omitting these unquantifiable values yields a bias assessment of biomarker effectiveness. Substituting a replacement value, 0, d /2, d ∕ 2 and d for unobservable data has been suggested as an alternative to provide a less biased assessment than omission [ 5 ]. It has been shown 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 [ 6 – 9 ]. Maximum likelihood estimators (MLE) of ROC measures have been developed based on univariate and multivariate likelihoods functions for these censored data sets and have shown to be efficient and asymptotically unbiased [ 8 – 12 ]. However, none of these methods have considered estimation based on data affected by both. In this paper, we assume the availability of normally distributed biomarker levels, in replicate from populations with and without disease, affected by ME and LOD. Based on these data, we propose in Section 2 to generate MLEs for μ X , σ X 2 , μ Y , σ Y 2 and σ ε 2 from a likelihood function that blends the general form of Lynn's bivariate normal developments for parameter estimates based on a sample with left censoring and Schisterman et al. , Resier and Faraggi's likelihood for replicates of normal biomarker levels with ME [ 3 , 4 , 13 ]. Subsequently, an MLE AÛC will be generated based on MLEs of normal parameter that governs the distributions of true biomarker levels. In Section 3, the asymptotic distribution of AÛC is considered and a confidence interval (CI) is developed. Point and CI estimators are evaluated via simulation in Section 4. An example in Section 5 using polychlorinated biphenyl levels to classify women with endometriosis is used to illustrate the impact of replacement and maximum likelihood techniques to correct for ME and LOD. We conclude with a brief discussion of issues surrounding estimation based on data affected by ME and LOD.

Methods that correct for ME due so through some form of repeated measurements in order to assess σ ε 2 . Here, correcting for ME is only half the battle but we will continue to assume the availability of repeated measurements on every individual in our data set. Subsequently, we can form estimators corrected for both sources of error using maximum likelihood. Consider p = 2 repeated measurements for every individual in a study population of individuals with and without disease. Let x ik = X i + ε X,ik and y jk = Y j + ε Y,jk denote the k th, k = 1 or 2, measurement of the i th and j th individuals, respectively, and x i = ( x i 1 , x i 2 ) and y j = ( y j 1 , y j 2 ) are bivariate vectors of observations. Using replicates in order to correct for ME was considered separately by Reiser and Faraggi and later by Schisterman [ 3 , 4 , 13 ]. The correlated replicates of the j th non-diseased individual follow the bivariate distribution y j ~ N 2 ( μ Y , Σ Y ) , where μ Y = ( μ Y , μ Y ) and Σ Y = [ σ Y 2 + σ ε 2 σ Y 2 σ Y 2 σ Y 2 + σ ε 2 ] with correlation ρ = σ Y 2 σ Y 2 + σ ε 2 Replicates of biomarker levels in diseased individuals are similarly distributed. MLE's are generated by maximizing the standard likelihood function for binormal data, accounting for the singular mean and specific covariance structure. Now let all measurement be subject to an LOD. Vectors z → X , i = ( z X , i 1 , z X , i 2 ) and z → Y , i = ( z Y , i 1 , z Y , i 2 ) are observed replications on diseased and health persons, respectively, after the transformation in Section 1. The MLE's developed previously do not account for this additional source of error. A maximum likelihood method does exist for multiple correlated normally distributed biomarkers affected by LODs in order to generate estimates of the parameters μ 1 , μ 2 , σ 1 , σ 2 and ρ [ 8 , 14 ]. As the replicates are bivariate normally distributed, we can correct for the ME and LOD here by adapting this general likelihood by letting μ 1 = μ Y , μ 2 = μ Y , σ 1 = σ Y 2 + σ ε 2 , σ 2 = σ Y 2 + σ ε 2 and again, ρ = σ Y 2 ∕ ( σ Y 2 + σ ε 2 ) . Given the two by n Y matrix, Z Y , the general likelihood in the presence of ME and LOD is (1) L ( μ → Y , Σ Y ; Z Y ) = ∏ j = 1 n Y [ 1 σ ε σ ε 2 + 2 σ Y 2 ϕ ( η Y 1 j ) ϕ ( η Y 2 ∣ Y 1 j ∗ ) ] I ( Z Y 1 j , Z Y 2 j ≠ N ∕ A ) × [ 1 σ ε 2 + σ Y 2 ϕ ( η Y 1 j ) Φ ( η d 2 ∣ Y 1 j ∗ ) ] I ( Z Y 1 j ≠ N ∕ A , Z Y 2 j = N ∕ A ) × [ 1 σ ε 2 + σ Y 2 ϕ ( η Y 2 j ) Φ ( η d 1 ∣ Y 2 j ∗ ) ] I ( Z Y 1 j = N ∕ A , Z Y 2 j ≠ N ∕ A ) × [ ∫ − ∞ d 1 1 σ ε 2 + σ Y 2 ϕ ( η Y 1 j ) Φ ( η d 2 ∣ Y 1 j ∗ ) d y 1 ] I ( Z Y 1 j , Z Y 2 j = N ∕ A ) where η Y 1 j = y j 1 − μ Y σ Y 2 + σ ε 2 and η Y 2 ∣ Y 1 j ∗ = η Y 2 j − ρ η Y 1 j 1 − ρ 2 for the non-dieased parameters. Maximizing the log of Equation (1) for both the diseased and non-diseased biomarker levels yields MLE's μ ^ X , μ ^ Y , σ ^ X 2 , σ ^ Y 2 and σ ^ ε 2 . Using these MLE's, we can construct the MLE A U ^ C = Φ ( μ ^ X − μ ^ Y σ ^ X 2 + σ ^ Y 2 ) that is corrected for both ME and an LOD, that is to say AÛC estimates the biomarkers true underlying ability to discriminate individuals with and without the disease.

Previous methods have corrected for ME and LOD individually. This is the first development that corrects the binormal ROC curve and AUC for both simultaneously. Through maximum likelihood methods we demonstrated both asymptotically unbiased point estimates and CIs with nominal coverage. This work could be extended to other summary measures of the ROC curve such that the Youden index and associated optimal cut-point and the partial area under the curve. Additionally, a natural extension to the p = 2 case we have demonstrated would be to perform the steps in the Appendix similarly to achieve estimates utilizing more than two replicates. However, the greatest benefit in terms of ME correction is realized from going from one observation, no ability to correct for ME, to two. Adding an additional replicate per person would improve estimates of error variability but is unlikely to be sufficient to justify increasing costs for data acquisition by 50 per cent over using two replicates. The methods we have developed here first estimate the normal distribution parameters for diseased and non-diseased populations accounting for both ME and LOD through maximizing the appropriate likelihood function for these data. We then use standard binormal techniques and these corrected normal distribution parameters to achieve estimates for the ROC and AUC corrected for ME and LOD. This technique has the benefit of using the familiar binormal ROC curve and AUC with all the properties attached such as being invariant to monotonic transformations, but also the restrictive requirement that the data are normally distributed. As with all parametric methods, these estimators should only be used when distributional assumptions are adequately assessed and reasonably satisfied. We assessed robustness by using gamma distributed data via simulation similar to that conducted in Section 4. The gamma distribution was chosen for its non-symmetric shape. The distribution of controls was fixed with α Y = 4 and β Y = 1, the distribution of cases was adjusted to achieve AUC = 0.6, 0.7, 0.8, 0.9. The results were encouraging with the bias as a percentage of AUC averaged less than 1 per cent across sample sizes n = 50, 100, 1000 displaying that reasonable departures from normality do not lead to overly biased estimation. Coverage probably held up as well, averaging 97 and 96 per cent for sample sizes n = 50, 100, respectively, but deteriorated greatly at n = 1000 where coverage averaged 75 per cent, mainly do to AUC = 0.9, 51 per cent coverage. These results are to be expected as when n = 1000 one should be able to adequate address distributional assumptions as reasonably satisfied so this decreased coverage is unlikely to occur in practice. It has also been shown elsewhere that similar estimators are robust to small departures from normality but are unreliable for gross departures from such assumptions both [ 3 , 9 , 19 ]. Correcting for ME and LOD may help identify potentially useful discriminators of disease status that might have been overlooked if only analyzed naively and without regard to the negative impact of ME and LOD. Biomarker showing potential despite measurement difficulties may warrant future investment of resources to improve measurement procedures and technologies to realize this ability.