Multivariate normally distributed biomarkers subject to limits of detection and receiver operating characteristic curve inference

Endometriosis is a gynecologic disease exclusive to menstruating species such as humans and primates and occurs primarily in reproductive-age women. Data from experimental animal, primate, and human studies suggests an association between dioxin and PCBs and endometriosis ( 18 ). Here, we chose estrogenic congeners PCBs 188, 153, and 180 with LODs d 188 = 0.020, d 153 = 0.200, and d 180 = 0.034 to illustrate our method, by exploring these environmental toxicants as potential indicators of endometriosis using biomarker levels on 28 cases and 51 controls. These LODs of PCBs 188, 153, and 180 resulted in censoring of 32%, 64%, and 0% of cases and 36%, 74%, and 11% of controls, respectively. Considering each biomarker alone, the appropriateness of univariate normal distributions wase assessed via the use of QQ plots ( Fig 2 ) on the log-transformed biomarkers, similar to Lyles et al ( 8 ). Diagonal points corresponding to the levels above the LODs shows the three biomarkers fit log normal distributions well based on the quantifiable samples. The horizontal points are essentially quantile placeholders for the observations below the LOD. Assuming that the data we do not see follow the data we do see, univariate log normal likelihoods led to the MLEs AÛC 188 = 0.518, AÛC 153 = 0.511, and AÛC 180 = 0.630 ( 11 ). PCB 180 has a markedly higher estimated ability to discriminate and might be the choice if a single marker had to be chosen. Individuals are, however, exposed to these biomarkers as mixtures and are closely linked biologically, leading to them being highly correlated. Thus, using multiple biomarkers jointly makes sense here. Current practice is simply to sum the congeners, producing a composite exposure but one that is overly simplified ( 19 ). We will use developments from Materials and Methods to construct a more effective combination for the outcome endometriosis. Applying our developments from Materials and Methods to the vector of log-transformed PCBs 188, 153, and 180, we maximize Equation 1 after grouping the Z⃗ X i and Z⃗ Y j by pattern of censoring, Q⃗ ( Z⃗ ), to generate MLEs While these mean and covariance estimates are informative on their own, we substitute them into Equation 3 and get the BLC β → ^ 0 T = ( 0.164 , - 1.496 , 1.598 ) . We could now apply these coefficients and generate a parametric ROC curve for the BLC, which we compare to a similar curve for PCB 180 in Figure 3 . It is obvious in Figure 3 that the estimated BLC ROC curve (solid) is completely and substantially higher than that of estimated PCB 180 (dashed) and implies a potential increase in discrimination. We can quantify this using the AÛC 0 by simply using Equation 4 to see what might be achieved using all three biomarkers jointly for separating women with and without endometriosis. Combining PCBs 188, 153, and 180 in this manner yields AÛC 0 = 0.731, an increase in estimated discriminatory ability from any single biomarker. Using Equations A2 and A5 , we can also generate a 95% CI of 0.677 to 0.779. While these estimates show a potential benefit of using three biomarkers over one, ROC analysis is a process in which superiority of one biomarker over another for practical purposes needs to be validated with independent data before a recommendation can be made. If we use conventional naïve replacement methods and substitute the d for values below the LOD, we can also obtain a naïve BLC. The linear combination will certainly not be “best” in the sense of maximizing AUC because the assumption of multivariate normality is clearly violated. However, we generated such an estimate, AÛC 0 R = 0.652 for the sake of comparison, which is significantly less than our AÛC 0 = 0.731.

The MLEs developed here are asymptotically unbiased. The small sample performance of these estimators was evaluated via simulation in the R programming language and through recorded bias and root mean square error (RMSE) over B = 2000 iterations. Coverage probabilities and width of 95% CIs for AUC 0 were recorded. Increasing sample sizes, n x = n Y = 50, 100, 200, were used to demonstrate estimator performance. Nondiseased values were generated from a p = 3 multivariate normal distribution with mean zero and unit diagonal covariance matrices with uniform correlation between bio-markers of ρ = 0.5. Additional simulations were conducted with a higher uniform correlation between biomarkers of ρ = 0.8. Diseased values were similarly obtained but with a mean μ⃗ X corresponding to AUC 0 = 0.6, 0.7, 0.8. For each simulated group of samples, μ → ^ X , Σ̂ X , μ → ^ Y , and Σ̂ Y were calculated, leading to β → ^ 0 T , AÛC 0 , and 95% CI for AUC 0 . Due to the complexity of the elements of Cov ( μ → ^ , ∑ ^ ) in Equation 2 , occasional numeric integration was necessary within some of the elements or portions of the elements. We calculated these estimators while applying various levels of censoring, where d values corresponded to the quantiles of the marginal distributions of the nondiseased population. First “even” censoring was invoked by applying a similar d to all resulting biomarkers. Levels of even censoring were set at the 10th, 25th, and 50th quantiles across all biomarkers. To assess uneven censoring, d values were chosen at the following quantile sets: (10th, 10th, 50th), (10th, 25th, 50th), and (10th, 50th, 75th). As expected, the bias and RMSE of AÛC 0 decrease as sample sizes increase. The relative bias, ( AÛC 0 − AUC 0 )/ AUC 0 , ranged from 0.010 to 0.152, 0.003 to 0.088, and −0.0001 to 0.050 for n x = n y = 50, 100, and 200, respectively, including even the highest levels of censoring. Ranges of the RMSE were 0.0433 to 0.1443, 0.0305 to 0.0892, and 0.0213 to 0.0546 for n x = n y = 50, 100, and 200, respectively. The majority of the relative biases and RMSEs tended to the lower end of these ranges in all sample sizes. Figure 1 depicts the relative bias and RMSE for the ρ = 0.5 and AUC 0 = 0.6 (triangles) and 0.8 (squares) cases across all scenarios of censoring and as sample size increases n x = n y = 50, 100, and 200 (corresponding to increasing point size). The first column of points of each plot is for the MLE with no censoring and serves as a reference for subsequent MLEs based on data with censoring. Clearly within each column both relative bias and RMSE decrease as function of AUC 0 (same-size squares lower than triangles) and sample size (increasing size of same shape). Also evident is that relative bias and RMSE increase as censoring increases, generally from left to right. Figure 1 shows that relative bias and RMSE are only slightly increased when comparing columns 2, 3, 5, and 6, which have variations of 10% and 25% censoring, to column 1 based on the full data. Columns 4 and 7 show that while relative bias and RMSE are markedly increased from column 1, estimation is still potentially useful, accounting for 50% and 75% censored biomarker data. This is especially true for n x = n y = 100 and 200 near the middle and bottom, respectively, of both plots. While not displayed here, results for AUC 0 = 0.7 generally lie in between points in Figure 1 , and results for ρ = 0.8 show precisely the same relationships with a small increase in magnitude of relative bias and RMSE. The coverage probability of 95% CI’s were nominal or near nominal coverage for all AUC 0 averaging 0.940, 0.945, and 0.945 for n x = n y = 50, 100, and 200, respectively. No discernible pattern exists regarding ρ = 0.5 and 0.8 or with regard to censoring patterns due to the LOD.

As exposure assessment continues to stray from a single cause to mixtures of chemicals, images, and biomarkers, multiple bio-markers are increasingly being assessed and often measured simultaneously. However, whether measuring single or multiple biomarkers at once, emerging or established limitations of laboratory sensitivity resulting in operational limits, both upper and lower, are nothing new. The methods developed here provide a solution to the very common real-world data problem of censored values below the LOD, allowing for improved assessment of the biomarkers’ joint distribution and countless secondary analyses like regression, Bayesian inference, and, as we have shown, the ROC curve. While most imaging data are not left censored, continuous imaging variables may be part of a combination along with other biomarkers that are as censored. Multivariate normality is a widely used distribution and is often a reasonable assumption. Here, as in univariate literature, we rely on the assumption that the distribution of the data that we see above the LOD appropriately characterizes the data below that we do not see. This assumption allows for the application of maximum likelihood methods, whose estimator’s asymptotic properties, consistency and efficiency, are well known. We used a likelihood function that properly accounted for missing data due to LODs for P ≥ 2 biomarkers to generate MLEs for the mean vector and covariance matrix and further used it to develop the Fisher information for such parameters. The developments and simulation study for AÛC 0 demonstrated how these MLEs could be used to generate secondary estimators, which are also consistent and use the Fisher information to create accompanying CIs that achieve nominal coverage probability for a sample size as small as 50 with up to 50% missing values due to LODs. In our simulations, setting the true variances to be equal displayed the performance corresponding to the common setting where a mean shift occurs based on differing disease status. The behavior of AÛC 0 , bias, and variability, will differ based on each biomarker scenario and unequal variances, where our simulation might be conservative. As with all parametric estimation, the appropriateness of distributional assumptions is critical to the reliability of the estimators generated. The literature 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 ( 11 , 20 ). Similar caution and diligence should be taken when using the methods developed here for μ ⇀ ^ X , μ ⇀ ^ Y , Σ̂ X , Σ̂ Y , and, subsequently, AÛC 0 . In our example, we looked at marginal QQ plots to assess the assumption of normality of each PCB and, subsequently, assumed multivariate normality. Regarding the latter assumption, there are cases where multivariate normality fails to follow marginally normally distributed variables, but, practically speaking, these cases are rare and obvious. Care should always be taken to plot and investigate bivariate and multivariate data before any analyses. However, when inference that relies on the distribution of the data below the LOD is of interest, assumptions regarding the behavior of the biomarker values must be made. Using a replacement value makes an assumption that is almost assuredly false by creating a point mass at the chosen value a that grossly misrepresents the biomarkers distribution below the LOD. The likelihood function developed here can also be used in other applications where biomarkers are measured with LODs and normality can be assumed. Specifically, the ratio of likelihoods might be used to create an indices resulting in a proper ROC curve, without a hook such as in Figure 3 , which has been shown to be important in the assessment of medical imaging ( 21 , 22 ). Additionally, our accompanying covariance developments lead to nominal coverage of CIs. A simpler, alternative might be to use the Hessian generated during the maximization of the likelihood in lieu of our more complex findings. However, a thorough comparison of these two techniques should be conducted in future work before recommending a particular method. Future developments are also necessary to determine if the increased performance of combined biomarkers over that of a single biomarker, similar to that in our example, is statistically significant. ROC curves are of ever-increasing importance in reader studies of radiologic devices and clinical diagnostics in general ( 23 , 24 ). These methods for estimating μ⃑ X , μ⃑ Y , Σ X , Σ Y , and an ROC curve with AUC 0 while properly accounting for missing data censored below an LOD will provide a consistent view of the potential discriminatory ability of the BLC of a set of biomarkers. Moreover, the binormal ROC curve of the BLC indices means that we can combine methods here for multiple biomarkers with LODs with other important ROC analysis techniques ( 3 , 24 , 25 ). Assessing these biomarkers and diagnostic data properly will lead to better assessments of exposure and outcome relations. The hope is that we may find biomarkers that were once deemed useless displaying promise that warrants additional resources to improve the measurement process and subsequently lead to improved diagnostic care.

Suppose that biomarker levels W⃗ follow a multivariate distribution g ( W⃗ ; θ⃗ ) but measurements Z⃗ are actually observed with fixed LODs d⃗ = ( d 1 , …, d p ) T . Clearly, in most cases a likelihood function for parameters θ⃗ developed for the random sample W would be inappropriate for sample Z. This likelihood would not account for the censored values, even when censored values are replaced by a constant. A likelihood function can be formed that does account for censoring and maximizing the appropriate likelihood results in MLEs that are consistent, asymptotically unbiased. To build this likelihood, let Q⃗ ( Z⃗ ) be a p -dimensional vector of indicator functions, Q l = I ( Z l ≥ d l ), that take the value 1 when the measured bio-marker is greater than or equal to the LOD and 0 when it is less than the LOD. This creates, 2 p different Q⃗ ( Z⃗ ) that are possible and can be observed if each biomarker is measured with LOD. The extremes Q⃗ ( Z⃗ ) = (1, 1, 1, …, 1), indicating no censoring, and Q⃗ ( Z⃗ ) = (0, 0, 0, …, 0), indicating that all of that individual’s measurements are censored, set the edges of a variety of mixed observed and censored vectors. Individuals with the same pattern of censoring, Q⃗ ( Z⃗ ), have the same likelihood form. Thus, we can use these different patterns of censoring to construct a likelihood for our entire sample to maximize. The likelihood for the bivariate normal distribution in Lyles et al ( 8 ) had 2 2 = 4 parts for p = 2. Here we are considering the case where the biomarkers, W⃗ are multivariate normally distributed with mean vector μ⃗ and covariance matrix Σ. Suppose that we only observe the measured Z⃗ . Let us start creating our likelihood by first looking at the Q⃗ ( Z⃗ ) = (1, 1, 1, …, 1), where the likelihood function is the standard result L ( μ⃗ , Σ ; Z⃗ ) = L 1 ( μ⃗ , Σ ; w 1 , w 2 …, w p ) = g ( w⃗; μ⃗ , Σ ). The other extreme, where all biomarkers are censored, Q⃗ ( Z⃗ ) = (0, 0, 0, …, 0), corresponds to a likelihood, L ( μ⃗ , Σ ; Z⃗ ) = L 2 p ( μ⃗ , Σ ; w 1 < d 1 , w 2 <; d 2 , …, w p < d p ), that is the simple probability that all biomarkers are less than the LODs. While the details are presented in the Appendix , the probability that all of an individual’s biomarker levels are below the LODs for a given μ⃗ and Σ will be denoted by Λ ( d⃗ ; μ⃗ , Σ ). For the likelihoods of all the cases of Q⃗ ( Z⃗ ) in between, L k ( μ⃗ , Σ; Z⃗ ), k = 2,3, …,2 p − 1; that is, those Z⃗ composed of both measured and censored levels, we extend ideas from explicit results for p = 2 to the general p case here ( 8 , 15 ). In that work, the bivariate normal case with Q⃗ ( Z⃗ ) = (1, 0) or (0, 1) was shown to lead to a likelihood that is the product of a marginal distribution for the biomarker that is observed and a conditional distribution for the biomarker that is censored. We will rely on the same rationale here in a multivariate framework. All mixed vectors can be reordered and partitioned to have the same general form such that Z⃗ = ( Z⃗ (1) , Z⃗ (2) ), where the leading set is the observed biomarkers and the latter is the censored. The resultant multivariate normal density is then factored into a marginal density, g (1) ( W⃗ (1) ; μ⃗ (1) , Σ 11 ), and a conditional density, g (2|1) ( W⃗ (2) ; μ⃗ ( Q ) , Σ Q ), with μ → ( Q ) = μ → ( 2 ) + ∑ 21 ∑ 11 - 1 ( W → ( 1 ) - μ → ( 1 ) ) and ∑ Q = ∑ 22 - ∑ 21 ∑ 11 - 1 ∑ 12 . Since the observations in Z⃗ (2) are censored below the LODs, their contribution to the likelihood is simply the probability that those biomarkers would be below those LODs. Each of the k = 2, 3, …, 2 p − 1 likelihood contributions has the same form So for a random sample Z of independent identically distributed p -dimensional vectors Z⃗ j , j = 1, …, n , following a multivariate normal distribution, the likelihood has the form Maximizing Equation 1 with respect to each element of μ⃗ and Σ results in MLEs μ → ^ and Σ̂ . This must be performed numerically due to the complexity of the likelihood for which a closed-form solution does not exist for even the univariate case. The details of this development are found in the Appendix . To accompany these point estimates, we used Taylor expansion in the Appendix to generate covariance matrices that will allow us to create CIs for the AUC. Of note is that it is well known that μ → ^ and Σ̂ are independent in the setting without LODs. We have lost that property thanks to the introduction of the conditional probability into the likelihood when some of the data are below the LOD. Suppose the multiple biomarkers measured with LODs are intended to be used as diagnostic indices of disease status for a particular outcome. The ROC curve is used to assess such diagnostic ability by plotting sensitivity [the probability of a positive test result given the individual is truly diseased, denoted by q ( c )] by 1 – specificity [the probability of a negative test result given the individual is truly healthy, denoted by p ( c )] across all possible cut points c . Often a positive or negative test is indicated by having a biomarker level above or below, respectively, a given c but can be similarly used when the opposite is true. The ROC curve can be used to assess discriminatory ability over all c , over a specific range of q ( c ) or p ( c ), and to identify cut points that result in a specific discriminatory ability ( 3 , 4 ). The AUC, or P( X > Y ), is the most commonly used summary measure, ranging from 0.5 to 1, with larger values indicating a greater ability to discriminate ( 3 , 4 ). Consideration has already been given to its estimation based on a single biomarker affected by LOD ( 11 ). Su and Liu have shown that a linear combination can achieve discrimination greater than using a single biomarker alone ( 16 ). Let p biomarkers for individuals with and without a particular disease follow multivariate normal distributions, X⃑ ~ N p ( μ⃑ X , Σ X ) and Y⃗ ~ N p ( μ⃑ Y , Σ Y ), respectively. Then, indices formed from the linear combinations of those biomarkers, U = β → T X → = ∑ j = 1 p β j X j and V = β → T Y → = ∑ j = 1 p β j Y j , are univariate normally distributed, U ~ N 1 ( β⃗ T μ⃗ X , β⃗ T Σ X β⃗ ) and V ~ N 1 ( β⃗ T ⇀ μ Y , β⃗ T Σ Y β⃗ ), respectively. Consequentially, an ROC curve can be formed from U and V that depends on the choice of β⃗ T . The choice of β⃗ T can be equally weighted or an alternative designated by expert knowledge. In this setting, Su and Liu developed the BLC of multiple biomarkers (2) β → 0 T ∝ ( μ → X - μ → Y ) T ( ∑ X + ∑ Y ) - 1 , leading to (3) AUC 0 = P ( U > V ) = Φ ( ( μ ⇀ X - μ ⇀ Y ) T ( ∑ X + ∑ Y ) - 1 ( μ ⇀ X - μ ⇀ Y ) ) , which is “best” with respect to maximizing AUC over all real β⃗ T ( 16 ). Using samples measured with LOD, Z x , and Z Y , we can maximize the likelihood in Equation 1 for both groups, generating MLEs μ → ^ X and Σ̂ X and μ ⇀ ^ Y and Σ̂ Y . Placing these MLEs in Equations 2 and 3 results in MLEs β → ^ 0 T and AÛC 0 . Without loss of generality, we assume that both diseased and nondiseased biomarker levels are affected by the censoring at the same level, d xl = d Yl = d l . The development of the CI for AUC 0 is shown in the Appendix ( 17 ).