|Home | About | Journals | Submit | Contact Us | Français|
Multiple alternative diagnostic tests for one disease are commonly available to clinicians. It’s important to use all the good diagnostic predictors simultaneously to establish a new predictor with higher statistical utility. Under the generalized linear model for binary outcomes, the linear combination of multiple predictors in the link function is proved optimal in the sense that the area under the receiver operating characteristic (ROC) curve of this combination is the largest among all possible linear combination. The result was applied to analysis of the data from the Study of Osteoporotic Fractures (SOF) with comparison to Su and Liu’s approach.
Multiple alternative diagnostic tests for one disease are commonly available to clinicians, and statistical methods are available for evaluating these alternatives. For instance, multivariate regression models can be used to determine the marginal advantage of additional tests (Richards, 1995), and this approach has been used in osteoporosis research to determine the relative risk and/or odds ratio of osteoporotic frature (Hans, 1999). While these parameters are important for determining future fracture risk or developing clinical predictor rules based on tests, multivariate regression does not directly provide information on test utility or optimizing test combinations.
Because sensitivity and specificity of a diagnostic test depend on the threshold used to define abnormal, the receiver operating characteristic (ROC) curve is often used to assess the utility of a diagnostic test (DeLong, 1988; Pepe, 1997). This methodology has been extended to multivariate tests by constructing a linear function and then evaluating a combined utility based on the ROC of the linear function. Su and Liu (1993) have shown that under a normal distribution assumption, the best linear combination of diagnostic tests to achieve maximum area under a ROC curve is the linear discriminant function. However, no literature covers the general situation.
In this article, under the generalized linear model, we provide the optimal linear combination among all possibilities under the criterion that the area under the corresponding ROC curve of the combination is maximized, which is named the ROC criterion. The main results are presented in Section 2. Section 3 considers the estimation of the optimal linear combination and its corresponding area under ROC curve. Section 4 demonstrates an application of the proposed method to data from the Study of Osteoporotic Fractures (SOF) with comparison to Su and Liu’s approach. Discussion and conclusions are presented in the last section.
In this section it is proved that the linear combination in the link function of a model for binary responses maximizes sensitivity uniformly at any given specificity under the generalized linear model. Here we only present results for a simple case with two predictors, which can be easily generalized to the situation of multiple diagnostic markers.
Let Z be the binary outcome, (X1, X2) be the random predictive variables. Suppose the generalized linear model for binary responses holds, that is to say,
where β = (β0, β1, β2)′ is the 3-vector of parameters and h is a known function that is bounded in the unit interval (0,1). Although a wide choice of link function is available, the most commonly used in practical are the three models as follows:
Let α1X1 + α2X2 be any linear combination, and their sensitivity and specificity are respectively Sn and Sp. First, we will show that the selected linear combination by the link function of the model, β1X1 + β2X2, dominates all the other possible linear combinations in the sense that it provides the highest sensitivity uniformly over the entire range of specificity and thus, it is the best linear combination that results in the largest area under the ROC curve (AUC) among all linear combinations. Then we give the mathematical formula to compute the corresponding largest AUC.
The following useful lemma was originally due to Chebyshev (see Hardy, 1959)
Let U be a one-dimensional non-degenerate random variable. If g is bounded, positive, and strictly decreasing, then
holds provide the expectation of U exists.
Suppose the generalized linear model (1) for binary responses holds and the function h is continuous and strictly monotone. If (X1, X2) is two-dimensional continuous variable with continuous probability density and its expectation exists, then for any given specificity Sp, the coefficients for the best linear combination that provides the highest sensitivity uniformly is (α1, α2) (β1, β2).
Theorem 1 may be easily extended to multiple-markers situation, which is presented in the following theorem without proof.
Suppose the generalized linear model for binary responses holds, i.e.,
where β = (β0, β1, ···, βk)′ is the (k + 1)-vector of parameters with k ≥ 2. Let α1X1 + αkXk be any linear combination. If h is a known function that is also continuous, strictly monotone and bounded in the unit interval (0,1), (X1, ···, Xk) is k-dimensional continuous variable with continuous probability density and its expectation exists, then for any given specificity, the coefficients for the optimal linear combination that provides the highest sensitivity uniformly is (α1, ···, αk) (β1, ···, βk).
Suppose the conditions in Theorem 2 hold. Then the best linear combination under the ROC criterion is also (α1, ···, αk) (β1, ···, βk). And if h is strictly increasing, then the largest AUC is
where P(Z = 1) = ∫h(x)dG(x), P(Z = 0) = 1 − P(Z = 1) and G(x) is the distribution function of X = β0 + β1X1 + ··· + βkXk, which can be derived from the joint distribution of (X1, ···, Xk).
For rare diseases such as Osteoporotic Fractures discussed in the example section, we have the following simple formula to calculate the largest area under the ROC curve of the best linear combination, which can be directly derived from simple computation.
If the prevalence of the disease, P(Z = 1), is small enough, for example less than 5%, then the largest AUC may be approximated as follows:
whose error is about P(Z = 1)/2, say less than 2.5% if P(Z = 1) ≤ 5%.
In practical application, we need first to fit the generalized linear model from data. Suppose that for the study subject i, the response Zi has the Bernoulli distribution
and θi is determined by the link function θi = h(β0 + β1xi1 +···+ βkxik), where β = (β0, β1, ···, βk)′ is the (k + 1) -vector of parameters with k ≥ 2, h is a known continuous and strictly monotone function bounded in the unit interval (0,1), and (xi1, ···,xik) are the values of the two diagnostic predictors included in the model for the ith subject (i = 1, ···,n).
It is common to make use of the likelihood function as follows to estimate the parameters β = (β0, β1, ···, βk)′:
In fact, we can obtain the maximum likelihood estimate of β by maximizing L(β; z). An iterative algorithm, such as the popular Newton-Raphson method implemented in standard software from SAS and S-Plus, may be used to compute numerically (Cox, 1989).
As for estimation of the largest AUC, we need further know the joint distribution of the diagnostic predictors. There are two approaches for this problem. One belongs to conventional parametric methods, among which the normal approximation is the most popular. Another is the standard nonparametric technique, which use the empirical distribution, based on the sample data, to directly estimate the population distribution. After that, we could derive the estimation of the distribution function of β0 + β1X1 +···+ βkXk, and the corresponding area under the ROC curve based on Theorem 3 in this paper.
In this section, we apply our method to the data obtained from the Study of Osteoporotic Fractures (SOF) for illustration and compare it to the approach of Su and Liu (1993). From 1986 to 1988, SOF recruited 9,704 white women aged 65 years or older from four areas of the United States. At baseline, bone mineral density (BMD) was measured at the calcaneus, distal radius and proximal radius using single photon absorptiometry (SPA). At the second visit (1988–89), surviving participants had BMD measurements of the posterior-anterior (PA) spine (L1–L4) and proximal femur (neck, trochanter, total hip regions of interest) using dual x-ray absorptiometry (DXA). Fractures of the hip were recorded for each subject at each visit. More details about the study design and the data have been published previously (Cummings, 1995).
We included 7,127 women from the study. All of these women had forearm, calcaneal, hip and spine bone mineral density (BMD) measurements. In addition to the BMD measurements, many other previously identified predictive variables at baseline were also investigated. Furthermore, they all had known 5-year hip fracture status: either they were followed for 5 years after visit 2 without hip fracture or they had hip fracture within five years after visit 2. Women who were lost to follow-up within 5 years without known hip fracture were excluded.
Although all these 43 variables, including patient demographic data, clinical BMD measured by DXA and SXA, medical history, X-ray for prevalent fracture and vertebral heights, functional status, vision test results, and nutrition, etc, have been identified as significant predictors of hip fracture risk and they may reflect different aspects of osteoporosis and aging and may help in understanding the etiology of osteoporosis and hip fracture, only a few of them are necessary to identify subjects with elevated fracture risks. Standard approaches to the analysis of binary data such as logistic regression and probit model show that a linear combination of age, the femoral neck BMD and loss of height could best predict the hip fracture. Our previous study also suggested that these 3 variables could be used to build a non-inferior classification rule to the optimum recursive partitioning rule (Jin, 2004). So here we consider these three predictors and want to find out the best linear combination under the ROC criterion.
Assume the generalized linear model (2) holds for the SOF data, where Z1 = 1 stands for hip fracture, X1, X2 and X3 denote age, femoral neck BMD and loss of height respectively. If we further assume , the standard software S-plus or SAS provides the following fitted logistic regression model:
It follows from our theorems 2 and 3 that X(l) = 0.075X1 − 8.90X2 + 0.100X3 is the optimal linear combination in the sense that it provides the highest sensitivity uniformly over the entire range of specificity and also under the ROC criterion. Hence, the estimated best coefficients are proportional to .
We need to estimate the joint distribution of (X1, X2, X3) in order to get the corresponding largest area under the ROC curve. If we further assume that (X1, X2, X3)τ ~ N μ,Σ, where τ means the transpose of a vector, the standard normal estimation gives = (71.06,0.65,3.36)τ and
which leads to −3.89 + 0.075X1−8.90X2 + 0.10X3 ~ N −4.05,1.22. It follows from Theorem 3 that the largest AUC under the logistic regression model is 0.798. The 250 times bootstrap standard error is 0.049.
We may obtain similar result under the probit model, that is to say, h(x) = Φ(x). Here, the fitted model of hip fracture risk is
which suggests X(p) = 0.036X1 −4.09X2 + 0.048X3 be the best linear combination under the probit model. Hence, the estimated best coefficients are proportional to , and the corresponding largest AUC be 0.805 (with standard error 0.050). By the way, if we use the approximate formula in Corollary 1, these two largest AUCs are estimated to be 0.788 and 0.795, just 1% lower than the true values.
We also can consider applying Su and Liu’s method to the hip fracture data. Under their approach, the two conditional distributions are estimated as
The general method in Section 3 of their paper may give the best coefficients as follows, with the estimated largest AUC 0.760.
For comparison of fit goodness of the three models, we directly estimate the ROC curves corresponding to the three coefficients and empirically, i.e. based on the SOF data itself, and then calculate the areas under the curves. It’s not surprising that the real AUCs are all 0.805 because of little difference between the three optimal linear combinations resulted by different models. Therefore, the logistic regression and probit model fit the SOF data very well while Su and Liu’s method underestimates the area under the ROC curve of the best linear combination.
In medical applications, it is important to use all the good diagnostic predictors of a disease simultaneously to establish a new predictor with higher statistical utility. The linear combinations of multiple predictors are of particular interest to us. In this paper, we consider a common issue where the generalized linear models are utilized to fit the data with binary outcomes. Now that the linear combination identified by the link function is proved optimal under the ROC criterion, we just need to estimate the best linear combination using the standard procedure once the generalized linear models pass the statistical test. Therefore it is easy for clinicians to get the optimal linear combination of multiple diagnostic predictors under the condition of generalized linear models.
One referee of Su and Liu’s paper pointed out, “Using logistic regression to identify tests that best predict presence or absence of disease is also common”. Although the logistic regression is usually less efficient than the normal discriminant analysis with the normal assumption holding (Efron, 1975; Ruiz-Velasco, 1991), the generalized linear models for binary data will be more robust than the latter, because choice and estimation of the best linear combination need no assumption of the joint distribution of the multiple predictors. Su and Liu admitted this feature, exposed clearly by Cox and Snell (1989) that “once a vector of explanatory variables is given, then the probability that this individual belongs to one of the two groups is determined.” Our SOF data could be a real example for illustration of that point.
There may be two directions for extension of this paper. On the one hand, we can extend our work from linear to non-linear models, which may explore a new concept--the ROC region to assess statistical utility of a diagnostic predictor instead of the ROC curve. On the other hand, we may extend our interest from binary data to more complex data with more than two outcomes, which was discussed by Yang and Carlin (2000) who used ROC surface approach. All these kinds of research are currently under investigation.
The study is supported by grants from the National Institutes of Health R01EB004079 and National Bureau of Statistics of China LX:2006B45
Without loss of generality, let α2 = β2. For simplicity, let α1 = α.
Suppose that β2 > 0 and h is strictly increasing. Let f(x1, x2) be the distribution density function of the two-dimensional continuous variable (X1, X2), and (x) = 1 − h(x), then the specificity of the linear combination, αXl + β2X2, can be expressed as
For any given specificity Sp, the equation (A1) defines c as a differentiable function of α, denoted as c(α) below. Furthermore, it follows from differentiating both sides of equation (A1) with respective to α that
Evaluating (A2) at the specific choice α = β1, and noting that the nonzero function (β0 + β1x + c(β1) − β1x) does not depend on x and so can be taken out of the above integral, leads to the equality:
where c′(α) denotes the derivative of c(α) with respect to α.
As for the sensitivity of X, it can be rewritten as
On the other hand, it’s easy to see that, for any fixed α, is the value of the density function of the linear combination αX1 + β2X2 at point c(α). Let be the density function of a random variable U, and g(x) = (β0 + β1x+ c(α) − αx). So, when α < βl, g is bounded, positive, and strictly decreasing. From Lemma 1, we have
which leads to the following inequality:
From (A2), we see that
Similarly, we can prove that S′(α) > 0 when α > β1. Therefore, combining with (A5), we finish the proof that S(α) has the absolute minimum, hence Sn(c(α)) has the absolute maximum, at α = β1 when β2 > 0 and h is strictly increasing.
In a similar way, it’s easy to prove that the coefficients for the best linear combination that provides the highest sensitivity uniformly is (α1, α2) (β1, β2) when β2 < 0 or h is strictly decreasing. Thus the proof of the theorem is complete.
Once again, we let (x) = 1 − h(x). It follows from (2) that
So we have P(Z = 1) = ∫h(x)dG(x), and the specificity of the best linear combination can be expressed as
Then , and
On the other hand, its sensitivity is
Then the largest area under the ROC curve is given by
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.